Functions for Movement Programs

##############################
# PITCHf/x Movement Function #
##############################

PFX_Mvt <- function(PITCH,y_loc){
    # Set gravity
    g <- -32.174;
    # Find the time for the pitch to reach y_loc feet
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from 40 feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the movement in the x- and z-directions (in inches)
    pfx_x <- 6*PITCH$ax*t_diff^2;
    pfx_z <- 6*(PITCH$az-g)*t_diff^2;
    # Return the PITCHf/x movement as a vector
    return(c(pfx_x,pfx_z));
}

#######################################
# PITCHf/x Movement with Drag Removed #
#######################################

PFX_Mvt_Drag <- function(PITCH,y_loc){
    # Set gravity
    g <- -32.174;
    G <- c(0,0,g);
    # Find the time for the pitch to reach 40 feet
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from y_loc feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the velocity at 40 feet
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vy_loc <- PITCH$vy0 + PITCH$ay*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the velocity at the front of the plate
    vx_plate <- PITCH$vx0 + PITCH$ax*t_plate;
    vy_plate <- PITCH$vy0 + PITCH$ay*t_plate;
    vz_plate <- PITCH$vz0 + PITCH$az*t_plate;
    # Find the average velocity
    vx_avg <- (vx_loc + vx_plate)/2;
    vy_avg <- (vy_loc + vy_plate)/2;
    vz_avg <- (vz_loc + vz_plate)/2;
    v_avg <- sqrt(vx_avg^2 + vy_avg^2 + vz_avg^2);
    # Set the average velocity and acceleration vectors
    V <- c(vx_avg,vy_avg,vz_avg);
    A <- c(PITCH$ax,PITCH$ay,PITCH$az);
    # Normalize the velocity vector
    V_norm <- V/sqrt(V%*%V);
    # Find the drag vector
    Drag <- -abs((A-G)%*%V_norm)*V_norm;
    # Find the acceleration with drag removed
    ax_d <- PITCH$ax - Drag[1];
    az_d <- PITCH$az - Drag[3];
    # Find the movement without drag
    m_x <- 6*ax_d*t_diff^2;
    m_z <- 6*(az_d-g)*t_diff^2;
    # Return the movement without drag as a vector
    return(c(m_x,m_z));
}

##########################
# Planar Movement (in w) #
##########################

W_Mvt <- function(PITCH,y_loc){
    # Set gravity
    G <- c(0,0,-32.174);
    # Find the time for the pitch to reach y_loc feet in (x,y,z)-space
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate in (x,y,z)-space
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from 40 feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the PITCHf/x basis
    BV <- PFX_Basis(PITCH);
    # Transpose the matrix
    BV_T <- t(BV);
    # Find the w-acceleration of the pitch
    aw <- BV_T[3,] %*% c(PITCH$ax,PITCH$ay,PITCH$az);
    # Find the contribution of gravity in the w-direction
    gw <- BV_T[3,] %*% G;
    # Return the in-plane movement
    return((6*(aw-gw)*t_diff^2)*c(BV_T[3,1],BV_T[3,3]));
}

############################################
# Planar Movement (in w) with Drag Removed #
############################################

W_Mvt_Drag <- function(PITCH,y_loc){
    # Set gravity
    G <- c(0,0,-32.174);
    # Find the time for the pitch to reach y_loc feet in (x,y,z)-space
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate in (x,y,z)-space
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from y_loc feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the velocity at 40 feet
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vy_loc <- PITCH$vy0 + PITCH$ay*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the velocity at the front of the plate
    vx_plate <- PITCH$vx0 + PITCH$ax*t_plate;
    vy_plate <- PITCH$vy0 + PITCH$ay*t_plate;
    vz_plate <- PITCH$vz0 + PITCH$az*t_plate;
    # Find the average velocity
    vx_avg <- (vx_loc + vx_plate)/2;
    vy_avg <- (vy_loc + vy_plate)/2;
    vz_avg <- (vz_loc + vz_plate)/2;
    v_avg <- sqrt(vx_avg^2 + vy_avg^2 + vz_avg^2);
    # Set the average velocity and acceleration vectors
    V <- c(vx_avg,vy_avg,vz_avg);
    A <- c(PITCH$ax,PITCH$ay,PITCH$az);
    # Normalize the velocity vector
    V_norm <- V/sqrt(V%*%V);
    # Find the drag vector
    Drag <- -abs((A-G)%*%V_norm)*V_norm;
    # Find the PITCHf/x basis
    BV <- PFX_Basis(PITCH);
    # Transpose the matrix
    BV_T <- t(BV);
    # Find the w-acceleration of the pitch
    aw <- BV_T[3,] %*% c(PITCH$ax,PITCH$ay,PITCH$az);
    # Find the w-drag
    drag_w <- BV_T[3,] %*% Drag;
    # Find the contribution of gravity in the w-direction
    gw <- BV_T[3,] %*% G;
    # Find the w-acceleration without drag
    aw_d <- aw - drag_w;
    # Return the in-plane movement without drag
    return((6*(aw_d-gw)*t_diff^2)*c(BV_T[3,1],BV_T[3,3]));
}

##########################
# PITCHf/x Movement Norm #
##########################

PFX_Mvt_Norm <- function(PITCH,y_loc){
    # Set gravity
    g <- -32.174;
    # Find the time for the pitch to reach y_loc feet
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from 40 feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the movement in the x- and z-directions (in inches)
    pfx_x <- 6*PITCH$ax*t_diff^2;
    pfx_z <- 6*(PITCH$az-g)*t_diff^2;
    # Return the PITCHf/x movement as a vector
    return(sqrt(pfx_x^2 + pfx_z^2));
}

#############################################
# PITCHf/x Movement Norm with Drag Removed  #
#############################################

PFX_Mvt_Drag_Norm <- function(PITCH,y_loc){
    # Set gravity
    g <- -32.174;
    G <- c(0,0,g);
    # Find the time for the pitch to reach 40 feet
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from y_loc feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the velocity at 40 feet
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vy_loc <- PITCH$vy0 + PITCH$ay*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the velocity at the front of the plate
    vx_plate <- PITCH$vx0 + PITCH$ax*t_plate;
    vy_plate <- PITCH$vy0 + PITCH$ay*t_plate;
    vz_plate <- PITCH$vz0 + PITCH$az*t_plate;
    # Find the average velocity
    vx_avg <- (vx_loc + vx_plate)/2;
    vy_avg <- (vy_loc + vy_plate)/2;
    vz_avg <- (vz_loc + vz_plate)/2;
    v_avg <- sqrt(vx_avg^2 + vy_avg^2 + vz_avg^2);
    # Set the average velocity and acceleration vectors
    V <- c(vx_avg,vy_avg,vz_avg);
    A <- c(PITCH$ax,PITCH$ay,PITCH$az);
    # Normalize the velocity vector
    V_norm <- V/sqrt(V%*%V);
    # Find the drag vector
    Drag <- -abs((A-G)%*%V_norm)*V_norm;
    # Find the acceleration with drag removed
    ax_d <- PITCH$ax - Drag[1];
    az_d <- PITCH$az - Drag[3];
    # Find the movement without drag
    m_x <- 6*ax_d*t_diff^2;
    m_z <- 6*(az_d-g)*t_diff^2;
    # Return the movement without drag as a vector
    return(sqrt(m_x^2 + m_z^2));
}

###############################
# Planar Movement Norm (in w) #
###############################

W_Mvt_Norm <- function(PITCH,y_loc){
    # Set gravity
    G <- c(0,0,-32.174);
    # Find the time for the pitch to reach y_loc feet in (x,y,z)-space
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate in (x,y,z)-space
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from 40 feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the PITCHf/x basis
    BV <- PFX_Basis(PITCH);
    # Transpose the matrix
    BV_T <- t(BV);
    # Find the w-acceleration of the pitch
    aw <- BV_T[3,] %*% c(PITCH$ax,PITCH$ay,PITCH$az);
    # Find the contribution of gravity in the w-direction
    gw <- BV_T[3,] %*% G;
    # Return the in-plane movement
    return(6*(aw-gw)*t_diff^2);
}

#################################################
# Planar Movement Norm (in w) with Drag Removed #
#################################################

W_Mvt_Drag_Norm <- function(PITCH,y_loc){
    # Set gravity
    G <- c(0,0,-32.174);
    # Find the time for the pitch to reach y_loc feet in (x,y,z)-space
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_loc;
    t_loc <- (-b - sqrt(b^2 - 4*a*c_loc))/(2*a);
    # Find the time for the pitch to reach the front of the plate in (x,y,z)-space
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the time from y_loc feet to the front of the plate
    t_diff <- t_plate - t_loc;
    # Find the velocity at 40 feet
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vy_loc <- PITCH$vy0 + PITCH$ay*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the velocity at the front of the plate
    vx_plate <- PITCH$vx0 + PITCH$ax*t_plate;
    vy_plate <- PITCH$vy0 + PITCH$ay*t_plate;
    vz_plate <- PITCH$vz0 + PITCH$az*t_plate;
    # Find the average velocity
    vx_avg <- (vx_loc + vx_plate)/2;
    vy_avg <- (vy_loc + vy_plate)/2;
    vz_avg <- (vz_loc + vz_plate)/2;
    v_avg <- sqrt(vx_avg^2 + vy_avg^2 + vz_avg^2);
    # Set the average velocity and acceleration vectors
    V <- c(vx_avg,vy_avg,vz_avg);
    A <- c(PITCH$ax,PITCH$ay,PITCH$az);
    # Normalize the velocity vector
    V_norm <- V/sqrt(V%*%V);
    # Find the drag vector
    Drag <- -abs((A-G)%*%V_norm)*V_norm;
    # Find the PITCHf/x basis
    BV <- PFX_Basis(PITCH);
    # Transpose the matrix
    BV_T <- t(BV);
    # Find the w-acceleration of the pitch
    aw <- BV_T[3,] %*% c(PITCH$ax,PITCH$ay,PITCH$az);
    # Find the w-drag
    drag_w <- BV_T[3,] %*% Drag;
    # Find the contribution of gravity in the w-direction
    gw <- BV_T[3,] %*% G;
    # Find the w-acceleration without drag
    aw_d <- aw - drag_w;
    # Return the in-plane movement without drag
    return(6*(aw_d-gw)*t_diff^2);
}

###################
# Binormal Vector #
###################

Binormal_Vector <- function(PITCH){
    # Set the velocity vector
    V <- c(PITCH$vx0,PITCH$vy0,PITCH$vz0);
    # Set the acceleration vector
    A <- c(PITCH$ax,PITCH$ay,PITCH$az);
    # Find the cross product of velocity and acceleration
    V_x_A <- c(V[2]*A[3] - V[3]*A[2],V[3]*A[1] - V[1]*A[3],V[1]*A[2] - V[2]*A[1]);
    # Find the magnitude of the cross product
    V_x_A_norm <- sqrt(V_x_A[1]^2 + V_x_A[2]^2 + V_x_A[3]^2);
    # Normalize the cross product to get the binormal vector
    return(V_x_A/V_x_A_norm);
}

##################
# PITCHf/x Basis #
##################

PFX_Basis <- function(PITCH){
    # Find the binormal vector
    V1 <- Binormal_Vector(PITCH);
    # Find a unit vector in the general direction toward home plate
    V2 <- sign(V1[1])*c(-V1[2],V1[1],0)/sqrt(V1[1]^2 + V1[2]^2);
    # Find a unit vector perpendicular to both previous vectors
    V3 <- sign(V1[1]*V2[2]-V1[2]*V2[1])*c(-V1[3]*V2[2],V1[3]*V2[1],V1[1]*V2[2]-V1[2]*V2[1]);
    V3 <- V3/sqrt(V3[1]^2 + V3[2]^2 + V3[3]^2);
    # Store the vectors in an array
    V <- array(0,dim=c(3,3));
    V[,1] <- V1;
    V[,2] <- V2;
    V[,3] <- V3;
    return(V);
}

