Proj1D: Code for a sample M.U.P.P.E.T. Program
PROGRAM Projectile1D; { proj1d.pas }
{***********************************************}
{* *}
{* Program to calculate motion of *}
{* a particle in 1D with gravity *}
{* and air resistance using RK2. *}
{* *}
{***********************************************}
USES
Crt,Dos,Graph,Printer,MUPPET;
CONST
numData : Integer = 200; {Number of points to plot}
g : Real = 9.8; {m/sec/sec}
VAR
t,x : DataVector; { time, position }
v,a : DataVector; { velocity, accel }
x0,v0 : Real; { initial conds. }
m : Real; { mass }
b : Real; { air resis. coeff. }
dt : Real; { time step }
i : Integer; { loop variable }
IC : Screen; { data screen }
act : Char; { control character }
{ The types 'DataVector" and "Screen" are defined }
{ inside the unit MUPPET. }
{*------------ Physics Procedures -------------*}
FUNCTION Force(x,v,t:Real) : Real;
BEGIN
Force := -m*g - b*v*abs(v)
END;
{*---------- Mathematics Procedures -----------*}
{ Second order Runge-Kutta routine for stepping }
{ from variables at time t (In variables) to }
{ variables at time t+dt (Out variables). }
PROCEDURE StepRK2(xIn, vIn, tIn, aIn,tStep:Real;
VAR xOut,vOut,tOut,aOut:Real);
VAR
xHalf,vHalf : Real;
tHalf,aHalf : Real;
BEGIN
tHalf := tIn + 0.5*tStep;
xHalf := xIn + 0.5*vIn*tStep;
vHalf := vIn + 0.5*aIn*tStep;
aHalf := Force(xHalf,vHalf,tHalf)/m;
tOut := tIn + tStep;
xOut := xIn + vHalf*tStep;
vOut := vIn + aHalf*tStep;
aOut := Force(xOut,vOut,tOut)/m;
END;
{*--------- Data Screen Procedures ------------*}
PROCEDURE MakeDataScreen;
BEGIN
DefineInputport(0,0.45,0,0.9);
_A[01]:='"M.U.P.P.E.T." ';
_A[02]:='"University of Maryland" ';
_A[03]:=' ';
_A[04]:='"PROJECTILE PROGRAM: 1D" ';
_A[05]:='"F = -mg - bv*abs(v)" ';
_A[06]:=' ';
_A[07]:='"PARAMETERS" ';
_A[08]:=' "Mass m = " 0.14++ "kg" ';
_A[09]:=' ';
_A[10]:=' "Air Resistance" ';
_A[11]:=' "Coefficient, b = " 0+++++ "kg/m"';
_A[12]:=' ';
_A[13]:=' "Time step, dt = " 0.050+ "sec" ';
_A[14]:=' ';
_A[15]:='"INITIAL CONDITIONS" ';
_A[16]:=' "Position: x0 = " 0++++ "m" ';
_A[17]:=' "Velocity: v0 = " 30+++ "m/sec"';
LoadScreen(IC,17);
END;
PROCEDURE GetScreenData(VAR m,b,x0,v0,dt:Real);
BEGIN
ClearMUPPETport;
Message('Press to plot, to quit');
Accept(IC); {displays screen}
m := Valu(IC,1); {puts 1st entry on IC into m}
b := Valu(IC,2); {puts 2nd entry on IC into b}
dt := Valu(IC,3);
x0 := Valu(IC,4); {etc...}
v0 := Valu(IC,5);
END;
{*----------- Graphics Procedures -------------*}
PROCEDURE GraphSetUp;
BEGIN
GraphBackColor:=DarkGray;
DefineViewport(1, 0.55,1, 0.5,0.9); {Define ViewPort 1}
DefineViewport(2, 0.55,1, 0.05,0.45); {ViewPort 2}
DefineScale(1, 0, 10, -75.0, 75); {Define Scale 1}
DefineScale(2, 0, 10, -75.0, 75); {Define Scale 2}
END;
PROCEDURE PlotIt(viewPort,color:Integer; x,y:DataVector;
nameLabel:BigStr);
BEGIN
Setcolor(color);
SelectScale(viewPort);
OpenViewport(viewPort);
Axis(0,0,1,20);
PlotData(x,y,numData);
PutLabel(Inside,nameLabel);
END;
{*--------------- Main Program ----------------*}
BEGIN
MUPPETinit;
MakeDataScreen;
GraphSetUp;
REPEAT
GetScreenData(m,b,x0,v0,dt);
IF EscapedFrom(IC) THEN
BEGIN
MUPPETdone;
EXIT
END;
t[1] := 0; {initializes first step}
x[1] := x0;
v[1] := v0;
a[1] := -g - b*v0*abs(v0)/m;
FOR i := 2 to numData DO {solve the equation}
StepRK2(x[i-1],v[i-1],t[i-1],a[i-1],dt,
x[i], v[i], t[i], a[i]);
Message('Press for new data, to quit');
PlotIt(1, lightGreen, t, x, 'X vs T');
PlotIt(2, lightRed, t, v, 'V vs T');
act := ReadKey;
UNTIL ord(act) = 27;
MUPPETdone;
END.
Some Links:
This page prepared 24. March 1995 by
Edward F. Redish
Department of Physics
University of Maryland
College Park, MD 20742
Phone: (301) 405-6120
Email: redish@quark.umd.edu