Calculator Programming Tutorial

Programming Building Blocks I

Examples

Ballistic Pumpkins

Simulation

Program Analysis

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;
while h > 0
do
    begin
    f := drag (h, v) – 53.38;
    a := f / m;
    hnew := .5 * a * Δt2 + v * Δt + h;
    vnew := a * Δt + v;
    t := t + Δt;
    v := vnew;
    h := hnew
    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;
vmax := 0;
hvmax := 0;
tvmax := 0;
while h > 0
do
    begin
    f := drag (h, v) – 53.38;
    a := f / 5.443;
    hnew := .5 * a * Δt2 + v * Δt + h;
    vnew := a * Δt + v;
    t := t + Δt;
    v := vnew;
    h := hnew;
    if |v| > vmax
    then
        begin
        vmax := |v|;
        hvmax := h;
        tvmax := t
        end
    end;
{ report t, v, vmax, hvmax, tvmax }

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;
vmax := 0;
hvmax := 0;
tvmax := 0;
amax := 0;
hamax := 0;
tamax := 0;
vamax := 0;
while h > 0
do
     begin
    f := drag (h, v) – 53.38;
    a := f / 5.443;
    if a > amax
    then
        begin
        amax := a;
        hamax := h;
        tamax := t;
        vamax := v
        end;
    hnew := .5 * a * Δt2 + v * Δt + h;
    vnew := a * Δt + v;
    t := t + Δt;
    v := vnew;
    h := hnew;
    if |v| > vmax
    then
        begin
        vmax := v;
        hvmax := h;
        tvmax := t
        end
    end;
{ report t, v, vmax, hvmax, tvmax, amax, hamax, tamax, vamax }

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;
vmax := 0;
hvmax := 0;
tvmax := 0;
amax := 0;
hamax := 0;
tamax := 0;
vamax := 0;
hhyp := 0;
thyp := 0;
vhyp := 0;
while h > 0
do
     begin
    f := drag (h, v) – 53.38;
    a := f / 5.443;
    if a > amax
    then
        begin
        amax := a;
        hamax := h;
        tamax := t;
        vamax := v
        end;
    hnew := .5 * a * Δt2 + v * Δt + h;
    vnew := a * Δt + v;
    t := t + Δt;
    v := vnew;
    h := hnew;
    if |v| > vmax
    then
        begin
        vmax := v;
        hvmax := h;
        tvmax := t
        end;
    if vhyp = 0
    then
        begin
        vsnd := vsound (h);
        if |v| > vsnd
        then
            begin
            vhyp := v;
            thyp := t;
            hhyp := h
            end
        end
    end;
{ report t, v, vmax, hvmax, tvmax, amax, hamax, tamax, vamax,
thyp, vhyp, hhyp }

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.

Calculator Code

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:

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:

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.

Visual Basic for Application Code

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.