###################################################
# Minimum Distance from a Point to a Line Segment #
###################################################

Dist_Seg <- function(P,P0,P1){
    # Form a vector between the endpoints of the segments
    v = P1 - P0;
    # Form another vector between the point and one end of the segment
    w = P - P0;
    # Find the dot product between the two vectors
    c1 = v[1]*w[1] + v[2]*w[2];
    # Find the squared length of the line segment
    c2 = v[1]*v[1] + v[2]*v[2];
    # If c1 is less than or equal to zero, then the closest point on the line is the first endpoint
    if ( c1 <= 0 ){
    # Find the distance between the points
        d <- Dist_Min(P,P0);
        return(d);
    # If c2 is less than or equal to c1, then the closest point on the line is the second endpoint
    } else if( c2 <= c1 ){
    # Find the distance between the points
        d <- Dist_Min(P,P1);
        return(d);
    # Else, the closest point on the line is in the middle
    } else {
        b = c1/c2;
        Pb = P0 + b*v;
    # Find the distance between the points
        d <- Dist_Min(P,Pb);
        return(d);
      }
}

##################################################
# Location of a Pitch Relative to a Line Segment #
##################################################

Loc_Seg <- function(P,P0,P1){
    # Form a vector between the endpoints of the segments
    v = P1 - P0;
    # Form another vector between the point and one end of the segment
    w = P - P0;
    # Find the dot product between the two vectors
    c1 = v[1]*w[1] + v[2]*w[2];
    # Find the squared length of the line segment
    c2 = v[1]*v[1] + v[2]*v[2];
    # If c1 is less than or equal to zero, then the closest point on the line is the first endpoint
    if ( c1 <= 0 ){
    # Find the distance between the points
        d <- Dist_Min(P,P0);
        return(list(d=d,Pc=P0));
    # If c2 is less than or equal to c1, then the closest point on the line is the second endpoint
    } else if( c2 <= c1 ){
    # Find the distance between the points
        d <- Dist_Min(P,P1);
        return(list(d=d,Pc=P1));
    # Else, the closest point on the line is in the middle
    } else {
        b = c1/c2;
        Pb = P0 + b*v;
    # Find the distance between the points
        d <- Dist_Min(P,Pb);
        return(list(d=d,Pc=Pb));
      }
}

###############################
# Distance Between Two Points #
###############################

Dist_Min <- function(P0,P1){
    d <- sqrt((P1[1]-P0[1])^2 + (P1[2]-P0[2])^2);
    return(d);
}

#########################################################################
# Minimum Distance Between and Point and the Outside of a Quadrilateral #
#########################################################################

Min_Dist_Quad <- function(P,UL,UR,LL,LR){
    # Distance from the top of the quadrilateral
    d_top <- Dist_Seg(P,UL,UR);
    m_top <- (UR[2]-UL[2])/(UR[1]-UL[1]);
    b_top <- UR[2] - m_top*UR[1];
    y_top <- m_top*P[1] + b_top;
    sign_top <- sign(P[2]-y_top);
    # Distance from the bottom of the quadrilateral
    d_bottom <- Dist_Seg(P,LL,LR);
    m_bottom <- (LR[2]-LL[2])/(LR[1]-LL[1]);
    b_bottom <- LR[2] - m_bottom*LR[1];
    y_bottom <- m_bottom*P[1] + b_bottom;
    sign_bottom <- sign(y_bottom-P[2]);
    # Distance from the left of the quadrilateral
    d_left <- Dist_Seg(P,UL,LL);
    m_left <- (UL[1]-LL[1])/(UL[2]-LL[2]);
    b_left <- UL[1] - m_left*UL[2];
    x_left <- m_left*P[2] + b_left;
    sign_left <- sign(x_left-P[1]);
    # Distance from the right of the quadrilateral
    d_right <- Dist_Seg(P,UR,LR);
    m_right <- (UR[1]-LR[1])/(UR[2]-LR[2]);
    b_right <- UR[1] - m_right*UR[2];
    x_right <- m_right*P[2] + b_right;
    sign_right <- sign(P[1]-x_right);
    # Find the minimum distance from the outside of the quadrilateral (Inside = 0)
    if ( (sign_top < 0) && (sign_bottom < 0) && (sign_left < 0) && (sign_right < 0) ){
        d_min <- 0;
    } else {
        d_min <- min(d_top,d_bottom,d_left,d_right);
    }
    return(d_min);
}

################################################################
# Direction Between a Point and the Outside of a Quadrilateral #
################################################################

Dir_Quad <- function(P,UL,UR,LL,LR){
    # Distance from the top of the quadrilateral
    P_top <- Loc_Seg(P,UL,UR);
    m_top <- (UR[2]-UL[2])/(UR[1]-UL[1]);
    b_top <- UR[2] - m_top*UR[1];
    y_top <- m_top*P[1] + b_top;
    sign_top <- sign(P[2]-y_top);
    P_dir <- P - P_top$Pc;
    P_dist <- P_top$d;
    # Distance from the bottom of the quadrilateral
    P_bottom <- Loc_Seg(P,LL,LR);
    m_bottom <- (LR[2]-LL[2])/(LR[1]-LL[1]);
    b_bottom <- LR[2] - m_bottom*LR[1];
    y_bottom <- m_bottom*P[1] + b_bottom;
    sign_bottom <- sign(y_bottom-P[2]);
    if (P_dist > P_bottom$d){
        P_dir <- P - P_bottom$Pc;
    }
    # Distance from the left of the quadrilateral
    P_left <- Loc_Seg(P,UL,LL);
    m_left <- (UL[1]-LL[1])/(UL[2]-LL[2]);
    b_left <- UL[1] - m_left*UL[2];
    x_left <- m_left*P[2] + b_left;
    sign_left <- sign(x_left-P[1]);
    if (P_dist > P_left$d){
        P_dir <- P - P_left$Pc;
    }
    # Distance from the right of the quadrilateral
    P_right <- Loc_Seg(P,UR,LR);
    m_right <- (UR[1]-LR[1])/(UR[2]-LR[2]);
    b_right <- UR[1] - m_right*UR[2];
    x_right <- m_right*P[2] + b_right;
    sign_right <- sign(P[1]-x_right);
    if (P_dist > P_right$d){
        P_dir <- P - P_right$Pc;
    }
    # Find the minimum distance from the outside of the quadrilateral (Inside = 0)
    if ( (sign_top < 0) && (sign_bottom < 0) && (sign_left < 0) && (sign_right < 0) ){
        P_dir <- c(0,0);
    }
    return(P_dir);
}

######################################
# Find the Zone Containing the Pitch #
######################################

Zone_Quad <- function(P,UL,UR,LL,LR){
    # Distance from the top of the quadrilateral
    d_top <- Dist_Seg(P,UL,UR);
    m_top <- (UR[2]-UL[2])/(UR[1]-UL[1]);
    b_top <- UR[2] - m_top*UR[1];
    y_top <- m_top*P[1] + b_top;
    sign_top <- sign(P[2]-y_top);
    # Distance from the bottom of the quadrilateral
    d_bottom <- Dist_Seg(P,LL,LR);
    m_bottom <- (LR[2]-LL[2])/(LR[1]-LL[1]);
    b_bottom <- LR[2] - m_bottom*LR[1];
    y_bottom <- m_bottom*P[1] + b_bottom;
    sign_bottom <- sign(y_bottom-P[2]);
    # Distance from the left of the quadrilateral
    d_left <- Dist_Seg(P,UL,LL);
    m_left <- (UL[1]-LL[1])/(UL[2]-LL[2]);
    b_left <- UL[1] - m_left*UL[2];
    x_left <- m_left*P[2] + b_left;
    sign_left <- sign(x_left-P[1]);
    # Distance from the right of the quadrilateral
    d_right <- Dist_Seg(P,UR,LR);
    m_right <- (UR[1]-LR[1])/(UR[2]-LR[2]);
    b_right <- UR[1] - m_right*UR[2];
    x_right <- m_right*P[2] + b_right;
    sign_right <- sign(P[1]-x_right);
    # Find the minimum distance from the outside of the quadrilateral (Inside = 0)
    if ( (sign_top < 0) && (sign_bottom < 0) && (sign_left < 0) && (sign_right < 0) ){
        d_min <- 0;
    } else {
        d_min <- min(d_top,d_bottom,d_left,d_right);
    }
    if (d_min <= (log(2)/4)^(1/4)){
        Zone <- 5;
        return(Zone)
    } else {
        if (sign_top > 0){
            Zone <- 0;
        } else if (sign_bottom > 0){
            Zone <- 6;
        } else {
            Zone <- 3;
        }
        if (sign_left > 0){
            Zone <- Zone + 1;
        } else if (sign_right > 0){
            Zone <- Zone + 3;
        } else {
            Zone <- Zone + 2;
        }
        return(Zone);
    }
}

##########################################
# Angle (in Degrees) Between Two Vectors #
##########################################

Dot_Prod_Angle <- function(V1,V2){
    # Find the dot product of the two vectors
    V1dotV2 <- V1 %*% V2;
    # Find the length of the two vectors
    V1_norm <- sqrt(V1 %*% V1);
    V2_norm <- sqrt(V2 %*% V2);
    # Return the angle, in degrees, between the two vectors
    return((180/pi)*acos(V1dotV2/(V1_norm*V2_norm)));
}

###################################################
# Strike Zone Probability Based on Location (LHB) #
###################################################

Dist_LHB_Quad <- function(P){
    # Define the four corners of the strike zone core (2016)
    UL <- c(-0.47,2.85);
    UR <- c(0.22,2.87);
    LL <- c(-0.5,2.04);
    LR <- c(0.24,2.16);
    # Find the distance from the pitch location to the core
    xz_dist <- Min_Dist_Quad(P,UL,UR,LL,LR);
    # Find the strike probability
    strike_prob <- exp(-4*xz_dist^4);
    return(strike_prob);
}

###################################################
# Strike Zone Probability Based on Location (RHB) #
###################################################

Dist_RHB_Quad <- function(P){
    # Define the four corners of the strike zone core (2016)
    UL <- c(-0.38,2.88);
    UR <- c(0.35,2.81);
    LL <- c(-0.46,2.1);
    LR <- c(0.35,2.09);
    # Find the distance from the pitch location to the core
    xz_dist <- Min_Dist_Quad(P,UL,UR,LL,LR);
    # Find the strike probability
    strike_prob <- exp(-4*xz_dist^4);
    return(strike_prob);
}

###################################################################################
# Track the projected strike probability with remaining PITCHf/x movement removed #
###################################################################################

Pitch_Track_PFX <- function(PITCHES,Ny,stand){
    g <- -32.174;
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    Proj_Prb <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Set a list for the 9 PITCHf/x parameters
        PITCH <- list(ax = ax[i], ay = ay[i], az = az[i], vx0 = vx0[i], vy0 = vy0[i], vz0 = vz0[i], x0 = x0[i], y0 = y0[i], z0 = z0[i]);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            p_xz <- PFX_Calc(PITCH,t_inc,t_plate);
            # Calculate the probability that the projected pitch is a strike
            if (stand == 'L'){
                strike_prob <- Dist_LHB_Quad(p_xz);
            }
            else {
                strike_prob <- Dist_RHB_Quad(p_xz);
            }
            Proj_Prb[i,Ny-j+1] <- strike_prob;
        }
    }
    return(Proj_Prb);
}

#####################################################################################################
# Track the projected strike probability with remaining PITCHf/x movement removed but drag added in #
#####################################################################################################

