c<html>
c<body>
c<pre>
program fall
implicit none
c
c Program to calculate the dynamics of a falling body
c
c Arrays:
c v - velocity at each time integration step
c z - height at each time integration step
c t - time for each corresponding v and z
c zreal - Actual height at time t(i) for comparison with computed z
c
c John Mahaffy 4/15/95
c
integer np
parameter (np=2000)
c
c In this program, I am using blank common just to save space in the
c executable program file (a.out). It only appears in the main program
c Named common "const" communicates between subroutines. By using the
c "common/const/..." statement in the main routine, I insure that some
c machines do not lose values in that common block when the first
c subroutine containing "common/const/..." is exited.
c
c<a name="com"><font color="FF0000">
common v(np),z(np),t(np), zreal(np)
c</font></a>
real v,z,t,zreal
common/const/dtmin,g
real dtmin,g
real dt
integer nsteps
c
call input(z,dt)
call odesolve(v,z,t,dt,np,nsteps)
call realans(t,z,nsteps,zreal)
call output (t,z,zreal,v,nsteps)
stop
end
c
c<a name="bdata"><font color="FF0000">
block data cons
c</font></a>
c
c Initialize various constants
c
c John Mahaffy 4/15/95
c
c g - gravitational acceleration
c dtmin - time step selected if 0.0 is entered
c
common/const/dtmin,g
real dtmin,g
data dtmin,g/.001,9.807/
c<a name="ebd"><font color="FF0000">
end block data cons
c</font></a>
subroutine input (z,dt)
c
c Obtain user input for initial height and solution time step
c
c John Mahaffy 4/15/95
c
c Output Arguments:
c z(1) - initial height
c dt - integration time step
c
real z(*)
common/const/dtmin,g
real dtmin,g
c
write(6,*) ' What is the initial height (meters)?'
read *, z(1)
write(6,*) ' What time step size do you want (seconds)?'
read *, dt
if(dt.le.0.) dt=dtmin
return
end
c
subroutine odesolve(v,z,t,dt,np,nsteps)
c
c Solve the Ordinary Differential Equation of motion for the body
c
c John Mahaffy 4/15/95
c
c Arguments:
c Input
c dt - timestep size
c np - size of v,z, and t arrays
c Output:
c v - velocity
c z - height
c t - time
c nsteps - last step in the integration
c
real v(*),z(*),t(*),dt,rdt
common/const/dtmin,g
real dtmin,g
c
c Solve the equation for a falling body
c
c 2
c d z
c ---2 = - g
c d t
c
c NOTE: When you try this you will see that solving this for of the
c equation will result in the amplification of machine round-off
c error. Never use this approach for this type equation.
c
c Set remaining initial conditions:
c
t(1)=0.
v(1)=0.
c
c For the first step we need to fold in the initial velocity condition
c v(1)=(z(2)-z(0))/(2*dt) giving z(0)=z(1) - 2*dt*v(1)
c Substituting this value of z(0) into the equation at t=0 gives:
c
z(2) = z(1) - .5*g*dt**2
t(2) = dt
c
c Now loop through time steps until z goes negative or we run out of space
c
rdt=.5/dt
do 100 i=2,np-1
z(i+1)= 2*z(i)-z(i-1)-dt**2*g
v(i)=(z(i+1) - z(i-1)) * rdt
t(i+1)=t(i)+dt
if(z(i).lt.0.) go to 200
100 continue
write(6,*) 'Ran out of space to continue integration'
write(6,*) ' Last height was ',z(np),' meters'
i=np
200 nsteps=i
c return
end
c
subroutine realans(t,z,nsteps,zreal)
c
c Get the values of the analytic solution to the differential equation
c for each time point to check the numerical accuracy.
c
c John Mahaffy 4/15/95
c
common/const/dtmin,g
real dtmin,g
real t(*),z(*),zreal(*)
c
do 10 i=1,nsteps
10 zreal(i)=z(1)-.5*g*t(i)**2
return
end
c
subroutine output(t,z,zreal,v,nsteps)
implicit none
c
c Outputs the full results of the time integration
c
c John Mahaffy 4/15/95
c
real v(*),z(*),t(*), zreal(*)
integer nsteps,i
print *, 'Look in the file fall.output for detailed results'
open (11,file='fall.output')
write(11,2000)
do 300 i=1,nsteps
write(11,2001) t(i),v(i),z(i),zreal(i)
300 continue
2000 format (33x,'Computed',8x,'Actual',/,
& 6x,'Time',9x,'Velocity', 8x,'Height',8x,'Height')
2001 format (1x,1p,4e15.7)
return
c<a name="1"><font color="FF0000">
end subroutine output
c</font></a>
c</pre>
c</body>
c</html>
Copyright ©2017 - STI - Todos os direitos reservados