转载请标明是引用于 http://blog.csdn.net/chenyujing1234
源码:
http://www.rayfile.com/zh-cn/files/e5f02f0a-8799-11e1-b6a2-0015c55db73d/
高尔夫球模型的算法是用Fortan 90写的，以前从来没接触过！！！！
FORTRAN，亦译为福传，是英文“FORmula TRANslator”的缩写，译为“公式翻译器”，它是世界上最早出现的计算机高级程序设计语言，广泛应用于科学和工程计算领域。FORTRAN语言以其特有的功能在数值、科学和工程计算领域发挥着重要作用。
从现在开始学习吧。
Fortran语言的最大特性是接近数学公式的自然描述，在计算机里具有很高的执行效率。 易学，语法严谨。 可以直接对矩阵和复数进行运算，这一点类似matlab。 自诞生以来广泛地应用于数值计算领域，积累了大量高效而可靠的源程序。 很多专用的大型数值运算计算机针对Fortran做了优化。 广泛地应用于并行计算和高性能计算领域。 Fortran90，Fortran95，Fortran2003的相继推出使Fortran语言具备了现代高级编程语言的一些特性。
Fortran 90 不是软件，是语言规范。
Compaq Visual Fortran 是国内常用的支持 Fortran 90 的编译器.
下载地址:
http://dc.whu.edu.cn/bencandy.php?fid=8&id=177
安装成功后的Visual Fortran与VC6.0的界面很像
=============================================================================================
由于要在Fotran中要画图，究竟怎么实现呢？
参考 http://hi.baidu.com/desiree100/blog/item/e769a7b78db046f531add19b.html
Visual Fortran中有提供绘图功能。Visual Fortran的绘图功能不完全是以扩充函数的类型存在，使用它的绘图功能必须在选择Project类型时，选择Standard Graphics或QuickWin模式。
Visual Fortran的QuickWin及Standard Graphics模式在简单的绘图使用上会比较方便，它的绘图函数功能比较多样，不过效率会比较差，而且不支持动画功能。
Standard Graphics和QuickWin模式在绘图方面的使用方法完全相同，它们都是调用相同的函数来绘图。差别在于Standard Graphics只能打开一个窗口来绘图、QuickWin模式则可以打开多个窗口来绘图。QuickWin模式下可以有菜单及对话窗的功能，Standard Graphics则不行。Standard Graphics模式的程序代码可以原封不动直接转换到QuickWin模式下使用，但是QuickWin的程序代码并不一定可以直接拿到Standard Graphics模式下使用。
我试着运行一个画y = sin(x)的工程。
新建一个工程->选择Fortran Standard Graphics or QuickWin Application。
注意: 工程路径不能含有中文。
新建sinGraph.f90文件
! sin函数的绘图范例
program Plot_Sine
use DFLIB
implicit none
integer, parameter :: lines=500 ! 用多少线段来画函数曲线
real(kind=8), parameter :: X_Start=-5.0 ! x轴最小范围
real(kind=8), parameter :: X_End=5.0 ! x轴最大范围
real(kind=8), parameter :: Y_Top=2.0 ! y轴最大范围
real(kind=8), parameter :: Y_Bottom=-2.0 ! y轴最小范围
integer :: result ! 取回绘图函数运行状态
integer(kind=2) :: color ! 设定颜色用
real(kind=8) :: step ! 循环的增量
real(kind=8) :: x,y ! 绘图时使用,每条小线段都连接
real(kind=8) :: NewX,NewY ! (x,y)及(NewX,NewY)
real(kind=8), external :: f ! 待绘图的函数
type(wxycoord) :: wt ! 返回上一次的虚拟坐标位置
type(xycoord) :: t ! 返回上一次的实际坐标位置
! 设定虚拟坐标范围大小
result=SetWindow( .true. , X_Start, Y_Top, X_End, Y_Bottom )
! 用索引值的方法来设定颜色
result=SetColor(2) ! 内定的2号是应该是绿色
call MoveTo(10,20,t) ! 移动画笔到窗口的(10,20)
call OutGText("f(x)=sin(x)") ! 写出内容
! 使用全彩RGB 0-255的256种色阶来设定颜色
color=RGBToInteger(255,0,0) ! 把控制RGB的三个值转换到color中
result=SetColorRGB(color) ! 利用color来设定颜色
call MoveTo_W(X_Start,0.0_8,wt) ! 画X轴
result=LineTo_W(X_End,0.0_8) !
call MoveTo_W(0.0_8,Y_Top,wt) ! 画Y轴
result=LineTo_W(0.0_8,Y_Bottom) !
step=(X_End-X_Start)/lines ! 计算小线段间的X间距
! 参数#FF0000是使用16进制的方法来表示一个整数
result=SetColorRGB(#FF0000)
! 开始绘制小线段们
do x=X_Start,X_End-step,step
y=f(x) ! 线段的左端点
NewX=x+step
NewY=f(NewX) ! 线段的右端点
call MoveTo_W(x,y,wt)
result=LineTo_W(NewX,NewY)
end do
! 设定程序结束后,窗口会继续保留
result=SetExitQQ(QWIN$EXITPERSIST)
end
! 所要绘图的函数
real(kind=8) function f(x)
implicit none
real(kind=8) :: x
f=sin(x)
return
end function f
马上就编译通过了.
运行结果为
============================================================================================================================
试了一个Demo后现在该上演自己的golfModem代码了.
! 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 degrees C (70 F) and 100KPa.
! That the average driver club head 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.
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=8)::one=1,zero=0,n_one=(-1),pi
!used to determine the value of pi
pi = (ACOS(n_one))
!Launch angle = launch angle of the golf ball in relation 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 club head 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
PRINT *, 'Welcome to the driver distance figure-outer'
PRINt *, ' The driver distance figure outer will ask you for your swing speedalong 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 club head 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 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 club head 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 to 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
! 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
!Separating 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.)
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)
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.ps/PS',1,1)
IER = PGBEG(0,'/xserve',1,1) !Naming the graph,选择graphics设备
IF (IER.NE.1) STOP
CALL PGENV(0,MAXVAL(REAL(x_yards(1:100000))),0,2*MAXVAL(REAL(y_yards(1:10000))),0,1)
!根据最大最小的数据组定义plot的scales
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)
!在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
编译后提示如下错误:
Linking...
GolfModem.obj : error LNK2001: unresolved external symbol _PGBEG@20
GolfModem.obj : error LNK2001: unresolved external symbol _PGENV@24
GolfModem.obj : error LNK2001: unresolved external symbol _PGLAB@24
GolfModem.obj : error LNK2001: unresolved external symbol _PGPT@16
GolfModem.obj : error LNK2001: unresolved external symbol _PGEND@0
Debug/GolfModem.exe : fatal error LNK1120: 5 unresolved externals
Error executing link.exe.
这是因为代码里的红色部分。它们用OPENGL库绘图，我没有OpenGL库，接下来准备改Visual Fortran自带的Visual Fortran的QuickWin及Standard Graphics来画图了。
以下是修改后的代码:红色为我添加的内容。
! 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 degrees C (70 F) and 100KPa.
! That the average driver club head 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.
PROGRAM golf_ball
use DFLIB
IMPLICIT NONE
type(wxycoord) :: wt ! 返回上一次的虚拟坐标位置
real(kind=8) :: x,y ! 绘图时使用,每条小线段都连接
real(kind=8) :: NewX,NewY ! (x,y)及(NewX,NewY)
integer :: result ! 取回绘图函数运行状态
real(kind=8), parameter :: X_Start=0 ! x轴最小范围
real(kind=8), parameter :: X_End=400 ! x轴最大范围
real(kind=8), parameter :: Y_Top=100 ! y轴最大范围
real(kind=8), parameter :: Y_Bottom=0 ! y轴最小范围
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=8)::one=1,zero=0,n_one=(-1),pi
!used to determine the value of pi
pi = (ACOS(n_one))
!Launch angle = launch angle of the golf ball in relation 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 club head 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
PRINT *, 'Welcome to the driver distance figure-outer'
PRINt *, ' The driver distance figure outer will ask you for your swing speedalong 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 club head 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 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 club head 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 to 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
! 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
!Separating 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)
! 设定虚拟坐标范围大小
result=SetWindow( .true. , X_Start, Y_Bottom, X_End, Y_Top )
call MoveTo_W(X_Start,0.0_8,wt) ! 画X轴
result=LineTo_W(X_End,0.0_8) !
call MoveTo_W(0.0_8,Y_Top,wt) ! 画Y轴
result=LineTo_W(0.0_8,Y_Bottom) !
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.)
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)
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
x = x_yards(i-1)
y = y_yards(i-1)
NewX = x_yards(i)
NewY = y_yards(i)
!PRINT *, 'now x is ', x,'y is',y
call MoveTo_W(x,y,wt)
result=LineTo_W(NewX,NewY)
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.ps/PS',1,1)
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!'
! 设定程序结束后,窗口会继续保留
result=SetExitQQ(QWIN$EXITPERSIST)
END PROGRAM golf_ball