Pitch_Track_PFX_Drag <- function(PITCHES,Ny,stand){
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    Proj_Prb <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the velocity at the front of the plate
        vx_plate <- ax[i]*t_plate + vx0[i];
        vy_plate <- ay[i]*t_plate + vy0[i];
        vz_plate <- az[i]*t_plate + vz0[i];
        v_plate <- list(x=vx_plate, y=vy_plate, z=vz_plate);
        # Set a list for the 9 PITCHf/x parameters
        PITCH <- list(ax = ax[i], ay = ay[i], az = az[i], vx0 = vx0[i], vy0 = vy0[i], vz0 = vz0[i], x0 = x0[i], y0 = y0[i], z0 = z0[i]);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            p_xz <- PFX_Drag_Calc(PITCH,t_inc,t_plate,v_plate);
            # Calcuate the probability that the projected pitch is a strike
            if (stand == 'L'){
                strike_prob <- Dist_LHB_Quad(p_xz);
            }
            else {
                strike_prob <- Dist_RHB_Quad(p_xz);
            }
            Proj_Prb[i,Ny-j+1] <- strike_prob;
        }
    }
    return(Proj_Prb);
}

###################################################################################
# Track the planar projected strike probability with remaining w-movement removed #
###################################################################################

Pitch_Track_W <- function(PITCHES,Ny,stand){
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    Proj_Prb <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the x, y, and z coordinates of the pitch at the front of the plate
        x_plate <- x0[i] + vx0[i]*t_plate + 0.5*ax[i]*t_plate^2;
        y_plate <- y0[i] + vy0[i]*t_plate + 0.5*ay[i]*t_plate^2;
        z_plate <- z0[i] + vz0[i]*t_plate + 0.5*az[i]*t_plate^2;
        # Set up a list containing the parameters of the pitch
        PITCH <- list(ax=ax[i], ay=ay[i], az=az[i], vx0=vx0[i], vy0=vy0[i], vz0=vz0[i], x0=x0[i], y0=y0[i], z0=z0[i]);
        BV <- PFX_Basis(PITCH);
        BV_T = t(BV);
        # Find the u and b coordinates of the pitch at the front of the plate
        b_plate <- BV_T[1,] %*% c(x_plate,y_plate,z_plate);
        u_plate <- BV_T[2,] %*% c(x_plate,y_plate,z_plate);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            p_xz <- W_Calc(PITCH,BV_T,t_inc,t_plate,b_plate,u_plate);
            if (stand == 'L'){
                strike_prob <- Dist_LHB_Quad(p_xz);
            }
            else {
                strike_prob <- Dist_RHB_Quad(p_xz);
            }
            Proj_Prb[i,Ny-j+1] <- strike_prob;
        }
    }
    return(Proj_Prb);
}

##########################################################################################################
# Track the planar projected strike probability with remaining planar movement removed but drag added in #
##########################################################################################################

Pitch_Track_W_Drag <- function(PITCHES,Ny,stand){
    g <- -32.174;
    G <- c(0,0,g);
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    Proj_Prb <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the x, y, and z coordinates of the pitch at the front of the plate
        x_plate <- x0[i] + vx0[i]*t_plate + 0.5*ax[i]*t_plate^2;
        y_plate <- y0[i] + vy0[i]*t_plate + 0.5*ay[i]*t_plate^2;
        z_plate <- z0[i] + vz0[i]*t_plate + 0.5*az[i]*t_plate^2;
        plate_xz <- c(x_plate,z_plate);
        # Find the velocity at the front of the plate
        vx_plate <- ax[i]*t_plate + vx0[i];
        vy_plate <- ay[i]*t_plate + vy0[i];
        vz_plate <- az[i]*t_plate + vz0[i];
        v_plate <- list(x=vx_plate, y=vy_plate, z=vz_plate);
        # Set up a list containing the parameters of the pitch
        PITCH <- list(ax=ax[i], ay=ay[i], az=az[i], vx0=vx0[i], vy0=vy0[i], vz0=vz0[i], x0=x0[i], y0=y0[i], z0=z0[i]);
        BV <- PFX_Basis(PITCH);
        BV_T = t(BV);
        # Find the u and b coordinates of the pitch at the front of the plate
        b_plate <- BV_T[1,] %*% c(x_plate,y_plate,z_plate);
        u_plate <- BV_T[2,] %*% c(x_plate,y_plate,z_plate);
        # Set the acceleration vector
        A <- c(ax[i],ay[i],az[i]);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            p_xz <- W_Drag_Calc(PITCH,BV_T,t_inc,t_plate,v_plate,b_plate,u_plate);
            if (stand == 'L'){
                strike_prob <- Dist_LHB_Quad(p_xz);
            }
            else {
                strike_prob <- Dist_RHB_Quad(p_xz);
            }
            Proj_Prb[i,Ny-j+1] <- strike_prob;
        }
    }
    return(Proj_Prb);
}

###############################################################
# Determine a ball or strike based on strike zone probability #
###############################################################

ball_strike_prob <- function(PITCHES,stand){
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    Str_Prb <- array(0,dim=Np);
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the x and z coordinates of the pitch at the front of the plate
        x_plate <- x0[i] + vx0[i]*t_plate + 0.5*ax[i]*t_plate^2;
        z_plate <- z0[i] + vz0[i]*t_plate + 0.5*az[i]*t_plate^2;
        plate_xz <- c(x_plate,z_plate);
        if (stand == 'L'){
            Str_Prb[i] <- Dist_LHB_Quad(plate_xz);
        } else {
            Str_Prb[i] <- Dist_RHB_Quad(plate_xz);
        }
        Str_Index <- which(Str_Prb >= 0.5);
        Ball_Index <- which(Str_Prb < 0.5);
        }
    return(list(S = Str_Index, B = Ball_Index, Prb = Str_Prb));
}

##################################################################
# Determine to xz-location of a pitch at the front of home plate #
##################################################################

