PROGRAM long_ball
IMPLICIT NONE

  ! This program will simulate the flight of a baseball after impact with a bat
  ! launch_v is the launch speed off the bat in m/s
  ! bat_v is the bat speed in mph
  ! pitch_v is the speed of the pitch in mph
  ! gravity is in m/s^2
  ! delta_t is the time increment in seconds
  ! i_angle is the initial trajectory angle in degrees
  ! ball_radius is the radius of the ball in m
  ! mall_mass is the mass of the ball in kg
  ! i_rad is the intitial trajectory angle in radians
  ! max_distance is the maximum distance the baseball travels in feet
  ! zk is the distance from the knob of the bat in m
  ! cm is the center of mass of the bat in m
  ! rg is the radius of gyration
  ! roe is the atmospheric pressure in kg/m^3
  ! bat_si is the bat speed in m/s
  ! pitch_si is the pitch speed in m/s
  ! Written 014-Feb-2007
  ! height and distance are in m
  ! the x and y velocities are in m/s
  ! angle is the angle of the bat in reference to the ground in radians
  ! spin is the spin of the ball in rpm
  ! distance_ft is the distance the ball travels in feet 
  ! height_ft is the height the ball travels in feet
  ! written 3/10/07
  ! by Lance Wheeler

!jc  Good inclusion of definitions of units and variables, though it would 
!jc  probably look nicer if it were separated from the rest of the header
!jc  comment.

REAL :: launch_v, bat_v, pitch_v, gravity, delta_t = .01, i_angle, lift_const, drag_const
REAL :: R, area, pi = 3.14159265, ball_radius = .0366, ball_mass = .145, bat_mass = .9, i_rad, max_distance
REAL :: zk = .564, cm = .564, cor = .6, rg = .23, roe =1.21, drag_coef = .5, bat_si, pitch_si 
REAL,DIMENSION(10000) :: height, distance, velocity_x, velocity_y, angle, v_total, spin, S, lift_coef
REAL,DIMENSION(10000) :: distance_ft, height_ft 
INTEGER :: i, j

PRINT *, 'Welcome to Home Run Derby'

PRINT *, 'Enter the speed of the pitch (20 to 105 mph)'
READ (*,*) pitch_v                                                       ! 1st user input speed of the pitch

 pitch_si = pitch_v * (1609.344 / 3600.)                                 ! converts input from mph to m/s

PRINT *, 'Enter the bat speed (20 to 120 mph)'                           ! 2nd user input is bat speed
READ (*,*) bat_v

 bat_si = bat_v * (1609.344 / 3600.)                                     ! converts input from mph to m/s

PRINT *, 'Enter the uppercut angle of the swing (10 to 70 degrees)'      ! 3rd user input is trajectory
READ (*,*) i_angle

!jc  So this is ball's initial angle?  You probably should have called it 
!jc  something else - it isn't really the angle of the swing of the bat
!jc  which is what this sounds like.

 i_rad = i_angle * (pi / 180.)                                           ! converts input from degrees to radians
  area = pi * (ball_radius ** 2.)                                        ! area of baseball

  R = (ball_mass/bat_mass) * ( 1.+((zk - (cm**2.)) / rg))                ! reynolds number

!jc  Where did this come from?  References in code are nice.
  launch_v = (((cor - R) / (1. + R)) * pitch_si) + (((cor + 1.) / (1. + R)) * bat_si)     ! speed of ball off bat

  lift_const = roe * area                                      ! calculates the constants in the force of lift

  drag_const = roe * area * drag_coef                          ! calculates the constants in the force of drag

  gravity = -(ball_mass * 9.81)                                ! calculates the force of gravity

  velocity_x(1) = (launch_v*(COS(i_rad)))                      ! sets initial x component of velocity
  velocity_y(1) = launch_v*(SIN(i_rad))                        ! sets initial y component of velocity
  v_total(1) = launch_v                                        ! sets initial velocity
  angle(1) = i_rad                                             ! sets initial trajectory angle
  distance(1) = 0                                              ! sets initial distance
  height(1) = 1.2192                                           ! height of bat from ground during swing
  spin(1) = 500                                                ! sets initial spin in rpm
  !jc  Why rpm ?
  lift_coef(1) = 0                                             ! sets intial amount of lift force

DO i = 2, 10000                                         ! loop to calculate position of the baseball in flight

 spin(i) = spin(i-1) - ((1./5.) * delta_t * spin(i-1))            ! determines the decay of spin during flight
 S(i) = (ball_radius * spin(i)) / velocity_y(1)                   ! variable used to determine value of lift coefficient

