You thought we would never get here, but we are now ready to do the simulation itself.

We will assign initial conditions to the pumpkin: height 100,000 feet, velocity 0. We will figure the force acting on the pumpkin: gravity will pull it down, air resistance will push it up. We will use this to figure the acceleration on the pumpkin. Then we will assume the forces will be constant for one-half second. We will figure the effect of the pumpkin’s velocity and acceleration on its position, and the effect of its acceleration on the velocity. Then, if the pumpkin has not yet impacted the ground, we will do it all again — until the pumpkin does impact the ground.

We can pseudocode this fairly easily:

h := 30480; t := 0; v := 0; Δt := 0.5;whileh > 0dobeginf := drag (h, v) – 53.38; a := f / m; h_{new}:= .5 * a * Δt^{2}+ v * Δt + h; v_{new}:= a * Δt + v; t := t + Δt; v := v_{new}; h := h_{new}end; { report t, v }

This is cleverly written so all the acceleration and velocity units are positive when directed upwards, and negative when directed downward. This is why we can compute the new position without making changes to the signs in the formulas.

The program reports the time taken for the pumpkin to fall, and its impact velocity. But we have one unknown: m, the mass of the pumpkin. The pumpkin weighs 12 pounds in a standard gravitational field, so its mass is about 5.443 kilograms.

Unfortunately, there are more questions to answer. We are also supposed to find the height and time after dropping of the maximum velocity. We can do this with a test in the middle of the loop, changing our pseudocode to this:

h := 30480; v := 0; t := 0; Δt := 0.5; v_{max}:= 0; h_{vmax}:= 0; t_{vmax}:= 0;whileh > 0dobeginf := drag (h, v) – 53.38; a := f / 5.443; h_{new}:= .5 * a * Δt^{2}+ v * Δt + h; v_{new}:= a * Δt + v; t := t + Δt; v := v_{new}; h := h_{new};if|v| > v_{max}thenbeginv_{max}:= |v|; h_{vmax}:= h; t_{vmax}:= tendend; { report t, v, v_{max}, h_{vmax}, t_{vmax}}

Note how we took the absolute value of the velocity: the velocity is directed downwards, and so is negative. Presumably we want the peak downward velocity, rather than the peak upward velocity (which will be 0 — unless, of course, the pumpkin bounces.)

We are also supposed to find the height, time after dropping, and velocity at which the pumpkin has maximum upwards acceleration. Another test and change to the pseudocode.

h := 30480; v := 0; t := 0; Δt := 0.5; v_{max}:= 0; h_{vmax}:= 0; t_{vmax}:= 0; a_{max}:= 0; h_{amax}:= 0; t_{amax}:= 0; v_{amax}:= 0;whileh > 0dobeginf := drag (h, v) – 53.38; a := f / 5.443;ifa > a_{max}thenbegina_{max}:= a; h_{amax}:= h; t_{amax}:= t; v_{amax}:= vend; h_{new}:= .5 * a * Δt^{2}+ v * Δt + h; v_{new}:= a * Δt + v; t := t + Δt; v := v_{new}; h := h_{new};if|v| > v_{max}thenbeginv_{max}:= v; h_{vmax}:= h; t_{vmax}:= tendend; { report t, v, v_{max}, h_{vmax}, t_{vmax}, a_{max}, h_{amax}, t_{amax}, v_{amax}}

And, finally, there is the question of whether the pumpkin exceeds the speed of sound, and time after dropping, height, and velocity of the first such instance. Another test and change to the pseudocode.

h := 30480; v := 0; t := 0; Δt := 0.5; v_{max}:= 0; h_{vmax}:= 0; t_{vmax}:= 0; a_{max}:= 0; h_{amax}:= 0; t_{amax}:= 0; v_{amax}:= 0; h_{hyp}:= 0; t_{hyp}:= 0; v_{hyp}:= 0;whileh > 0dobeginf := drag (h, v) – 53.38; a := f / 5.443;ifa > a_{max}thenbegina_{max}:= a; h_{amax}:= h; t_{amax}:= t; v_{amax}:= vend; h_{new}:= .5 * a * Δt^{2}+ v * Δt + h; v_{new}:= a * Δt + v; t := t + Δt; v := v_{new}; h := h_{new};if|v| > v_{max}thenbeginv_{max}:= v; h_{vmax}:= h; t_{vmax}:= tend;ifv_{hyp}= 0thenbeginv_{snd}:= vsound (h);if|v| > v_{snd}thenbeginv_{hyp}:= v; t_{hyp}:= t; h_{hyp}:= hendendend; { report t, v, v_{max}, h_{vmax}, t_{vmax}, a_{max}, h_{amax}, t_{amax}, v_{amax}, t_{hyp}, v_{hyp}, h_{hyp}}