xz_plate_loc <- function(PITCH){
    x0 <- PITCH$x0;
    y0 <- PITCH$y0;
    z0 <- PITCH$z0;
    vx0 <- PITCH$vx0;
    vy0 <- PITCH$vy0;
    vz0 <- PITCH$vz0;
    ax <- PITCH$ax;
    ay <- PITCH$ay;
    az <- PITCH$az;
    # Find the time to the front of home plate
    a <- 0.5*ay;
    b <- vy0;
    c_plate <- y0-(17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the x and z coordinates of the pitch at the front of the plate
    x_plate <- x0 + vx0*t_plate + 0.5*ax*t_plate^2;
    z_plate <- z0 + vz0*t_plate + 0.5*az*t_plate^2;
    return(c(x_plate,z_plate));
}

#########################################################
# Determine the mean angle from a collection of vectors #
#########################################################

mean_angle <- function(x,z,N){
    x_norm <- array(0,dim=N);
    z_norm <- array(0,dim=N);
    # Find the length of each vector and normalize
    for (i in 1:N){
        xz_norm <- sqrt(x[i]^2 + z[i]^2);
        x_norm[i] <- x[i]/xz_norm;
        z_norm[i] <- z[i]/xz_norm;
    }
    # Find the mean of each component
    x_mean <- mean(x_norm);
    z_mean <- mean(z_norm);
    # Return the mean angle
    theta <- atan2(z_mean,x_mean);
    return(theta);
}

##############################################################################################
# Display the perctange of angles in groups of 30 degrees and distances from the strike zone #
##############################################################################################

Angle_Pct_Dist <- function(Pitch_Angle,SZ_Distance){
    # Find the number of pitches and indices for various ranges of angles from 0 to 180 degrees
    N <- length(Pitch_Angle);
    I_30 <- which(Pitch_Angle <= 30);
    N_30 <- length(I_30);
    I_60 <- which((Pitch_Angle > 30) & (Pitch_Angle <= 60));
    N_60 <- length(I_60);
    I_90 <- which((Pitch_Angle > 60) & (Pitch_Angle <= 90));
    N_90 <- length(I_90);
    I_120 <- which((Pitch_Angle > 90) & (Pitch_Angle <= 120));
    N_120 <- length(I_120);
    I_150 <- which((Pitch_Angle > 120) & (Pitch_Angle <= 150));
    N_150 <- length(I_150);
    I_180 <- which((Pitch_Angle > 150) & (Pitch_Angle <= 180));
    N_180 <- length(I_180);
    # Find the percentage of pitches in each range of angles
    Pcts <- round(c(N_30/N,N_60/N,N_90/N,N_120/N,N_150/N,N_180/N)*100,digits=1);
    # Find the average outside distance in each range of angles
    Avg_Dist <- round(c(mean(SZ_Distance[I_30]),mean(SZ_Distance[I_60]),mean(SZ_Distance[I_90]),mean(SZ_Distance[I_120]),mean(SZ_Distance[I_150]),mean(SZ_Distance[I_180])),digits=2);
    # Print the results
    print(paste("0-30: ",Pcts[1],"% / 30-60: ",Pcts[2],"% / 60-90: ",Pcts[3],"%",sep=""));
    print(paste("90-120: ",Pcts[4],"% / 120-150: ",Pcts[5],"% / 150-180: ",Pcts[6],"%",sep=""));
    print(paste("0-30: ",Avg_Dist[1],"ft / 30-60: ",Avg_Dist[2],"ft / 60-90: ",Avg_Dist[3],"ft",sep=""));
    print(paste("90-120: ",Avg_Dist[4],"ft / 120-150: ",Avg_Dist[5],"ft / 150-180: ",Avg_Dist[6],"ft",sep=""));
}

###########################################################################################
# Find the perctange of angles in groups of 30 degrees and distances from the strike zone #
###########################################################################################

Angle_Pct_Dist_Store <- function(Pitch_Angle,SZ_Distance){
    # Find the number of pitches and indices for various ranges of angles from 0 to 180 degrees
    N <- length(Pitch_Angle);
    I_30 <- which(Pitch_Angle <= 30);
    N_30 <- length(I_30);
    I_60 <- which((Pitch_Angle > 30) & (Pitch_Angle <= 60));
    N_60 <- length(I_60);
    I_90 <- which((Pitch_Angle > 60) & (Pitch_Angle <= 90));
    N_90 <- length(I_90);
    I_120 <- which((Pitch_Angle > 90) & (Pitch_Angle <= 120));
    N_120 <- length(I_120);
    I_150 <- which((Pitch_Angle > 120) & (Pitch_Angle <= 150));
    N_150 <- length(I_150);
    I_180 <- which((Pitch_Angle > 150) & (Pitch_Angle <= 180));
    N_180 <- length(I_180);
    # Find the percentage of pitches in each range of angles
    Pcts <- round(c(N_30/N,N_60/N,N_90/N,N_120/N,N_150/N,N_180/N)*100,digits=1);
    # Find the average outside distance in each range of angles
    Avg_Dist <- array(0,dim=6);
    if (N_30 != 0){
        Avg_Dist[1] <- round(mean(SZ_Distance[I_30]),digits=2);
    }
    if (N_60 != 0){
        Avg_Dist[2] <- round(mean(SZ_Distance[I_60]),digits=2);
    }
    if (N_90 != 0){
        Avg_Dist[3] <- round(mean(SZ_Distance[I_90]),digits=2);
    }
    if (N_120 != 0){
        Avg_Dist[4] <- round(mean(SZ_Distance[I_120]),digits=2);
    }
    if (N_150 != 0){
        Avg_Dist[5] <- round(mean(SZ_Distance[I_150]),digits=2);
    }
    if (N_180 != 0){
        Avg_Dist[6] <- round(mean(SZ_Distance[I_180]),digits=2);
    }
    # Store and return the data
    return(c(Pcts,Avg_Dist));
}

########################################################
# Find the perctange of angles in groups of 30 degrees #
########################################################

Angle_Pct_By_Case <- function(Pitch_Angle_C,Pitch_Angle_M,Pitch_Angle_T){
    # Find the number of pitches with contact and indices for various ranges of angles from 0 to 180 degrees
    IC_30 <- which(Pitch_Angle_C <= 30);
    C_30 <- length(IC_30);
    IC_60 <- which((Pitch_Angle_C > 30) & (Pitch_Angle_C <= 60));
    C_60 <- length(IC_60);
    IC_90 <- which((Pitch_Angle_C > 60) & (Pitch_Angle_C <= 90));
    C_90 <- length(IC_90);
    IC_120 <- which((Pitch_Angle_C > 90) & (Pitch_Angle_C <= 120));
    C_120 <- length(IC_120);
    IC_150 <- which((Pitch_Angle_C > 120) & (Pitch_Angle_C <= 150));
    C_150 <- length(IC_150);
    IC_180 <- which((Pitch_Angle_C > 150) & (Pitch_Angle_C <= 180));
    C_180 <- length(IC_180);
    # Find the number of pitches missed and indices for various ranges of angles from 0 to 180 degrees
    IM_30 <- which(Pitch_Angle_M <= 30);
    M_30 <- length(IM_30);
    IM_60 <- which((Pitch_Angle_M > 30) & (Pitch_Angle_M <= 60));
    M_60 <- length(IM_60);
    IM_90 <- which((Pitch_Angle_M > 60) & (Pitch_Angle_M <= 90));
    M_90 <- length(IM_90);
    IM_120 <- which((Pitch_Angle_M > 90) & (Pitch_Angle_M <= 120));
    M_120 <- length(IM_120);
    IM_150 <- which((Pitch_Angle_M > 120) & (Pitch_Angle_M <= 150));
    M_150 <- length(IM_150);
    IM_180 <- which((Pitch_Angle_M > 150) & (Pitch_Angle_M <= 180));
    M_180 <- length(IM_180);
    # Find the number of pitches taken and indices for various ranges of angles from 0 to 180 degrees
    IT_30 <- which(Pitch_Angle_T <= 30);
    T_30 <- length(IT_30);
    IT_60 <- which((Pitch_Angle_T > 30) & (Pitch_Angle_T <= 60));
    T_60 <- length(IT_60);
    IT_90 <- which((Pitch_Angle_T > 60) & (Pitch_Angle_T <= 90));
    T_90 <- length(IT_90);
    IT_120 <- which((Pitch_Angle_T > 90) & (Pitch_Angle_T <= 120));
    T_120 <- length(IT_120);
    IT_150 <- which((Pitch_Angle_T > 120) & (Pitch_Angle_T <= 150));
    T_150 <- length(IT_150);
    IT_180 <- which((Pitch_Angle_T > 150) & (Pitch_Angle_T <= 180));
    T_180 <- length(IT_180);
    # Sum the counts
    N_30 <- C_30 + M_30 + T_30;
    N_60 <- C_60 + M_60 + T_60;
    N_90 <- C_90 + M_90 + T_90;
    N_120 <- C_120 + M_120 + T_120;
    N_150 <- C_150 + M_150 + T_150;
    N_180 <- C_180 + M_180 + T_180;
    # Find the percentage of pitches in each range of angles
    C_Pcts <- round(c(C_30/N_30,C_60/N_60,C_90/N_90,C_120/N_120,C_150/N_150,C_180/N_180)*100,digits=1);
    M_Pcts <- round(c(M_30/N_30,M_60/N_60,M_90/N_90,M_120/N_120,M_150/N_150,M_180/N_180)*100,digits=1);
    T_Pcts <- round(c(T_30/N_30,T_60/N_60,T_90/N_90,T_120/N_120,T_150/N_150,T_180/N_180)*100,digits=1);
    return(list(C = C_Pcts, M = M_Pcts, T = T_Pcts));
}

##############################################################################
# Find the angle between the movement and the direction from the strike zone #
##############################################################################

Angle_Groupings <- function(PITCHES,QUAD,stand){
    # Find the indices of the pitches that are balls and are strikes
    Index <- ball_strike_prob(PITCHES,stand);
    # Set the indices of pitches outside the strike zone
    B_Ind <- Index$B;
    # Find the number of pitches outside the strike zone
    N_B <- length(B_Ind);
    # Set arrays for each set of angles for each type of movement
    Angle_PFX_B <- array(0,dim=N_B);
    Angle_PFX_Drag_B <- array(0,dim=N_B);
    Angle_W_B <- array(0,dim=N_B);
    Angle_W_Drag_B <- array(0,dim=N_B);
    # Loop over the pitches outside the zone and find the angle between the movement and direction of the pitch
    k <- 1;
    for (i in B_Ind){
        # Set the individual pitch data
        PITCH <- list(ax = PITCHES$ax[i], ay = PITCHES$ay[i], az = PITCHES$az[i], vx0 = PITCHES$vx0[i], vy0 = PITCHES$vy0[i], vz0 = PITCHES$vz0[i], x0 = PITCHES$x0[i], y0 = PITCHES$y0[i], z0 = PITCHES$z0[i]);
        # Set the location of the pitch at the front of home plate
        xz_loc <- xz_plate_loc(PITCH);
        # Find the four different types of movement
        PFX_B <- PFX_Mvt(PITCH,55);
        PFX_Drag_B <- PFX_Mvt_Drag(PITCH,55);
        W_B <- W_Mvt(PITCH,55);
        W_Drag_B <- W_Mvt_Drag(PITCH,55);
        # Find the direction of the pitch, relative to the strike zone
        xz_dir <- Dir_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the angle between the direction of the pitch and the movement
        Angle_PFX_B[k] <- Dot_Prod_Angle(xz_dir,PFX_B);
        Angle_PFX_Drag_B[k] <- Dot_Prod_Angle(xz_dir,PFX_Drag_B);
        Angle_W_B[k] <- Dot_Prod_Angle(xz_dir,W_B);
        Angle_W_Drag_B[k] <- Dot_Prod_Angle(xz_dir,W_Drag_B);
        k <- k+1;
    }
    return(list(PFX = Angle_PFX_B, PFX_D = Angle_PFX_Drag_B, W = Angle_W_B, W_D = Angle_W_Drag_B));
}

###################################################################################################
# Find the angle between the movement and the direction from the strike zone within some distance #
###################################################################################################

Angle_Groupings_Dist <- function(PITCHES,QUAD,stand){
    # Find the indices of the pitches that are balls and are strikes
    Index <- ball_strike_prob(PITCHES,stand);
    # Set the indices of pitches outside the strike zone
    B_Ind <- Index$B;
    # Find the number of pitches outside the strike zone
    N_B <- length(B_Ind);
    # Set arrays for each set of angles for each type of movement
    Angle_PFX_B <- array(0,dim=N_B);
    Angle_PFX_Drag_B <- array(0,dim=N_B);
    Angle_W_B <- array(0,dim=N_B);
    Angle_W_Drag_B <- array(0,dim=N_B);
    # Set an array for the distances from the strike zone
    SZ_Dist_B <- array(0,dim=N_B);
    # Loop over the pitches outside the zone and find the angle between the movement and direction of the pitch
    k <- 1;
    for (i in B_Ind){
        # Set the individual pitch data
        PITCH <- list(ax = PITCHES$ax[i], ay = PITCHES$ay[i], az = PITCHES$az[i], vx0 = PITCHES$vx0[i], vy0 = PITCHES$vy0[i], vz0 = PITCHES$vz0[i], x0 = PITCHES$x0[i], y0 = PITCHES$y0[i], z0 = PITCHES$z0[i]);
        # Set the location of the pitch at the front of home plate
        xz_loc <- xz_plate_loc(PITCH);
        # Find the four different types of movement
        PFX_B <- PFX_Mvt(PITCH,55);
        PFX_Drag_B <- PFX_Mvt_Drag(PITCH,55);
        W_B <- W_Mvt(PITCH,55);
        W_Drag_B <- W_Mvt_Drag(PITCH,55);
        # Find the direction of the pitch, relative to the strike zone
        xz_dir <- Dir_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the distance between the pitch and the core quadrilateral
        xz_sz_dist <- Min_Dist_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the angle between the direction of the pitch and the movement
        Angle_PFX_B[k] <- Dot_Prod_Angle(xz_dir,PFX_B);
        Angle_PFX_Drag_B[k] <- Dot_Prod_Angle(xz_dir,PFX_Drag_B);
        Angle_W_B[k] <- Dot_Prod_Angle(xz_dir,W_B);
        Angle_W_Drag_B[k] <- Dot_Prod_Angle(xz_dir,W_Drag_B);
        # Set an array for the distances from the strike zone
        SZ_Dist_B[k] <- xz_sz_dist - (log(2)/4)^(1/4);
        k <- k+1;
    }
    # Find the indices of the array where the pitches are within d feet of the strike zone
    d <- 1;
    Index_d <- which(SZ_Dist_B <= 1);
    return(list(PFX = Angle_PFX_B[Index_d], PFX_D = Angle_PFX_Drag_B[Index_d], W = Angle_W_B[Index_d], W_D = Angle_W_Drag_B[Index_d]));
}

##############################################
# Sort Percentages by Contact, Miss, or Take #
##############################################

Ct_Ms_Tk_Table <- function(CONTACT,MISS,TAKE,QUAD,stand){
    Angle_C <- Angle_Groupings(CONTACT,QUAD,stand);
    Angle_M <- Angle_Groupings(MISS,QUAD,stand);
    Angle_T <- Angle_Groupings(TAKE,QUAD,stand);
    PFX_Pct <- Angle_Pct_By_Case(Angle_C$PFX,Angle_M$PFX,Angle_T$PFX);
    PFX_D_Pct <- Angle_Pct_By_Case(Angle_C$PFX_D,Angle_M$PFX_D,Angle_T$PFX_D);
    W_Pct <- Angle_Pct_By_Case(Angle_C$W,Angle_M$W,Angle_T$W);
    W_D_Pct <- Angle_Pct_By_Case(Angle_C$W_D,Angle_M$W_D,Angle_T$W_D);
    df_ct_ms_tk <- data.frame(PFX_Pct$C,PFX_Pct$M,PFX_Pct$T,PFX_D_Pct$C,PFX_D_Pct$M,PFX_D_Pct$T,
                              W_Pct$C,W_Pct$M,W_Pct$T,W_D_Pct$C,W_D_Pct$M,W_D_Pct$T);
    print(df_ct_ms_tk);
}

###################################################################
# Sort Percentages by Contact, Miss, or Take For a Given Distance #
###################################################################

Ct_Ms_Tk_Dist_Table <- function(CONTACT,MISS,TAKE,QUAD,stand){
    Angle_C <- Angle_Groupings_Dist(CONTACT,QUAD,stand);
    Angle_M <- Angle_Groupings_Dist(MISS,QUAD,stand);
    Angle_T <- Angle_Groupings_Dist(TAKE,QUAD,stand);
    PFX_Pct <- Angle_Pct_By_Case(Angle_C$PFX,Angle_M$PFX,Angle_T$PFX);
    PFX_D_Pct <- Angle_Pct_By_Case(Angle_C$PFX_D,Angle_M$PFX_D,Angle_T$PFX_D);
    W_Pct <- Angle_Pct_By_Case(Angle_C$W,Angle_M$W,Angle_T$W);
    W_D_Pct <- Angle_Pct_By_Case(Angle_C$W_D,Angle_M$W_D,Angle_T$W_D);
    df_ct_ms_tk <- data.frame(PFX_Pct$C,PFX_Pct$M,PFX_Pct$T,PFX_D_Pct$C,PFX_D_Pct$M,PFX_D_Pct$T,
                              W_Pct$C,W_Pct$M,W_Pct$T,W_D_Pct$C,W_D_Pct$M,W_D_Pct$T);
    print(df_ct_ms_tk);
}

###################################################################################
# Track the projected strike probability with remaining PITCHf/x movement removed #
###################################################################################

Pitch_Project_PFX <- function(PITCHES,Ny,stand){
    g <- -32.174;
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Set a list for the PITCHf/x data
        PITCH <- list(ax = ax[i], ay = ay[i], az = az[i], vx0 = vx0[i], vy0 = vy0[i], vz0 = vz0[i], x0 = x0[i], y0 = y0[i], z0 = z0[i]);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the projected location at the front of the plate
            xz_PL <- PFX_Calc(PITCH,t_inc,t_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];
        }
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc));
}

