Now that we have a function for air density as a function of height, we can (finally) compute drag.
Recall the drag equation:
Fd = ½CDAρv2
CD, the coefficient of drag, we are taking as 0.18. A, the cross-sectional area presented to the direction of travel through the air, is the area of a circle one foot in diameter. This is π(6 in)2 or 0.07297 m2. ρ, the air density, we can compute based on height. v, the velocity through the air, we simply need to be supplied. Substituting in the known quantities, we have
Fd = (0.5)(0.18)(0.07297 m2)ρv2
= 0.0065673ρv2 m2
So there are effectively two inputs to our computation of the drag force: the altitude and velocity. Once we have those, however, we can fairly easily compute the drag.
Now, testing this routine will be an interesting challenge. We do not have any known results with which to compare the output of this routine; we cannot straightforwardly graph its output, as it takes two variables.
After a little thought, however, we can identify two graphs we can make for this routine. We can make graphs of its output as a function of height while holding velocity fixed — and we can do this on the same graph for several velocities. And we can compute and graph the terminal velocity of the pumpkin as a function of height. We can perhaps judge the reasonableness of the output from these.
The first unit test, then, simply graphs the drag as a function of height, for velocities of 50, 100, 150, 200, and 250 miles per hour. The second unit test can figure out terminal velocity — the speed at which drag balances weight — as a function of height.
Since drag decreases exponentially with height, the terminal velocity also increases exponentially. Looking forward to the results of the first unit test, we can see that our second unit test, the one that computes terminal velocity, should not bother with heights greater than about 40,000 feet. (This also gives us a hint that the odds are pretty good that the pumpkin’s velocity will exceed the speed of sound in the high regions of the atmosphere — there is essentially no drag, and the pumpkin will build up considerable speed while in essentially free fall.)
We need to figure out how to get the terminal velocity. At any height, the terminal velocity will be somewhat greater than the terminal velocity at lower heights; so we can just start testing one velocity after another, going up by any convenient amount, until the drag exceeds 12 pounds, the weight of the pumpkin. This is an inefficient and inaccurate way to do it, of course, but on the other hand this is a unit test, not a production program, and the purpose of the number is to graph it, not compute with it. The resolution of the graph is at most one part in 63, or about 2%, so our accuracy does not need to be high. Were we to need to compute with the terminal velocity, we could use a zero finding method such as bisection to figure it out to any precision desired. However, we can suspect the accuracy of the approximations in use is such that perhaps the 2% resolution of the graph is about right.
We can pseudocode the terminal velocity finding portion of the unit test, then, as follows:
v := 170; for h := 0 to 40000 step 1000 do begin f := drag (h, v); while f < 12 pounds do begin v := v + 1; f := drag (h, v) end; plot (h, v) end;
(Remember we decided that the concept of a terminal velocity did not seem particularly useful above about 40,000 feet, and we estimated the terminal velocity at height 0 to be 185. We are starting the velocity search at 170.)
We can relatively easily code these unit tests, so let us do so.
Variable allocation is starting to get interesting in our calculator routines. Our philosophy so far has been to pack all variables at the start of the alphabet, and to use the first variables to pass parameters and results. As programs call other programs, they need to use the first variables to communicate with the subprograms. Any state that programs need to save across calls to other programs therefore needs to be saved in variables just after the start of the alphabet. The density function DENSH, for example, takes a single parameter (the height) and returns a single result (the density). In the absence of a need to save internal state, it would require only one variable, A, for this communication. However, DENSH also consumes variables B and C. The first variable available for the drag function’s use is actually D. It is for considerations such as that that program documentation — in particular, a list of the resources used by the program — is so critical. (The little trick of using a list to save the caller’s state is looking better and better, isn’t it?)
At any rate, here are the programs.
Program DRAG (43 bytes):
; Program DRAG ; ; Compute the drag force in Newtons acting on a medium pumpkin as a ; function of the pumpkin's height in meters and velocity in meters ; per second ; ; Variables: ; A is height [on entry, input to DENSH], air density [output from ; DENSH], drag [on exit] ; B is velocity [on entry], destroyed by DENSH ; C is destroyed by DENSH ; D is saved velocity ; ; Symbols: ; -> is assignment arrow ; x^2 is x squared operator B->D ; Save velocity in D Prog "DENSH" ; Get density as a function of height in A .0065673ADx^2->A ; Compute drag into A
Program UTDRAG1 (133 bytes):
; Program UTDRAG1 ; ; Graph the drag in pounds against the height in feet for velocities ; of 50 to 250 miles per hour ; ; Variables: ; A is height in meters [input to DRAG], drag in newtons [output from ; DRAG] ; B is velocity in meters/second [input to DRAG] ; C-D are destroyed by DRAG ; E is velocity in miles/hour ; F is height in feet ; ; Symbols: ; -> is assignment arrow ; <> is not equal to relational ViewWindow 0,100000,10000,0,25,5 ; Set the graph view window For 50->E To 250 Step 50 ; Loop through velocities For 0->F To 100000 Step 1000 ; Loop through heights .3048F->A ; Get height in meters .447E->B ; Get velocity in meters/second Prog "DRAG" ; Compute the drag .2248A->A ; Convert drag to pounds Plot F,A ; Plot the drag against height If F<>0 ; Test if this is the first point Then Line ; If not, connect the dots IfEnd ; [End of test if this is the first point] Next ; [End of loop through heights] Next ; [End of loop through velocities]
Running this unit test gives us the following graph:
This graph indicates that, for a given velocity, drag decreases exponentially with height — which is exactly what we would expect, given that air density decreases exponentially with height. The vertical axis is pounds, with tick marks every five pounds; the graph also indicates that at height 0, a velocity of 150 miles per hour yields a drag of about 8 pounds, and one of 200 miles per hour yields a drag of about 14 pounds. We would expect the terminal velocity of the (12 pound) pumpkin at height 0 to be between these two values, probably about 185 miles per hour.
Program UTDRAG2 (149 bytes):
; Program UTDRAG2 ; ; Find and plot the terminal velocity of the 12-pound pumpkin as a ; function of height ; ; Variables: ; A is altitude [input to DRAG], drag [output from DRAG] ; B is velocity [input to DRAG] ; C-D are destroyed by DRAG ; E is current velocity estimate in miles/hour ; F is current altitude in feet ; ; Symbols: ; -> is assignment arrow ; >= is greater than or equal to relational ; <> is not equal to relational ViewWindow 0,40000,10000,0,400,100 ; Set the view window 170->E ; Start terminal velocity at 170 (we know For 0->F To 40000 Step 1000 ; Loop through all altitudes of interest While 1 ; Start a loop to search for velocity .3048F->A ; Get altitude in meters .447E->B ; Get current velocity in meters/second Prog "DRAG" ; Compute drag .2248A->A ; Convert drag into pounds If A>=12 ; If the drag is at least 12 pounds Then Break ; then exit the search loop IfEnd ; [End of test of drag] E+1->E ; Try 1 mile/hour more than the current ; guess at the velocity WhileEnd ; [End of loop to search for velocity] Plot F,E ; Plot the point If F<>0 ; If this is not the first point Then Line ; then connect the dots IfEnd ; [End of test if the first point] Next ; [End of altitude loop]
When we run this program, we get the following interesting graph of terminal velocity of the pumpkin as a function of height:
The X axis is height, from 0 to 40,000 feet, and the Y axis is terminal velocity, from 0 to 400 miles per hour. The graph indicates the terminal velocity of the pumpkin at the surface is about 185 miles per hour. We vaguely remember the terminal velocity of a skydiver as about 150 or 200 miles per hour, so the terminal velocity of a pumpkin being about 185 seems reasonable.
We can code the drag function and the two unit tests as follows.
Private Function Drag(ByVal Meters As Single, _ ByVal MetersPerSecond As Single) As Single ' FUNCTIONAL DESCRIPTION: ' ' Computes the drag in Newtons on the medium pumpkin as a ' function of the pumpkin's height in meters and velocity in ' meters per second. ' ' FORMAL PARAMETERS: ' ' Height.rf.v - The pumpkin's height in meters. ' Speed.rf.v - The pumpkin's velocity in meters/second. ' ' IMPLICIT INPUTS: ' ' None. ' ' IMPLICIT OUTPUTS: ' ' None. ' ' RETURN VALUE: ' ' The drag on the pumpkin in Newtons. ' ' SIDE EFFECTS: ' ' None. ' Dim sngDensity As Single sngDensity = Density(Meters) Drag = 0.0065673 * sngDensity * MetersPerSecond ^ 2 End Function
Private Sub UTDrag1() ' FUNCTIONAL DESCRIPTION: ' ' Unit tests the Drag routine. ' ' Computes and graphs the drag of a function of height for ' velocities of 50, 100, 150, 200, and 250 miles per hour. ' ' 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 Feet As Single Dim Pounds As Single Dim MPH As Single plt.ViewWindow 0, 100000, 10000, 0, 25, 5 For MPH = 50 To 250 Step 50 For Feet = 0 To 100000 Step 1000 Pounds = 0.2248 * Drag(0.3048 * Feet, 0.447 * MPH) plt.Plot Feet, Pounds If Feet <> 0 Then plt.PlotLine End If DoEvents Next Feet Next MPH WriteLn "Drag unit test 1 run." End Sub
Private Sub UTDrag2() ' FUNCTIONAL DESCRIPTION: ' ' Unit tests the Drag routine. ' ' Computes and graphs the terminal velocity of a pumpkin at ' heights from 0 to 100,000 feet in increments of 1,000 feet. ' ' 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 Feet As Single Dim MPH As Single Dim Pounds As Single plt.ViewWindow 0, 40000, 10000, 0, 400, 100 MPH = 170 For Feet = 0 To 40000 Step 1000 Do Pounds = 0.2248 * Drag(0.3048 * Feet, 0.447 * MPH) If Pounds >= 12 Then Exit Do End If MPH = MPH + 1 Loop plt.Plot Feet, MPH If Feet <> 0 Then plt.PlotLine End If DoEvents Next Feet WriteLn "Drag unit test 2 run." End Sub
Running these two unit tests in Microsoft Word gives the same results as the calculator unit tests.
Copyright © 2002 Brian Hetrick
Page last updated 10 February 2002.
Building Blocks I
Control Flow II
A First Program
Data Structures I
Building Blocks II
Data Structures II