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.

Returns

U. of Md. Computers in Physics Education HomePage

Return to part 1 of the M.U.P.P.E.T. paper by Redish and Wilson.