This seems to address all the needs.

We can code this up, now, and also graph the simulated velocity as a function of simulated height. This display will also let us know if there is anything very strange going on.

Before coding, we note that we have to keep track of a great many variables. It would be wise to draw up a little equivalence table showing where we are putting all these variables in the calculator variables. So we draw up the following list:

Var |
Use |
Var |
Use |
Var |
Use |
---|---|---|---|---|---|

A |
Subroutines |
K |
v |
U |
(unused) |

B |
Subroutines |
L |
h |
V |
(unused) |

C |
Subroutines |
M |
t |
W |
(unused) |

D |
Subroutines |
N |
a |
X |
(unused) |

E |
t |
O |
h |
Y |
(unused) |

F |
h |
P |
t |
Z |
(unused) |

G |
v |
Q |
v |
r |
(unused) |

H |
Δt |
R |
h |
θ |
(unused) |

I |
f, a, v |
S |
t |
||

J |
h |
T |
v |

Now we translate the pseudocode into a calculator program, being very careful to translate the variables in the pseudocode based on the table above.

Program DROP (456 bytes):

; Program: DROP ; ; Simulates the dropping of a 12 pound, 1 foot diameter pumpkin from ; 100,000 feet. ; ; Variables: ; A-D are used for communication with subprograms ; E is is the elapsed time since the pumpkin was dropped ; F is the current height of the pumpkin (meters) ; G is the current velocity of the pumpkin (meters/second) ; H is the delta t for the simulation (seconds) ; I is the force upon the pumpkin in Newtons, the accelleration of ; the pumpkin in m/sec^2, the velocity after the period of time, ; and the velocity of sound ; J is the height after the period of time ; K is the maximum velocity yet attained (m/s) ; L is the height at which the maximum velocity was attained (m) ; M is the time at which the maximum velocity was attained (s) ; N is the maximum upwards acceleration of the pumpkin (m/s^2) ; O is the height at which the maximum upwards acceleration was ; achieved (m) ; P is the time at which the maximum upwards acceleration was ; achieved (s) ; Q is the velocity at which the maximum upwards acceleration was ; achieved (m/s) ; R is the height at which hypersonic speed was first attained (m) ; S is the time at which hypersonic speed was first attained (s) ; T is the velocity at which hypersonic speed was first attained (m/s) ; ; Symbols: ; -> is the assignment arrow ; / is the division operator ; x^2 is the x squared operator ; <> is the not equals relational ; _ is the display triangle ViewWindow 0,100000,10000,0,800,100 ; Set the view window for the graph 0->A~Z ; Zero all variables 30480->F ; Initialize the height to 100,000 feet .5->H ; Initialize delta t to 0.5 seconds While F>0 ; Loop while the height is positive F->A ; Get the height into A G->B ; Get the velocity into B Prog "DRAG" ; Compute the drag into A A-53.38->I ; Compute the net force into I I/5.443->I ; Compute the net acceleration into I If I>N ; Test whether acceleration is greatest seen Then I->N ; If so, save acceleration into N F->O ; Save height into O E->P ; Save time into P G->Q ; Save velocity into Q IfEnd ; [End of test on greatest acceleration] .5IHx^2+GH+F->J ; Compute new height into J IH+G->I ; Compute new velocity into I E+H->E ; Update time in E I->G ; Update velocity in G J->F ; Update height in F If Abs G>K ; See if velocity is greatest seen Then Abs G->K ; If so, save velocity into K F->L ; Save height into L E->M ; Save time into M IfEnd ; [End of test on greatest velocity] If T=0 ; Test if hypersonic speed already achieved Then F->A ; If not, get height into A Prog "VSOUNDH" ; Compute speed of sound into A If Abs G>A ; Test if current velocity is hypersonic Then F->R ; If so, save height in R E->S ; Save time in S Abs G->T ; Save velocity in T IfEnd ; [End of test if hypersonic] IfEnd ; [End of test if hypersonic already achieved] Plot 3.281F,2.237Abs G ; Plot height and velocity If E<>0 ; If this is not the first point Then Line ; then connect the dots IfEnd ; [End of test of first point] WhileEnd ; [End of loop: splat!] "TIME SEC" ; Report time to impact E_ ; in seconds "V MPH" ; Report velocity at impact 2.237G_ ; in miles/hour "VMAX" ; Report maximum velocity 2.237K_ ; in miles/hour "HVMAX FT" ; Report height of maximum velocity 3.281L_ ; in feet "TVMAX" ; Report time until maximum velocity M_ ; in seconds "AMAX MPHPS" ; Report maximum acceleration 2.2447N_ ; in MPH/sec "HAMAX" ; Report height of maximum acceleration 3.281O_ ; in feet "TAMAX" ; Report time until maximum acceleration P_ ; in seconds "VAMAX" ; Report velocity at maximum acceleration 2.237Q_ ; in miles/hour "THYP" ; Report time until hypersonic S_ ; in seconds "HHYP" ; Report height at hypersonic 3.281R_ ; in feet "VHYP" ; Report velocity at hypersonic 2.237T_ ; in miles/hour "DONE" ; Report completion