!jc What does this do?  Why the variate in lift_coef ?
  IF (S(i) <= .1) THEN                               ! statement gives appropriate values to lift coeffient
   lift_coef(i) = 1.5 * S(i)
  ELSE
   lift_coef(i) = .09 + .6 * S(i)
  END IF

 velocity_x(i)=(velocity_x(i-1) + (delta_t*(((-lift_const*lift_coef(i)*(v_total(i-1)**2.))*SIN(angle(i-1))) &
               +((-drag_const*(v_total(i-1)**2.))*COS(angle(i-1)))) / (2.*ball_mass)))        ! x velocity component

 velocity_y(i)=(velocity_y(i-1)+(delta_t*((((lift_const*lift_coef(i)*(v_total(i-1)**2.))/2.)*COS(angle(i-1))) &
               +(((-drag_const*(v_total(i-1)**2.))/2.)*SIN(angle(i-1)))+gravity) / (ball_mass)))  ! y velocity component

 v_total(i) = SQRT((velocity_x(i)**2.) + (velocity_y(i)**2.))          ! velocity change after launch

 angle(i) = ATAN((velocity_y(i)) / (velocity_x(i)))                    ! changing angle after launch

 distance(i) = velocity_x(i)*delta_t + distance(i-1)                   ! x position of baseball
 distance_ft(i) = distance(i) * (1. / .305)                            ! x position converted from meters to feet

 height(i) = velocity_y(i)*delta_t + height(i-1)                       ! y position of baseball
 height_ft(i) = height(i) * (1. / .305)                                ! y position converted from meters to feet

  IF (height(i) <= 0) THEN                                     ! stops loop when ball hits the ground
   max_distance = distance_ft(i-1)                             ! determines max distance when ball hits ground
   EXIT
  END IF

END DO

 WRITE(*,'(A20, x ,F4.0, x, A4)') 'You hit the baseball', max_distance,'feet'       ! prints max distance of baseball

 CALL flight_plot(distance_ft, height_ft, i)                     ! calls flight_plot subroutine

END PROGRAM long_ball

SUBROUTINE flight_plot(distance_ft, height_ft, i)                ! subroutine to plot the flight of a baseball
IMPLICIT NONE

INTEGER,INTENT(IN) :: i
REAL,DIMENSION(1000),INTENT(IN) :: distance_ft, height_ft
INTEGER :: j, ier, pgbeg
REAL,DIMENSION(11) :: fence_y = (/ (j, j=2, 12) /)
REAL,DIMENSION(11) :: fence_x = 385

   ier = pgbeg(0, '/xserve',1,1)                                      ! Initializes plot

  CALL pgsci(7)                                                       ! sets color to yellow
  CALL pgenv(0., 550., 0., 400., 0, 1)                                ! Sets scale for plot
  CALL pglab('distance(feet)','height(feet)','Home Run Derby')        ! labels axes and title of the plot
  CALL pgsci(2)                                                       ! changes color to red
  CALL pgpt(i, REAL(distance_ft(1:i)), REAL(height_ft(1:i)), 1)	      ! Plot the points
  CALL pgtext(-40., 10., 'Batter')                                    ! labels where the batter stands
  CALL pgsci(10)                                                      ! changes color to green
  CALL pgtext(353., -10., 'Fence')                                    ! labels the fence
  CALL Pgpt(11, fence_x, fence_y, 6)                                  ! plots the fence
  CALL pgline(11, fence_x, fence_y)                                   ! makes line for points on fence
   IF (distance_ft(i) > 385.) THEN                                    ! determines if ball went home run distance
    CALL pgsci(1)                                                     ! changes color to white
    CALL pgsch (5.)                                                   ! enlarges font by 5x
    CALL pgtext(100., 300., 'HOME RUN!')                              ! prints if ball went 385 feet
   ELSE 
    CALL pgsci(1)                                                     ! changes color to white
    CALL pgsch (5.)                                                   ! enlarges font by 5x
    CALL pgtext(200., 300., 'OUT!')                                   ! prints if ball did not go home run distance
   END IF
  CALL pgend                                                          ! closes the plot

!jc  Fancy plot
END SUBROUTINE flight_plot

!jc  Nice program.  There are some places where the comments should be more 
!jc  descriptive, so it is a little more clear where things are coming from.