#####################################################################################################
# Track the projected strike probability with remaining PITCHf/x movement removed but drag added in #
#####################################################################################################

Pitch_Project_PFX_Drag <- function(PITCHES,Ny,stand){
    g <- -32.174;
    G <- c(0,0,g);
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the velocity at the front of the plate
        vx_plate <- ax[i]*t_plate + vx0[i];
        vy_plate <- ay[i]*t_plate + vy0[i];
        vz_plate <- az[i]*t_plate + vz0[i];
        v_plate <- list(x = vx_plate, y = vy_plate, z = vz_plate);
        # Set a list for the PITCHf/x data
        PITCH <- list(ax = ax[i], ay = ay[i], az = az[i], vx0 = vx0[i], vy0 = vy0[i], vz0 = vz0[i], x0 = x0[i], y0 = y0[i], z0 = z0[i]);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the projected location at the front of the plate
            xz_PL <- PFX_Drag_Calc(PITCH,t_inc,t_plate,v_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];

        }
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc));
}

#######################################################################################
# Track the planar projected strike probability with the remaining w-movement removed #
#######################################################################################

Pitch_Project_W <- function(PITCHES,Ny,stand){
    g <- -32.174;
    G <- c(0,0,g);
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the x, y, and z coordinates of the pitch at the front of the plate
        x_plate <- x0[i] + vx0[i]*t_plate + 0.5*ax[i]*t_plate^2;
        y_plate <- y0[i] + vy0[i]*t_plate + 0.5*ay[i]*t_plate^2;
        z_plate <- z0[i] + vz0[i]*t_plate + 0.5*az[i]*t_plate^2;
        # Set up a list containing the parameters of the pitch
        PITCH <- list(ax=ax[i], ay=ay[i], az=az[i], vx0=vx0[i], vy0=vy0[i], vz0=vz0[i], x0=x0[i], y0=y0[i], z0=z0[i]);
        BV <- PFX_Basis(PITCH);
        BV_T = t(BV);
        # Find the u and b coordinates of the pitch at the front of the plate
        b_plate <- BV_T[1,] %*% c(x_plate,y_plate,z_plate);
        u_plate <- BV_T[2,] %*% c(x_plate,y_plate,z_plate);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the projected location at the front of the plate
            xz_PL <- W_Calc(PITCH,BV_T,t_inc,t_plate,b_plate,u_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];
        }
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc));
}

#########################################################################################################
# Track the planar projected strike probability with the remaining w-movement removed but drag added in #
#########################################################################################################

Pitch_Project_W_Drag <- function(PITCHES,Ny,stand){
    g <- -32.174;
    G <- c(0,0,g);
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the x, y, and z coordinates of the pitch at the front of the plate
        x_plate <- x0[i] + vx0[i]*t_plate + 0.5*ax[i]*t_plate^2;
        y_plate <- y0[i] + vy0[i]*t_plate + 0.5*ay[i]*t_plate^2;
        z_plate <- z0[i] + vz0[i]*t_plate + 0.5*az[i]*t_plate^2;
        # Find the velocity at the front of the plate
        vx_plate <- ax[i]*t_plate + vx0[i];
        vy_plate <- ay[i]*t_plate + vy0[i];
        vz_plate <- az[i]*t_plate + vz0[i];
        v_plate <- list(x = vx_plate, y = vy_plate, z = vz_plate);
        # Set up a list containing the parameters of the pitch
        PITCH <- list(ax=ax[i], ay=ay[i], az=az[i], vx0=vx0[i], vy0=vy0[i], vz0=vz0[i], x0=x0[i], y0=y0[i], z0=z0[i]);
        BV <- PFX_Basis(PITCH);
        BV_T = t(BV);
        # Find the u and b coordinates of the pitch at the front of the plate
        b_plate <- BV_T[1,] %*% c(x_plate,y_plate,z_plate);
        u_plate <- BV_T[2,] %*% c(x_plate,y_plate,z_plate);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the projected location at the front of the plate
            xz_PL <- W_Drag_Calc(PITCH,BV_T,t_inc,t_plate,v_plate,b_plate,u_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];
        }
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc));
}

################################################################################
# Track the projected pitches and movement direction with PFX movement removed #
################################################################################

Pitch_Project_Dir_PFX <- function(PITCHES,QUAD,Ny,stand){
    # Set the arrays
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    Zone <- array(0,dim=Np);
    x_Dir <- array(0,dim=c(Np,Ny));
    z_Dir <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Compute the time to the plate
        c_plate <- y0[i]-(17/12);
        b <- vy0[i];
        a <- 0.5*ay[i];
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the location of the pitch at the plate
        PITCH <- list(ax = ax[i], ay = ay[i], az = az[i], vx0 = vx0[i], vy0 = vy0[i], vz0 = vz0[i], x0 = x0[i], y0 = y0[i], z0 = z0[i]);
        xz_loc <- xz_plate_loc(PITCH);
        # Find the zone (1-9) of the pitch at the plate
        Zone[i] <- Zone_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Loop over the y-slices from release to the plate
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the projected location at the front of the plate
            xz_PL <- PFX_Calc(PITCH,t_inc,t_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];
            # Find where the pitch projects from one inch ahead of y_inc
            c_F <- y0[i]-y_seq[j] + (1/12);
            t_F <- (-b - sqrt(b^2 - 4*a*c_F))/(2*a);
            xz_F <- PFX_Calc(PITCH,t_F,t_plate);
            # Find where the pitch projects from one inch behind of y_inc
            c_B <- y0[i]-y_seq[j] - (1/12);
            t_B <- (-b - sqrt(b^2 - 4*a*c_B))/(2*a);
            xz_B <- PFX_Calc(PITCH,t_B,t_plate);
            x_Dir[i,j] <- (xz_F[1] - xz_B[1])/(1/6);
            z_Dir[i,j] <- (xz_F[2] - xz_B[2])/(1/6);
        }
        x_Dir[i,Ny] <- (x_Proj_Loc[i,Ny] - xz_B[1])/(1/12);
        z_Dir[i,Ny] <- (z_Proj_Loc[i,Ny] - xz_B[2])/(1/12);
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc, x_Dir = x_Dir, z_Dir = z_Dir, Zone = Zone));
}

#####################################################################################################################
# Track the projected pitches and movement direction with the remaining PITCHf/x movement removed but drag added in #
#####################################################################################################################

Pitch_Project_Dir_PFX_Drag <- function(PITCHES,QUAD,Ny,stand){
    # Set the arrays
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    Zone <- array(0,dim=Np);
    x_Dir <- array(0,dim=c(Np,Ny));
    z_Dir <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Compute the time to the plate
        c_plate <- y0[i]-(17/12);
        b <- vy0[i];
        a <- 0.5*ay[i];
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the velocity at the front of the plate
        vx_plate <- ax[i]*t_plate + vx0[i];
        vy_plate <- ay[i]*t_plate + vy0[i];
        vz_plate <- az[i]*t_plate + vz0[i];
        v_plate <- list(x = vx_plate, y = vy_plate, z = vz_plate);
        # Find the location of the pitch at the plate
        PITCH <- list(ax = ax[i], ay = ay[i], az = az[i], vx0 = vx0[i], vy0 = vy0[i], vz0 = vz0[i], x0 = x0[i], y0 = y0[i], z0 = z0[i]);
        xz_loc <- xz_plate_loc(PITCH);
        # Find the zone (1-9) of the pitch at the plate
        Zone[i] <- Zone_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Loop over the y-slices from release to the plate
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the projected location at the front of the plate
            xz_PL <- PFX_Drag_Calc(PITCH,t_inc,t_plate,v_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];
            # Find where the pitch projects from one inch ahead of y_inc
            c_F <- y0[i]-y_seq[j] + (1/12);
            t_F <- (-b - sqrt(b^2 - 4*a*c_F))/(2*a);
            xz_F <- PFX_Drag_Calc(PITCH,t_F,t_plate,v_plate);
            # Find where the pitch projects from one inch behind of y_inc
            c_B <- y0[i]-y_seq[j] - (1/12);
            t_B <- (-b - sqrt(b^2 - 4*a*c_B))/(2*a);
            xz_B <- PFX_Drag_Calc(PITCH,t_B,t_plate,v_plate);
            x_Dir[i,j] <- (xz_F[1] - xz_B[1])/(1/6);
            z_Dir[i,j] <- (xz_F[2] - xz_B[2])/(1/6);
        }
        x_Dir[i,Ny] <- (x_Proj_Loc[i,Ny] - xz_B[1])/(1/12);
        z_Dir[i,Ny] <- (z_Proj_Loc[i,Ny] - xz_B[2])/(1/12);
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc, x_Dir = x_Dir, z_Dir = z_Dir, Zone = Zone));
}


#################################################################################################
# Track the projected pitches and movement direction with the remaining planar movement removed #
#################################################################################################

Pitch_Project_Dir_W <- function(PITCHES,QUAD,Ny,stand){
    g <- -32.174;
    G <- c(0,0,g);
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    Zone <- array(0,dim=Np);
    x_Dir <- array(0,dim=c(Np,Ny));
    z_Dir <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the x, y, and z coordinates of the pitch at the front of the plate
        x_plate <- x0[i] + vx0[i]*t_plate + 0.5*ax[i]*t_plate^2;
        y_plate <- y0[i] + vy0[i]*t_plate + 0.5*ay[i]*t_plate^2;
        z_plate <- z0[i] + vz0[i]*t_plate + 0.5*az[i]*t_plate^2;
        # Set up a list containing the parameters of the pitch
        PITCH <- list(ax=ax[i], ay=ay[i], az=az[i], vx0=vx0[i], vy0=vy0[i], vz0=vz0[i], x0=x0[i], y0=y0[i], z0=z0[i]);
        # Find the location of the pitch at the plate
        xz_loc <- c(x_plate,z_plate);
        # Find the zone (1-9) of the pitch at the plate
        Zone[i] <- Zone_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the coordinate transformation
        BV <- PFX_Basis(PITCH);
        BV_T = t(BV);
        # Find the u and b coordinates of the pitch at the front of the plate
        b_plate <- BV_T[1,] %*% c(x_plate,y_plate,z_plate);
        u_plate <- BV_T[2,] %*% c(x_plate,y_plate,z_plate);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the location at the front of the plate
            xz_PL <- W_Calc(PITCH,BV_T,t_inc,t_plate,b_plate,u_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];
            # Find where the pitch projects from one inch ahead of y_inc
            c_F <- y0[i]-y_seq[j] + (1/12);
            t_F <- (-b - sqrt(b^2 - 4*a*c_F))/(2*a);
            xz_F <- W_Calc(PITCH,BV_T,t_F,t_plate,b_plate,u_plate);
            # Find where the pitch projects from one inch behind of y_inc
            c_B <- y0[i]-y_seq[j] - (1/12);
            t_B <- (-b - sqrt(b^2 - 4*a*c_B))/(2*a);
            xz_B <- W_Calc(PITCH,BV_T,t_B,t_plate,b_plate,u_plate);
            x_Dir[i,j] <- (xz_F[1] - xz_B[1])/(1/6);
            z_Dir[i,j] <- (xz_F[2] - xz_B[2])/(1/6);
        }
        x_Dir[i,Ny] <- (x_Proj_Loc[i,Ny] - xz_B[1])/(1/12);
        z_Dir[i,Ny] <- (z_Proj_Loc[i,Ny] - xz_B[2])/(1/12);
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc, x_Dir = x_Dir, z_Dir = z_Dir, Zone = Zone));
}