Running the program gives us the following results:

Item |
Value |
---|---|

Time until impact |
184.5 sec |

Velocity at impact |
-184.8 MPH |

Maximum velocity achieved |
745.4 MPH |

Height of maximum velocity |
69,182 ft |

Time after dropping of maximum velocity |
46.5 sec |

Maximum upwards acceleration |
12.6 MPH/sec |

Height of maximum upwards acceleration |
50,959 ft |

Time after dropping of maximum upwards acceleration |
64.5 sec |

Velocity at maximum upwards acceleration |
597.0 MPH |

Time after dropping of first exceeding the speed of sound |
34 seconds |

Height at first exceeding the speed of sound |
82,327 ft |

Velocity at first exceeding the speed of sound |
661.4 MPH |

Note that the maximum upwards acceleration is 12.6 MPH/sec, over one-half a standard gravity. This is a net acceleration, so the one-half a standard gravity is in addition to the one gravity needed to balance the Earth’s gravitational field. The pumpkin is experiencing 1.6 gravities of drag. That is some air resistance!

And, as a bonus, we get the following graph of velocity as a function of height. Remember the pumpkin starts on the right of the graph, at 100,000 feet, and drops to lower altitudes, at the left of the graph:

Note that the pumpkin impacts the ground at roughly terminal velocity, about 185 mile per hour, even though it achieved a top speed of 745 miles per hour.

Changing the ViewWindow and Plot statements, and running the program again, gives us these graphs of the pumpkin’s height, velocity, and acceleration against the time since the pumpkin was dropped.

One final screen shot combines the graphs from the terminal velocity program (the second unit test of the drag function), extended to compute the terminal velocity for all heights up to 100,000 feet; the speed of sound unit test; and the simulation. The graphs show the region in which the pumpkin exceeds the speed of sound, and how it approaches, overshoots, then settles down to the terminal velocity. The color capability of the CFX-9850G is used to distinguish the various curves.

All the programs in this example are available here on these pages, or in a text file with .CAT file contents.

We can translate the pseudocode above directly into Visual Basic for Applications, choosing meaningful variable names:

Private Sub Drop() ' FUNCTIONAL DESCRIPTION: ' ' The pumpkin drop simulation. ' ' Simulates the forces acting upon a pumpkin initially released ' with zero velocity at 100,000 feet of altitude. Assumes the ' forces will be constant for half a second, computes the effect ' upon the pumpkin of the forces acting for half a second, and ' recomputes the forces. Terminates when the pumpkin's height ' becomes non-positive. ' ' Graphs the pumpkin's velocity against its height. ' ' Accumulates a number of observations of interest and reports ' them upon completion. ' ' FORMAL PARAMETERS: ' ' None. ' ' IMPLICIT INPUTS: ' ' None. ' ' IMPLICIT OUTPUTS: ' ' plt - The CalculatorGraph control. ' The text of the current document. ' ' RETURN VALUE: ' ' None. ' ' SIDE EFFECTS: ' ' None. ' Dim Meters As Single Dim MetersPerSecond As Single Dim Newtons As Single Dim MetersPerSecondSquared As Single Dim DeltaSeconds As Single Dim Seconds As Single Dim MaxAccel As Single Dim MaxAccelHeight As Single Dim MaxAccelVelocity As Single Dim MaxAccelTime As Single Dim NewMeters As Single Dim NewMetersPerSecond As Single Dim MaxVelocity As Single Dim MaxVelocityHeight As Single Dim MaxVelocityTime As Single Dim HypersonicVelocity As Single Dim HypersonicHeight As Single Dim HypersonicTime As Single plt.ViewWindow 0, 100000, 10000, 0, 800, 100 Meters = 30480 DeltaSeconds = 0.5 plt.Plot 3.281 * Meters, 2.237 * Abs(MetersPerSecond) While Meters > 0 Newtons = Drag(Meters, MetersPerSecond) - 53.38 MetersPerSecondSquared = Newtons / 5.443 If MetersPerSecondSquared > MaxAccel Then MaxAccel = MetersPerSecondSquared MaxAccelHeight = Meters MaxAccelVelocity = MetersPerSecond MaxAccelTime = Seconds End If NewMeters = 0.5 * MetersPerSecondSquared * DeltaSeconds ^ 2 _ + MetersPerSecond * DeltaSeconds + Meters NewMetersPerSecond = MetersPerSecond + _ MetersPerSecondSquared * DeltaSeconds Seconds = Seconds + DeltaSeconds MetersPerSecond = NewMetersPerSecond Meters = NewMeters If Abs(MetersPerSecond) > MaxVelocity Then MaxVelocity = Abs(MetersPerSecond) MaxVelocityHeight = Meters MaxVelocityTime = Seconds End If If HypersonicVelocity = 0 Then If Abs(MetersPerSecond) > VSoundHeight(Meters) Then HypersonicVelocity = Abs(MetersPerSecond) HypersonicHeight = Meters HypersonicTime = Seconds End If End If plt.Plot 3.281 * Meters, 2.237 * Abs(MetersPerSecond) plt.PlotLine DoEvents Wend WriteLn "Pumpkin drop simulation results:" WriteLn "Time (seconds): " & Format(Seconds, "#,0.0") WriteLn "Impact velocity (MPH): " & _ Format(2.237 * Abs(MetersPerSecond), "#,0.0") WriteLn "Maximum velocity: " & _ Format(2.237 * MaxVelocity, "#,0.0") WriteLn "Height of maximum velocity (ft): " & _ Format(3.281 * MaxVelocityHeight, "#,0.0") WriteLn "Time of maximum velocity: " & _ Format(MaxVelocityTime, "#,0.0") WriteLn "Maximum acceleration (MPH/sec): " & _ Format(2.2447 * MaxAccel, "#,0.0") WriteLn "Height of maximum acceleration: " & _ Format(3.281 * MaxAccelHeight, "#,0.0") WriteLn "Time of maximum acceleration: " & _ Format(MaxAccelTime, "#,0.0") WriteLn "Velocity of maximum acceleration: " & _ Format(2.237 * MaxAccelVelocity, "#,0.0") WriteLn "Hypersonic velocity: " & _ Format(2.237 * HypersonicVelocity, "#,0.0") WriteLn "Height of hypersonic velocity: " & _ Format(3.281 * HypersonicHeight, "#,0.0") WriteLn "Time of hypersonic velocity: " & _ Format(HypersonicTime, "#,0.0") WriteLn "Pumpkin drop simulation run." End Sub

Running this program gives the same results as did the calculator. The final result looks like this:

Not bad for a word processor, actually.

[ Previous page | Top of page | Next page ]

Copyright © 2002 Brian Hetrick

Page last updated 10 February 2002.

Tutorial

Building Blocks I

Control Flow II

Basic I/O

Algorithms

A First Program

Examples

Pumpkins

Simulation

Answers

Modularization

Data Structures I

Recursion

Program Attributes

Building Blocks II

Algorithm Analysis

Structuring

Data Structures II

Abstract Types

Objects

Problem Analysis

Reference Card