! Final Project for Fortran - Andrew Werner - March 13, 2007 - Flight of a Golf Ball ! Design a project that allows you to investigate a physics problem in depth throw the use of a Fortran program. ! The problem I chose to investigate is the flight path of a golf ball. ! Assumptions being used(made) ! *The golf shot is being made at standard ambient temperature and pressure of 25 degs C (70 F) and 100KPa. ! That the average driver clubhead mass is 200 grams. ! The drag coefficient of a golf ball is 0.21 ! The lift coefficient of a golf ball is 0.14 ! The radius of a golf ball is 21.335 mm ! The mass of a golf ball is 45.93 grams ! The coefficient of restitution of the club being used is at the USGA limit of 0.83 ! The average mass of a driver club head is 200 grams ! *The golf ball maintains a constant spin rate throughout the flight of the shot. ! *The launch angle of the shot is 3 degrees more than the loft of the driver. ! *The rotation of the golf ball is only in the x-direction and does not have any (what would be) z-components ! That the force of gravity is a constant, would actually depend upon elevation => distance from the center of the earth. !jc Good to list assumtions here. PROGRAM golf_ball IMPLICIT NONE REAL(KIND=8)::launch_angle,swing_speed,loft_club,ball_speed,delta_t=.001,lC=0.14,dC=0.21,air_density REAL(KIND=8)::mass_clubhead, drag_force, lift_force, g=9.8,radius,ball_mass REAL(KIND=8)::area_cross, swing_speed_mph, launch_angle_deg, last_point REAL,DIMENSION(100000)::a_drag,a_lift,a_x,v_x,x_m,a_y,v_y,y_m,theta,V ,y_yards,x_yards INTEGER::i INTEGER:: PGBEG, IER, max_distance_m, max_distance_yds !LOGICAL:: y_y REAL(KIND=10)::one=1, zero=0, n_one=(-1), pi !used to determine the value of pi pi = (ACOS(n_one)) !Lanuch angle = launch angle of the golf ball in reltation to the ground, figure from loft of club, in degrees !swing_speed = the speed of the particular golfer, prompt for mph then convert to m/s !loft_club = the loft of the club being used in degrees !ball_speed = the speed of the golf ball in m/s !t = time in seconds !delta_t = the increment of time, being used to create a number of heights at certain t's !lC = Lift coefficient !dC = drag coefficient !air_density = density of the air being traveled through, in g/m^3 !radius = golf ball radius in mm !ball_mass = mass of the golf ball in grams !mass_clubhead = mass of the clubhead of the golf club !area_cross = the cross sectional area of a golf ball, pi*r^2 !drag_force = the drag force felt by the golf ball !lift_force = the lift force felt by the golf ball !jcc Forces in N with conversions below ? PRINT *, 'Welcome to the driver distance figure-outer.' PRINt *, ' The driver distance figure outer will ask you for your swing speed along with the loft of your driver,' PRINT *, ' and in return figure out how far your shot will travel.' PRINT *, "So let's get started." PRINT *,' ' PRINT *, 'What is your swing speed in mph?' READ *, swing_speed_mph PRINT *, ' ' PRINT *, 'What is the loft of you driver (in degrees)?' READ *, loft_club !To decrease complications in this project, I will use the case of standard temperature and pressure (25 C and 100KPA) !At standard ambient temperature and pressure the density of air is 1.168 kg/(m^3 => 1168 g/(m^3) air_density = 1.168 !kg/m^3 !The properties of the golf ball of radius and mass will be those as defined by the USGA in the rules of golf. !The minimum diameter of a golf ball is specified as 42.67 mm, radius =>21.335 mm. radius = .021335 !m !The maximum mass of the golf ball is specified as 45.93 grams by the USGA in the rules of golf. ball_mass = .04593 !kg !The average mass of the driver clubhead is said to be 200 grams. mass_clubhead = .200 !kg !Here I am going to convert the ball_speed_mph from mph into meters/second !5280 ft = 1 mile, 3600 secs = 1 hour, .304800610 m = 1 ft swing_speed = swing_speed_mph * (5280./3600.) * .30480061 !Here is the relationship between the swing speed and the resulting ball speed. !By the conservation of momentum, we know m1v1 = m2v2 !ball_speed = (swing_speed * ball_mass) / clubhead_mass This equation would seem to be of the correct form, !However one of the hot topics in golf is the coefficient restitution, !which is a measure of the efficiency of the kinetic energy transfer between club and ball ! The USGA has also put a limit on this attribute of the golf club, which is 0.83 ! Applying this this attribute to our conservation of momentum equation we get ball_speed = swing_speed * ( (1.+0.83)/(1.0+ (ball_mass/mass_clubhead) ) ) !For the launch angle, I'm going to take into consideration that most golfers catch the ball on the slight ! upswing with their driver, adding a few degrees on top of their clubhead loft for their launch angle launch_angle_deg = loft_club + 3.0 launch_angle = (launch_angle_deg / 180) * pi !Converting the launch angle into radians in order tu use intrinsic functions. !In order to be congruent with the forms of equations I am using from my derivations I will define ! the launch angle to = theta, and the angle of (90 - theta) to = phi !theta = launch_angle !phi = 90. - theta !For the cross_sectional area of the golf ball, we get the form of area_cross = pi * (radius**2.) !For the value of pi, above I used a method from one of my previous homeworks !Now that we have the initial velocity of the ball and the launch angle, we can go to town on the ! equations that determine the flight of the ball !jc good descritption of the forces. ! The forces on the golf ball are the drag force, the lift force, and the force of gravity. ! While the force of gravity is always downward (-y), the drag force is always parallel to the direction of motion, ! while the lift force is always perpendicular to the direction of motion, and in the case of a golf ball, ! since there is always backspin imparted on the golf ball the lift force will always be upward (+y) !From these thoughts we can determine that the both the lift force and the drag force will have horizontal and !vertical components. !Using the form of F=ma, and Ftotal = Fdrag + Flift + mg = ma, we can come up with the forms of ! a = (Flift + Fdrag + mg) / ball_mass !Seperating this form of the acceleration of the golf ball into x and y components yields ! x direction: a_x = ( lift_force*COS(phi) + drag_force*SIN(theta) ) / ball_mass ! y direction: a_y = ( lift_force*SIN(phi) + drag_force*COS(theta) + ball_mass*g ) / ball_mass ! After taking two derivatives of y''=a, y'=at + v, y=(1/2)at^2 + vt + y, and the same for x, I come up with ! x = .5*a_x*(t**2) + t*(ball_speed*SIN(theta)) ! y = .5*a_y*(t**2) + t*(ball_speed*COS(theta)) !As equations of the y and x ball positions as a function of time (and currently theta) !Now for the drag and lift forces. !lift_force = .5*air_density*(ball_speed**2)*cross_area*lC !drag_force = .5*air_density*(ball_speed**2)*cross_area*dC !lift_force = .5*air_density*(ball_speed**2)*area_cross*lC !drag_force = .5*air_density*(ball_speed**2)*area_cross*dC !a_x = ( -lift_force*COS(phi) - drag_force*SIN(theta) ) / ball_mass !a_y = ( lift_force*SIN(phi) - drag_force*COS(theta) - ball_mass*g ) / ball_mass theta(1) = launch_angle !With V(i) being the magnitude of the velocity of the golf ball, we can us V = v_x^2 + v_y^2 for the magnitude v_x(1) = ball_speed * COS(launch_angle) v_y(1) = ball_speed * SIN(launch_angle) V(1) = SQRT( v_x(1)**2. + v_y(1)**2. ) a_x(1) = 0 a_y(1) = 0 x_m(1) = 0 y_m(1) = 0 a_lift(1) = 0 a_drag(1) = 0 !PRINT *, swing_speed, ball_speed !PRINT *, a_y(1:5) !PRINT *, v_y(1:5) !PRINT *, y(1:5) !PRINT *,' ' !PRINT *, air_density, area_cross, dC, ball_mass, lC, V(1) DO i=2,100000 a_drag(i) = ( (.5*air_density*area_cross*dC)/ball_mass) * (V(i-1)**2.) a_lift(i) = ( (.5*air_density*area_cross*lC)/ball_mass) * (V(i-1)**2.) !jc There is a problem here - you take -COS(theta(i)), but !jc theta(i) it evaluated below - this means that you probably !jc always get from this COS. The SIN at the end seems to be !jc taken correctly. a_x(i) = a_drag(i)*(-COS(theta(i))) + a_lift(i)*(-SIN(theta(i-1))) v_x(i) = v_x(i-1) + a_x(i)*delta_t x_m(i) = x_m(i-1) + v_x(I)*delta_t x_yards(i) = x_m(i)*(1./0.9144) !jc The same problem with theta shows up here as does above. !jc The consequence seems to be that there is alway a drag force !jc in the x-direction, regrdless of the angle of flight. a_y(i) = -g + a_drag(i)*(-SIN(theta(i))) + a_lift(i)*(COS(theta(i-1))) v_y(i) = v_y(i-1) + a_y(i)*delta_t y_m(i) = y_m(i-1) + v_y(I)*delta_t y_yards(i) = y_m(i)*(1./0.9144) !There are 0.9144 meters in a yard. V(i) = SQRT( v_x(i)**2. + v_y(i)**2. ) theta(i) = ATAN( v_y(i) / v_x(i) ) last_point = 0 IF ( (y_m(i)<0) .AND. (last_point == 0)) THEN max_distance_m = x_m(i) EXIT END IF END DO !PRINT *, a_drag(1:5) Used for debugging the program. !PRINT *, a_lift(1:5) !PRINT *, a_x(1:5) !PRINT *, v_x(1:5) !PRINT *, x(1:5) !IER = PGBEG(0,'ballflight.png/PNG',1,1) !IER = PGBEG(0,'ballflight.ps/PS',1,1) IER = PGBEG(0,'/xserve',1,1) !Naming the graph, selecting the graphics device. IF (IER.NE.1) STOP CALL PGENV(0,MAXVAL(REAL(x_yards(1:100000))),0,2*MAXVAL(REAL(y_yards(1:10000))),0,1) !Defining plot scales according to the min and max array values. CALL PGLAB('Distance (yds)','Height (yds)','Flight Path of a Golf Ball') !Creating labels for the graph CALL PGPT(20000, REAL(x_yards(1:20000)), REAL(y_yards(1:20000)),20) !Creating the point markers in the graph. CALL PGEND max_distance_yds = max_distance_m * (1./ 0.9144) !There are 0.9144 meters in a yard. PRINT *,' ' PRINT *, 'The length of you shot was', max_distance_yds,'yards' PRINT *, 'Wow, great shot!' PRINT *, ' ' PRINT *, 'Have a great day!' END PROGRAM golf_ball !jc _Very_ nicely commented. The problem with theta(i) seems to invalidate !jc your results, but it would be easy enough to fix.