#######################################################################################################################
# Track the planar projected location and movement direction with remaining planar movement removed but drag added in #
#######################################################################################################################

Pitch_Project_Dir_W_Drag <- function(PITCHES,QUAD,Ny,stand){
    g <- -32.174;
    G <- c(0,0,g);
    Np <- length(PITCHES$x0);
    x0 <- PITCHES$x0;
    y0 <- PITCHES$y0;
    z0 <- PITCHES$z0;
    vx0 <- PITCHES$vx0;
    vy0 <- PITCHES$vy0;
    vz0 <- PITCHES$vz0;
    ax <- PITCHES$ax;
    ay <- PITCHES$ay;
    az <- PITCHES$az;
    y_seq <- seq(55,17/12,length=Ny);
    x_Proj_Loc <- array(0,dim=c(Np,Ny));
    z_Proj_Loc <- array(0,dim=c(Np,Ny));
    Zone <- array(0,dim=Np);
    x_Dir <- array(0,dim=c(Np,Ny));
    z_Dir <- array(0,dim=c(Np,Ny));
    for (i in 1:Np){
        # Find the time to the front of home plate
        a <- 0.5*ay[i];
        b <- vy0[i];
        c_plate <- y0[i]-(17/12);
        t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
        # Find the x, y, and z coordinates of the pitch at the front of the plate
        x_plate <- x0[i] + vx0[i]*t_plate + 0.5*ax[i]*t_plate^2;
        y_plate <- y0[i] + vy0[i]*t_plate + 0.5*ay[i]*t_plate^2;
        z_plate <- z0[i] + vz0[i]*t_plate + 0.5*az[i]*t_plate^2;
        # Find the velocity at the front of the plate
        vx_plate <- ax[i]*t_plate + vx0[i];
        vy_plate <- ay[i]*t_plate + vy0[i];
        vz_plate <- az[i]*t_plate + vz0[i];
        v_plate <- list(x = vx_plate, y = vy_plate, z = vz_plate);
        # Set up a list containing the parameters of the pitch
        PITCH <- list(ax=ax[i], ay=ay[i], az=az[i], vx0=vx0[i], vy0=vy0[i], vz0=vz0[i], x0=x0[i], y0=y0[i], z0=z0[i]);
        # Find the location of the pitch at the plate
        xz_loc <- c(x_plate,z_plate);
        # Find the zone (1-9) of the pitch at the plate
        Zone[i] <- Zone_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the coordinate transformation
        BV <- PFX_Basis(PITCH);
        BV_T = t(BV);
        # Find the u and b coordinates of the pitch at the front of the plate
        b_plate <- BV_T[1,] %*% c(x_plate,y_plate,z_plate);
        u_plate <- BV_T[2,] %*% c(x_plate,y_plate,z_plate);
        # Loop over the y-slices
        for (j in 1:Ny){
            c_inc <- y0[i]-y_seq[j];
            t_inc <- (-b - sqrt(b^2 - 4*a*c_inc))/(2*a);
            # Find the projected location at the front of the plate
            xz_PL <- W_Drag_Calc(PITCH,BV_T,t_inc,t_plate,v_plate,b_plate,u_plate);
            x_Proj_Loc[i,j] <- xz_PL[1];
            z_Proj_Loc[i,j] <- xz_PL[2];
            # Find where the pitch projects from one inch ahead of y_inc
            c_F <- y0[i]-y_seq[j] + (1/12);
            t_F <- (-b - sqrt(b^2 - 4*a*c_F))/(2*a);
            xz_F <- W_Drag_Calc(PITCH,BV_T,t_F,t_plate,v_plate,b_plate,u_plate);
            # Find where the pitch projects from one inch behind of y_inc
            c_B <- y0[i]-y_seq[j] - (1/12);
            t_B <- (-b - sqrt(b^2 - 4*a*c_B))/(2*a);
            xz_B <- W_Drag_Calc(PITCH,BV_T,t_B,t_plate,v_plate,b_plate,u_plate);
            x_Dir[i,j] <- (xz_F[1] - xz_B[1])/(1/6);
            z_Dir[i,j] <- (xz_F[2] - xz_B[2])/(1/6);
        }
        x_Dir[i,Ny] <- (x_Proj_Loc[i,Ny] - xz_B[1])/(1/12);
        z_Dir[i,Ny] <- (z_Proj_Loc[i,Ny] - xz_B[2])/(1/12);
    }
    return(list(x_PL = x_Proj_Loc, z_PL = z_Proj_Loc, x_Dir = x_Dir, z_Dir = z_Dir, Zone = Zone));
}

#################################################
# Find the projection with PFX movement removed #
#################################################

PFX_Calc <- function(PITCH,t_loc,t_plate){
    # Set gravity
    g <- -32.174;
    # Compute x and z location in the quadratic part of the curve
    x_loc <- PITCH$x0 + PITCH$vx0*t_loc + 0.5*PITCH$ax*t_loc^2;
    z_loc <- PITCH$z0 + PITCH$vz0*t_loc + 0.5*PITCH$az*t_loc^2;
    # Find the velocity of the pitch in x and z at t_loc
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the projected location of the pitch at the front of home plate
    x_PL <- x_loc + vx_loc*(t_plate - t_loc);
    z_PL <- z_loc + vz_loc*(t_plate - t_loc) + 0.5*g*(t_plate - t_loc)^2;
    return(c(x_PL,z_PL));
}

###################################################################
# Find the projection with PFX movement removed and drag added in #
###################################################################

PFX_Drag_Calc <- function(PITCH,t_loc,t_plate,v_plate){
    # Set gravity
    g <- -32.174;
    G <- c(0,0,g);
    # Compute x and z location in the quadratic part of the curve
    x_loc <- PITCH$x0 + PITCH$vx0*t_loc + 0.5*PITCH$ax*t_loc^2;
    z_loc <- PITCH$z0 + PITCH$vz0*t_loc + 0.5*PITCH$az*t_loc^2;
    # Find the velocity of the pitch in x and z at t_inc
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vy_loc <- PITCH$vy0 + PITCH$ay*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the average velocity over the remaining distance to the plate
    vx_avg <- (vx_loc + v_plate$x)/2;
    vy_avg <- (vy_loc + v_plate$y)/2;
    vz_avg <- (vz_loc + v_plate$z)/2;
    v_avg <- sqrt(vx_avg^2 + vy_avg^2 + vz_avg^2);
    # Set the average velocity and acceleration vectors
    V <- c(vx_avg,vy_avg,vz_avg);
    A <- c(PITCH$ax,PITCH$ay,PITCH$az);
    # Normalize the velocity vector
    V_norm <- V/sqrt(V%*%V);
    # Find the drag vector
    Drag <- -abs((A-G)%*%V_norm)*V_norm;
    # Find the projected location at the front of the plate
    x_PL <- x_loc + vx_loc*(t_plate-t_loc) + 0.5*Drag[1]*(t_plate-t_loc)^2;
    z_PL <- z_loc + vz_loc*(t_plate-t_loc) + 0.5*(g + Drag[3])*(t_plate-t_loc)^2;
    return(c(x_PL,z_PL));
}

####################################################
# Find the projection with planar movement removed #
####################################################

W_Calc <- function(PITCH,BV_T,t_loc,t_plate,b_plate,u_plate){
    # Set gravity
    G <- c(0,0,-32.174);
    x_loc <- PITCH$x0 + PITCH$vx0*t_loc + 0.5*PITCH$ax*t_loc^2;
    y_loc <- PITCH$y0 + PITCH$vy0*t_loc + 0.5*PITCH$ay*t_loc^2;
    z_loc <- PITCH$z0 + PITCH$vz0*t_loc + 0.5*PITCH$az*t_loc^2;
    # Find the coordinates of the point in (u,w,b)-space
    w_loc <- BV_T[3,] %*% c(x_loc,y_loc,z_loc);
    # Find the velocity at y feet
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vy_loc <- PITCH$vy0 + PITCH$ay*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the velocity in (u,w,b)-space
    vw_loc <- BV_T[3,] %*% c(vx_loc,vy_loc,vz_loc);
    # Find the projection of gravity in the w direction
    gw <- BV_T[3,] %*% G;
    # Find the projected location of the pitch at the front of home plate
    pw <- w_loc + vw_loc*(t_plate-t_loc) + 0.5*gw*(t_plate - t_loc)^2;
    x_PL <- BV_T[,1] %*% c(b_plate,u_plate,pw);
    z_PL <- BV_T[,3] %*% c(b_plate,u_plate,pw);
    return(c(x_PL,z_PL));
}

####################################################################
# Find the projection of planar movement removed and drag added in #
####################################################################

W_Drag_Calc <- function(PITCH,BV_T,t_loc,t_plate,v_plate,b_plate,u_plate){
    # Set gravity
    G <- c(0,0,-32.174);
    x_loc <- PITCH$x0 + PITCH$vx0*t_loc + 0.5*PITCH$ax*t_loc^2;
    y_loc <- PITCH$y0 + PITCH$vy0*t_loc + 0.5*PITCH$ay*t_loc^2;
    z_loc <- PITCH$z0 + PITCH$vz0*t_loc + 0.5*PITCH$az*t_loc^2;
    # Find the coordinates of the point in (u,w,b)-space
    w_loc <- BV_T[3,] %*% c(x_loc,y_loc,z_loc);
    # Find the velocity at y feet
    vx_loc <- PITCH$vx0 + PITCH$ax*t_loc;
    vy_loc <- PITCH$vy0 + PITCH$ay*t_loc;
    vz_loc <- PITCH$vz0 + PITCH$az*t_loc;
    # Find the velocity in (u,w,b)-space
    vw_loc <- BV_T[3,] %*% c(vx_loc,vy_loc,vz_loc);
    # Find the average velocity over the remaining distance to the plate
    vx_avg <- (vx_loc + v_plate$x)/2;
    vy_avg <- (vy_loc + v_plate$y)/2;
    vz_avg <- (vz_loc + v_plate$z)/2;
    v_avg <- sqrt(vx_avg^2 + vy_avg^2 + vz_avg^2);
    # Set the average velocity and acceleration vectors
    V <- c(vx_avg,vy_avg,vz_avg);
    A <- c(PITCH$ax,PITCH$ay,PITCH$az);
    # Normalize the velocity vector
    V_norm <- V/sqrt(V%*%V);
    # Find the drag vector
    Drag <- -abs((A-G)%*%V_norm)*V_norm;
    # Find the drag in the w-direction
    Drag_w <- BV_T[3,] %*% Drag;
    # Find the projection of gravity in the w direction
    gw <- BV_T[3,] %*% G;
    # Find the point on the tangent line passing through the plane of home plate
    pw <- w_loc + vw_loc*(t_plate-t_loc) + 0.5*(gw + Drag_w)*(t_plate - t_loc)^2;
    x_PL <- BV_T[,1] %*% c(b_plate,u_plate,pw);
    z_PL <- BV_T[,3] %*% c(b_plate,u_plate,pw);
    return(c(x_PL,z_PL));
}

#########################################################################################
# Plot the pitches outside the strike zone based on angle vs. movement for all pitchers #
#########################################################################################

Angle_Dist_Table_Plot_All <- function(PITCHES,QUAD,throws,stand,pitch_type,result,res_color){
    # Find the indices of the pitches that are balls and are strikes
    Index <- ball_strike_prob(PITCHES,stand);
    # Set the indices of pitches outside the strike zone
    B_Ind <- Index$B;
    # Find the number of pitches outside the strike zone
    N_B <- length(B_Ind);
    # Set arrays for each set of angles for each type of movement
    Angle_PFX_B <- array(0,dim=N_B);
    Angle_PFX_Drag_B <- array(0,dim=N_B);
    Angle_W_B <- array(0,dim=N_B);
    Angle_W_Drag_B <- array(0,dim=N_B);
    # Set an array for the distances from the strike zone
    SZ_Dist_B <- array(0,dim=N_B);
    # Loop over the pitches outside the zone and find the angle between the movement and direction of the pitch
    k <- 1;
    for (i in B_Ind){
        # Set the individual pitch data
        PITCH <- list(ax = PITCHES$ax[i], ay = PITCHES$ay[i], az = PITCHES$az[i], vx0 = PITCHES$vx0[i], vy0 = PITCHES$vy0[i], vz0 = PITCHES$vz0[i], x0 = PITCHES$x0[i], y0 = PITCHES$y0[i], z0 = PITCHES$z0[i]);
        # Set the location of the pitch at the front of home plate
        xz_loc <- xz_plate_loc(PITCH);
        # Find the four different types of movement
        PFX_B <- PFX_Mvt(PITCH,55);
        PFX_Drag_B <- PFX_Mvt_Drag(PITCH,55);
        W_B <- W_Mvt(PITCH,55);
        W_Drag_B <- W_Mvt_Drag(PITCH,55);
        # Find the direction of the pitch, relative to the strike zone
        xz_dir <- Dir_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the distance between the pitch and the core quadrilateral
        xz_sz_dist <- Min_Dist_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the angle between the direction of the pitch and the movement
        Angle_PFX_B[k] <- Dot_Prod_Angle(xz_dir,PFX_B);
        Angle_PFX_Drag_B[k] <- Dot_Prod_Angle(xz_dir,PFX_Drag_B);
        Angle_W_B[k] <- Dot_Prod_Angle(xz_dir,W_B);
        Angle_W_Drag_B[k] <- Dot_Prod_Angle(xz_dir,W_Drag_B);
        # Set an array for the distances from the strike zone
        SZ_Dist_B[k] <- xz_sz_dist - (log(2)/4)^(1/4);
        k <- k+1;
    }
    Angle_Pct_Dist(Angle_PFX_B,SZ_Dist_B);
    Angle_Pct_Dist(Angle_PFX_Drag_B,SZ_Dist_B);
    Angle_Pct_Dist(Angle_W_B,SZ_Dist_B);
    Angle_Pct_Dist(Angle_W_Drag_B,SZ_Dist_B);
    # Plot the data for PITCHf/x movement
    windows();
    PFX_Title <- paste("2016: ",pitch_name," ",result," Outside (",throws,"HP v. ",stand,"HB) PFX",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_PFX_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(PFX_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
    # Plot the data for PITCHf/x movement with drag removed
    windows();
    PFX_Drag_Title <- paste("2016: ",pitch_name," ",result," Outside (",throws,"HP v. ",stand,"HB) PFX D",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_PFX_Drag_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(PFX_Drag_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
    # Plot the data for planar movement
    windows();
    W_Title <- paste("2016: ",pitch_name," ",result," Outside (",throws,"HP v. ",stand,"HB) W",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_W_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(W_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
    # Plot the data for planar movement with drag removed
    windows();
    W_Drag_Title <- paste("2016: ",pitch_name," ",result," Outside (",throws,"HP v. ",stand,"HB) W D",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_W_Drag_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(W_Drag_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
}

######################################################################################
# Plot the pitches outside the strike zone based on angle vs. movement for a pitcher #
######################################################################################

Angle_Dist_Table_Plot <- function(PITCHES,QUAD,first,last,stand,pitch_type,result,res_color){
    # Find the indices of the pitches that are balls and are strikes
    Index <- ball_strike_prob(PITCHES,stand);
    # Set the indices of pitches outside the strike zone
    B_Ind <- Index$B;
    # Find the number of pitches outside the strike zone
    N_B <- length(B_Ind);
    # Set arrays for each set of angles for each type of movement
    Angle_PFX_B <- array(0,dim=N_B);
    Angle_PFX_Drag_B <- array(0,dim=N_B);
    Angle_W_B <- array(0,dim=N_B);
    Angle_W_Drag_B <- array(0,dim=N_B);
    # Set an array for the distances from the strike zone
    SZ_Dist_B <- array(0,dim=N_B);
    # Loop over the pitches outside the zone and find the angle between the movement and direction of the pitch
    k <- 1;
    for (i in B_Ind){
        # Set the individual pitch data
        PITCH <- list(ax = PITCHES$ax[i], ay = PITCHES$ay[i], az = PITCHES$az[i], vx0 = PITCHES$vx0[i], vy0 = PITCHES$vy0[i], vz0 = PITCHES$vz0[i], x0 = PITCHES$x0[i], y0 = PITCHES$y0[i], z0 = PITCHES$z0[i]);
        # Set the location of the pitch at the front of home plate
        xz_loc <- xz_plate_loc(PITCH);
        # Find the four different types of movement
        PFX_B <- PFX_Mvt(PITCH,55);
        PFX_Drag_B <- PFX_Mvt_Drag(PITCH,55);
        W_B <- W_Mvt(PITCH,55);
        W_Drag_B <- W_Mvt_Drag(PITCH,55);
        # Find the direction of the pitch, relative to the strike zone
        xz_dir <- Dir_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the distance between the pitch and the core quadrilateral
        xz_sz_dist <- Min_Dist_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the angle between the direction of the pitch and the movement
        Angle_PFX_B[k] <- Dot_Prod_Angle(xz_dir,PFX_B);
        Angle_PFX_Drag_B[k] <- Dot_Prod_Angle(xz_dir,PFX_Drag_B);
        Angle_W_B[k] <- Dot_Prod_Angle(xz_dir,W_B);
        Angle_W_Drag_B[k] <- Dot_Prod_Angle(xz_dir,W_Drag_B);
        # Set an array for the distances from the strike zone
        SZ_Dist_B[k] <- xz_sz_dist - (log(2)/4)^(1/4);
        k <- k+1;
    }
    Angle_Pct_Dist(Angle_PFX_B,SZ_Dist_B);
    Angle_Pct_Dist(Angle_PFX_Drag_B,SZ_Dist_B);
    Angle_Pct_Dist(Angle_W_B,SZ_Dist_B);
    Angle_Pct_Dist(Angle_W_Drag_B,SZ_Dist_B);
    # Plot the data for PITCHf/x movement
    windows();
    PFX_Title <- paste("2016: ",pitch_name," ",result," Outside (",first," ",last," v. ",stand,"HB) PFX",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_PFX_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(PFX_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
    # Plot the data for PITCHf/x movement with drag removed
    windows();
    PFX_Drag_Title <- paste("2016: ",pitch_name," ",result," Outside (",first," ",last," v. ",stand,"HB) PFX D",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_PFX_Drag_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(PFX_Drag_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
    # Plot the data for planar movement
    windows();
    W_Title <- paste("2016: ",pitch_name," ",result," Outside (",first," ",last," v. ",stand,"HB) W",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_W_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(W_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
    # Plot the data for planar movement with drag removed
    windows();
    W_Drag_Title <- paste("2016: ",pitch_name," ",result," Outside (",first," ",last," v. ",stand,"HB) W D",sep="");
    data.df <- data.frame(Dist = SZ_Dist_B,Angle = Angle_W_Drag_B);
    print(ggplot(data.df,aes(x=Dist,y=Angle)) + geom_hex(binwidth=c(0.25,15)) +
          scale_x_continuous(breaks=seq(0,3,0.5),limits=c(-0.125,3.125)) +
          scale_y_continuous(breaks=seq(0,180,30),limits=c(-7.5,187.5)) +
          xlab("Distance From Strike Zone (Feet)") + ylab("Angle (Degrees)") +
          scale_fill_gradient(low="white", high=res_color) + ggtitle(W_Drag_Title) +
      theme(panel.background = element_rect(fill = "black", color = "black", size = 1)));
}

#########################################################################################
# Calculate the angle of the pitches outside the strike zone vs. movement for a pitcher #
#########################################################################################

Angle_Dist_Store <- function(PITCHES,QUAD,stand,pitch_type){
    # Find the indices of the pitches that are balls and are strikes
    Index <- ball_strike_prob(PITCHES,stand);
    # Set the indices of pitches outside the strike zone
    B_Ind <- Index$B;
    # Find the number of pitches outside the strike zone
    N_B <- length(B_Ind);
    # Set arrays for each set of angles for each type of movement
    Angle_PFX_B <- array(0,dim=N_B);
    Angle_PFX_Drag_B <- array(0,dim=N_B);
    Angle_W_B <- array(0,dim=N_B);
    Angle_W_Drag_B <- array(0,dim=N_B);
    # Set an array for the distances from the strike zone
    SZ_Dist_B <- array(0,dim=N_B);
    # Loop over the pitches outside the zone and find the angle between the movement and direction of the pitch
    k <- 1;
    for (i in B_Ind){
        # Set the individual pitch data
        PITCH <- list(ax = PITCHES$ax[i], ay = PITCHES$ay[i], az = PITCHES$az[i], vx0 = PITCHES$vx0[i], vy0 = PITCHES$vy0[i], vz0 = PITCHES$vz0[i], x0 = PITCHES$x0[i], y0 = PITCHES$y0[i], z0 = PITCHES$z0[i]);
        # Set the location of the pitch at the front of home plate
        xz_loc <- xz_plate_loc(PITCH);
        # Find the four different types of movement
        PFX_B <- PFX_Mvt(PITCH,55);
        PFX_Drag_B <- PFX_Mvt_Drag(PITCH,55);
        W_B <- W_Mvt(PITCH,55);
        W_Drag_B <- W_Mvt_Drag(PITCH,55);
        # Find the direction of the pitch, relative to the strike zone
        xz_dir <- Dir_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the distance between the pitch and the core quadrilateral
        xz_sz_dist <- Min_Dist_Quad(xz_loc,QUAD$UL,QUAD$UR,QUAD$LL,QUAD$LR);
        # Find the angle between the direction of the pitch and the movement
        Angle_PFX_B[k] <- Dot_Prod_Angle(xz_dir,PFX_B);
        Angle_PFX_Drag_B[k] <- Dot_Prod_Angle(xz_dir,PFX_Drag_B);
        Angle_W_B[k] <- Dot_Prod_Angle(xz_dir,W_B);
        Angle_W_Drag_B[k] <- Dot_Prod_Angle(xz_dir,W_Drag_B);
        # Set an array for the distances from the strike zone
        SZ_Dist_B[k] <- xz_sz_dist - (log(2)/4)^(1/4);
        k <- k+1;
    }
    PFX_Store <- Angle_Pct_Dist_Store(Angle_PFX_B,SZ_Dist_B);
    PFX_Drag_Store <- Angle_Pct_Dist_Store(Angle_PFX_Drag_B,SZ_Dist_B);
    W_Store <- Angle_Pct_Dist_Store(Angle_W_B,SZ_Dist_B);
    W_Drag_Store <- Angle_Pct_Dist_Store(Angle_W_Drag_B,SZ_Dist_B);
    return(c(PFX_Store,PFX_Drag_Store,W_Store,W_Drag_Store,N_B));
}

#####################################################
# Store the angle and distance data in a data frame #
#####################################################

Angle_Dist_to_DF <- function(First,Last,SR_Out,MR_Out,Angle_Dist){
    # Separate out the angles
    PFX_30_A <- Angle_Dist[,1];
    PFX_60_A <- Angle_Dist[,2];
    PFX_90_A <- Angle_Dist[,3];
    PFX_120_A <- Angle_Dist[,4];
    PFX_150_A <- Angle_Dist[,5];
    PFX_180_A <- Angle_Dist[,6];
    PFX_Drag_30_A <- Angle_Dist[,13];
    PFX_Drag_60_A <- Angle_Dist[,14];
    PFX_Drag_90_A <- Angle_Dist[,15];
    PFX_Drag_120_A <- Angle_Dist[,16];
    PFX_Drag_150_A <- Angle_Dist[,17];
    PFX_Drag_180_A <- Angle_Dist[,18];
    W_30_A <- Angle_Dist[,25];
    W_60_A <- Angle_Dist[,26];
    W_90_A <- Angle_Dist[,27];
    W_120_A <- Angle_Dist[,28];
    W_150_A <- Angle_Dist[,29];
    W_180_A <- Angle_Dist[,30];
    W_Drag_30_A <- Angle_Dist[,37];
    W_Drag_60_A <- Angle_Dist[,38];
    W_Drag_90_A <- Angle_Dist[,39];
    W_Drag_120_A <- Angle_Dist[,40];
    W_Drag_150_A <- Angle_Dist[,41];
    W_Drag_180_A <- Angle_Dist[,42];
    # Separate out the distances
    PFX_30_D <- Angle_Dist[,7];
    PFX_60_D <- Angle_Dist[,8];
    PFX_90_D <- Angle_Dist[,9];
    PFX_120_D <- Angle_Dist[,10];
    PFX_150_D <- Angle_Dist[,11];
    PFX_180_D <- Angle_Dist[,12];
    PFX_Drag_30_D <- Angle_Dist[,19];
    PFX_Drag_60_D <- Angle_Dist[,20];
    PFX_Drag_90_D <- Angle_Dist[,21];
    PFX_Drag_120_D <- Angle_Dist[,22];
    PFX_Drag_150_D <- Angle_Dist[,23];
    PFX_Drag_180_D <- Angle_Dist[,24];
    W_30_D <- Angle_Dist[,31];
    W_60_D <- Angle_Dist[,32];
    W_90_D <- Angle_Dist[,33];
    W_120_D <- Angle_Dist[,34];
    W_150_D <- Angle_Dist[,35];
    W_180_D <- Angle_Dist[,36];
    W_Drag_30_D <- Angle_Dist[,43];
    W_Drag_60_D <- Angle_Dist[,44];
    W_Drag_90_D <- Angle_Dist[,45];
    W_Drag_120_D <- Angle_Dist[,46];
    W_Drag_150_D <- Angle_Dist[,47];
    W_Drag_180_D <- Angle_Dist[,48];
    # Store the angles and distances in a data frame
    return(data.frame(First,Last,SR_Out,MR_Out,PFX_30_A,PFX_60_A,PFX_90_A,PFX_120_A,PFX_150_A,PFX_180_A,PFX_Drag_30_A,PFX_Drag_60_A,PFX_Drag_90_A,PFX_Drag_120_A,PFX_Drag_150_A,PFX_Drag_180_A,
                            W_30_A,W_60_A,W_90_A,W_120_A,W_150_A,W_180_A,W_Drag_30_A,W_Drag_60_A,W_Drag_90_A,W_Drag_120_A,W_Drag_150_A,W_Drag_180_A,
                            PFX_30_D,PFX_60_D,PFX_90_D,PFX_120_D,PFX_150_D,PFX_180_D,PFX_Drag_30_D,PFX_Drag_60_D,PFX_Drag_90_D,PFX_Drag_120_D,PFX_Drag_150_D,PFX_Drag_180_D,
                            W_30_D,W_60_D,W_90_D,W_120_D,W_150_D,W_180_D,W_Drag_30_D,W_Drag_60_D,W_Drag_90_D,W_Drag_120_D,W_Drag_150_D,W_Drag_180_D));
}

##################################
# PITCHf/x Late Break (3 inches) #
##################################

Late_Break_PFX <- function(PITCH){
    # Set gravity
    g <- -32.174;
    # Find the time for the pitch to reach the plate in y
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the PFX movement on the pitch
    pfx_loc <- PFX_Mvt_Norm(PITCH,677/24);
    if (pfx_loc >= 3){
        t_rem <- sqrt(0.5/sqrt(PITCH$ax^2 + (PITCH$az - g)^2));
        t_break <- t_plate - t_rem;
        a_arc <- (PITCH$ax)^2 + (PITCH$ay)^2 + (PITCH$az)^2;
        b_arc <- 2*(PITCH$ax*PITCH$vx0 + PITCH$ay*PITCH$vy0 + PITCH$az*PITCH$vz0);
        c_arc <- (PITCH$vx0)^2 + (PITCH$vy0)^2 + (PITCH$vz0)^2;
        arc_rem <- Arc_Length(a_arc,b_arc,c_arc,t_break,t_plate);
    } else {
        t_rem <- 0;
        arc_rem <- 0;
    }
    LB_PFX <- list(pfx_loc = pfx_loc, t_rem = t_rem, arc_rem = arc_rem);
    return(LB_PFX);
}

####################################################
# PITCHf/x Late Break with Drag Removed (3 inches) #
####################################################

Late_Break_PFX_Drag <- function(PITCH){
    # Set gravity
    g <- -32.174;
    # Find the time for the pitch to reach the plate in y
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the PFX movement with drag removed on the pitch
    pfx_drag_loc <- PFX_Mvt_Drag_Norm(PITCH,677/24);
    if (pfx_drag_loc >= 3){
        # Solve for the time of a 3-inch displacement iteratively
        t_break <- PFX_Drag_LB_IT(PITCH,677/24,17/12);
        t_rem <- t_plate - t_break;
        a_arc <- (PITCH$ax)^2 + (PITCH$ay)^2 + (PITCH$az)^2;
        b_arc <- 2*(PITCH$ax*PITCH$vx0 + PITCH$ay*PITCH$vy0 + PITCH$az*PITCH$vz0);
        c_arc <- (PITCH$vx0)^2 + (PITCH$vy0)^2 + (PITCH$vz0)^2;
        arc_rem <- Arc_Length(a_arc,b_arc,c_arc,t_break,t_plate);
    } else {
        t_rem <- 0;
        arc_rem <- 0;
    }
    LB_PFX_Drag <- list(pfx_drag_loc = pfx_drag_loc, t_rem = t_rem, arc_rem = arc_rem);
    return(LB_PFX_Drag);
}

#####################################
# w-Direction Late Break (3 inches) #
#####################################

Late_Break_W <- function(PITCH){
    # Find the time for the pitch to reach the plate in y
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the in-plane movement on the pitch
    w_loc <- abs(W_Mvt_Norm(PITCH,677/24));
    if (w_loc >= 3){
        # Find the PITCHf/x basis
        BV <- PFX_Basis(PITCH);
        BV_T <- t(BV);
        # Find the w-acceleration of the pitch
        aw <- BV_T[3,] %*% c(PITCH$ax,PITCH$ay,PITCH$az);
        gw <- BV_T[3,] %*% c(0,0,-32.174);
        # Solve for the time of a 3-inch displacement
        t_rem <- sqrt(0.5/abs(aw - gw));
        t_break <- t_plate - t_rem;
        a_arc <- (PITCH$ax)^2 + (PITCH$ay)^2 + (PITCH$az)^2;
        b_arc <- 2*(PITCH$ax*PITCH$vx0 + PITCH$ay*PITCH$vy0 + PITCH$az*PITCH$vz0);
        c_arc <- (PITCH$vx0)^2 + (PITCH$vy0)^2 + (PITCH$vz0)^2;
        arc_rem <- Arc_Length(a_arc,b_arc,c_arc,t_break,t_plate);
    } else {
        t_rem <- 0;
        arc_rem <- 0;
    }
    LB_W <- list(w_loc = w_loc, t_rem = t_rem, arc_rem = arc_rem);
    return(LB_W);
}

#######################################################
# w-Direction Late Break with Drag Removed (3 inches) #
#######################################################

Late_Break_W_Drag <- function(PITCH){
    # Find the time for the pitch to reach the plate in y
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_plate <- PITCH$y0 - (17/12);
    t_plate <- (-b - sqrt(b^2 - 4*a*c_plate))/(2*a);
    # Find the in-plane movement without drag on the pitch
    w_drag_loc <- abs(W_Mvt_Drag_Norm(PITCH,677/24));
    if  (w_drag_loc >= 3){
        t_break <- W_Drag_LB_IT(PITCH,677/24,17/12);
        t_rem <- t_plate - t_break;
        a_arc <- (PITCH$ax)^2 + (PITCH$ay)^2 + (PITCH$az)^2;
        b_arc <- 2*(PITCH$ax*PITCH$vx0 + PITCH$ay*PITCH$vy0 + PITCH$az*PITCH$vz0);
        c_arc <- (PITCH$vx0)^2 + (PITCH$vy0)^2 + (PITCH$vz0)^2;
        arc_rem <- Arc_Length(a_arc,b_arc,c_arc,t_break,t_plate);
    } else {
        t_rem <- 0;
        arc_rem <- 0;
    }
    LB_W_Drag <- list(w_drag_loc = w_drag_loc, t_rem = t_rem, arc_rem = arc_rem);
    return(LB_W_Drag);
}

#######################
# Arc Length Function #
#######################

Arc_Func <- function(a,b,c,t){
    af <- (2*sqrt(a)*(b + 2*a*t)*sqrt(c+t*(b+a*t))-
          (b^2-4*a*c)*log(b+2*a*t+2*sqrt(a)*sqrt(c+t*(b+a*t))))/(8*a^(3/2));
    return(af);
}

################################
# Arc Length Between Two Times #
################################

Arc_Length <- function(a,b,c,t_0,t_1){
    arc_0 <- Arc_Func(a,b,c,t_0);
    arc_1 <- Arc_Func(a,b,c,t_1);
    arc <- arc_1 - arc_0;
    return(arc);
}

########################################
# Three-Inch PFX Drag Iterative Solver #
########################################

PFX_Drag_LB_IT <- function(PITCH,y_pitcher,y_plate){
    pfx_d <- PFX_Mvt_Drag_Norm(PITCH,y_pitcher);
    y_new <- y_pitcher;
    while(abs(pfx_d - 3) > 0.001){
        # Find the midpoint in y for the pitch
        y_new <- (y_pitcher + y_plate)/2;
        # Find the movement at this point
        pfx_d <- PFX_Mvt_Drag_Norm(PITCH,y_new);
        if (pfx_d > 3){
            y_pitcher <- y_new;
        } else {
            y_plate <- y_new;
        }
    }
    # Find the time associated with y
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_new;
    return((-b - sqrt(b^2 - 4*a*c_loc))/(2*a));
}

######################################
# Three-Inch W Drag Iterative Solver #
######################################

W_Drag_LB_IT <- function(PITCH,y_pitcher,y_plate){
    w_d <- abs(W_Mvt_Drag_Norm(PITCH,y_pitcher));
    y_new <- y_pitcher;
    while(abs(w_d - 3) > 0.001){
        # Find the midpoint in y for the pitch
        y_new <- (y_pitcher + y_plate)/2;
        # Find the movement at this point
        w_d <- abs(W_Mvt_Drag_Norm(PITCH,y_new));
        if (w_d > 3){
            y_pitcher <- y_new;
        } else {
            y_plate <- y_new;
        }
    }
    # Find the time associated with y
    a <- 0.5*PITCH$ay;
    b <- PITCH$vy0;
    c_loc <- PITCH$y0 - y_new;
    return((-b - sqrt(b^2 - 4*a*c_loc))/(2*a));
}