Wrote a README file that clearly describes the simulation. Also updated the simulation code to be clear and consistent with the documentation. Fixes #36

trick_models/parachute/include/Parachutist.hh

@@ -14,24 +14,14 @@ LIBRARY DEPENDENCIES:
 class Parachutist {
 
     public:
-        double position ;                    /* m     xyz-position */
+        double altitude ;                    /* m     xyz-position */
         double velocity ;                    /* m/s   xyz-velocity */
         double acceleration ;                /* m/s2  xyz-acceleration  */
 
-        double crossSectionalArea;           /* m2 */
+        double area;                         /* m2 */
-        double crossSectionalAreaRate;       /* m2/s */
+        double Cd;                           /* -- */
-        double cooefficientOfDrag;           /* -- */
-        double cooefficientOfDragRate;       /* -- */
         double mass;                         /* kg */
 
-        double parachuteDeploymentStartTime; /* s */
-        double parachuteDeploymentDuration;  /* s */
-
-        static const double crossSectionalArea_freefall;  /* m2 */
-        static const double crossSectionalArea_parachute; /* m2 */
-        static const double cooefficientOfDrag_freefall;  /* -- */
-        static const double cooefficientOfDrag_parachute; /* -- */
-
         bool touchDown;                      /* -- */
 
         REGULA_FALSI rf ;
@@ -40,10 +30,7 @@ class Parachutist {
         int state_init();
         int state_deriv();
         int state_integ();
-        int parachute_control();
+        double touch_down(double groundAltitude);
-        double touch_down();
-    private:
-        double parachuteDeploymentEndTime;   /* s */
 };
 
 #endif

trick_models/parachute/src/Parachutist.cpp

@@ -13,7 +13,6 @@ PROGRAMMERS:
 // http://www.engineeringtoolbox.com/standard-atmosphere-d_604.html
 // U.S Standard Atmosphere Air Properties in SI Units
 #define NUM_ELEMENTS 21
-
 // Units = meters (above sea level).
 const double altitude_array[NUM_ELEMENTS] = {
  -1000.0,     0.0,  1000.0,  2000.0,  3000.0,  4000.0,  5000.0,  6000.0,
@@ -33,78 +32,43 @@ const double gravity_array[NUM_ELEMENTS] = {
 9.684, 9.654, 9.624, 9.594, 9.564
 };
 
-const double Parachutist::crossSectionalArea_freefall  =  0.75 ;
-const double Parachutist::crossSectionalArea_parachute = 30.00;
-const double Parachutist::cooefficientOfDrag_freefall  =  0.75;
-const double Parachutist::cooefficientOfDrag_parachute =  1.30;
-
 int Parachutist::default_data() {
 
-    position = 38969.3; //38969.3 meters = 127852 feet
+    altitude = 38969.3; //38969.3 meters = 127852 feet
     velocity = 0.0;
-    acceleration = 0.0;
+    area = 0.75;
-    parachuteDeploymentStartTime = 259; /* 4 minutes, 19 seconds*/
+    Cd = 0.75;
-    parachuteDeploymentDuration  = 3;
-    crossSectionalArea = crossSectionalArea_freefall;
-    cooefficientOfDrag = cooefficientOfDrag_freefall;
     touchDown = false;
     mass = 82.0;
-
     return (0);
 }
 
 int Parachutist::state_init() {
-    parachuteDeploymentEndTime = parachuteDeploymentStartTime + parachuteDeploymentDuration;
-    return (0);
-}
-
-#include "sim_services/Executive/include/exec_proto.h"
-int Parachutist::parachute_control() {
-
-    double currentTime =  exec_get_sim_time();
-
-    if ((currentTime >  parachuteDeploymentStartTime) &&
-        (currentTime <= parachuteDeploymentEndTime)) {
-        cooefficientOfDragRate = (cooefficientOfDrag_parachute - cooefficientOfDrag_freefall)
-                / parachuteDeploymentDuration;
-        crossSectionalAreaRate = (crossSectionalArea_parachute - crossSectionalArea_freefall)
-                / parachuteDeploymentDuration;
-    } else {
-        cooefficientOfDragRate = 0.0;
-        crossSectionalAreaRate = 0.0;
-    }
     return (0);
 }
 
 int Parachutist::state_deriv() {
 
-   // Calculate the force of gravity.
+// Calculate the forces and acceleration.
-#if 1
-   double g = 9.81;
-#else
-   double g = interpolate( position, altitude_array, gravity_array, NUM_ELEMENTS );
-#endif
 
+   // Calculate Force of gravity.
+   double g = interpolate( altitude, altitude_array, gravity_array, NUM_ELEMENTS );
    double Force_gravity = mass * -g;
 
-   // Calculate the force of drag.
+   // Calculate Force of drag.
-   double air_density = interpolate( position, altitude_array, air_density_array, NUM_ELEMENTS );
+   double air_density = interpolate( altitude, altitude_array, air_density_array, NUM_ELEMENTS );
-   double Force_drag = cooefficientOfDrag * 0.5 * air_density * velocity * velocity * crossSectionalArea;
+   double Force_drag = Cd * 0.5 * air_density * velocity * velocity * area;
 
-   // Sum the forces and calculate acceleration.
+   // Calculate Total Force.
    double Force_total;
-
+   if ( !touchDown ) {
-   // If skydiver has touched down then set state derivatives to zero.
+       Force_total = Force_gravity + Force_drag ;
-   if ( touchDown ) {
+       acceleration = Force_total / mass;
+   } else {
        Force_total = 0.0;
        velocity = 0.0;
        acceleration = 0.0;
-   } else {
-       Force_total = Force_gravity + Force_drag ;
-       acceleration = Force_total / mass;
    }
-
-   /* RETURN: -- Always return zero. */
    return(0);
 }
 
@@ -112,39 +76,20 @@ int Parachutist::state_deriv() {
 int Parachutist::state_integ() {
 
    int integration_step;
-
+   load_state( &altitude, &velocity, (double*)0);
-   load_state( &position,
+   load_deriv ( &velocity, &acceleration, (double*)0);
-               &velocity,
-               &cooefficientOfDrag,
-               &crossSectionalArea,
-               (double*)0
-             );
-
-   load_deriv ( &velocity,
-                &acceleration,
-                &cooefficientOfDragRate,
-                &crossSectionalAreaRate,
-                (double*)0
-              );
-
    integration_step = integrate();
-
+   unload_state( &altitude, &velocity, (double*)0);
-   unload_state( &position,
-                 &velocity,
-                 &cooefficientOfDrag,
-                 &crossSectionalArea,
-                 (double*)0
-               );
 
    return(integration_step);
 }
 
-double Parachutist::touch_down() {
+double Parachutist::touch_down(double groundAltitude) {
     double tgo ;
     double now ;
 
     /* error function: how far the skydiver is above ground */
-    rf.error = position ;
+    rf.error = altitude - groundAltitude ;
 
     now = get_integ_time() ;
     tgo = regula_falsi( now, &rf ) ;

trick_sims/SIM_parachute/DP_Product/DP_parachutist.xml

@@ -10,7 +10,7 @@
             <title>Plot</title>
             <curve>
                 <var units="s">sys.exec.out.time</var>
-                <var units="m">dyn.parachutist.position</var>
+                <var units="m">dyn.parachutist.altitude</var>
             </curve>
         </plot>
     </page>

trick_sims/SIM_parachute/Modified_data/parachutist.dr

@@ -11,10 +11,10 @@ drg.append(trick.DRBinary("parachutist"))
 drg[DR_GROUP_ID].set_freq(trick.DR_Always)
 drg[DR_GROUP_ID].set_cycle(0.1)
 drg[DR_GROUP_ID].set_single_prec_only(False)
-drg[DR_GROUP_ID].add_variable("dyn.parachutist.position")
+drg[DR_GROUP_ID].add_variable("dyn.parachutist.altitude")
 drg[DR_GROUP_ID].add_variable("dyn.parachutist.velocity")
 drg[DR_GROUP_ID].add_variable("dyn.parachutist.acceleration")
-drg[DR_GROUP_ID].add_variable("dyn.parachutist.crossSectionalArea")
+drg[DR_GROUP_ID].add_variable("dyn.parachutist.area")
-drg[DR_GROUP_ID].add_variable("dyn.parachutist.cooefficientOfDrag")
+drg[DR_GROUP_ID].add_variable("dyn.parachutist.Cd")
 trick.add_data_record_group(drg[DR_GROUP_ID], trick.DR_Buffer)
 drg[DR_GROUP_ID].enable()

trick_sims/SIM_parachute/README.md

@@ -0,0 +1,106 @@
+### Background
+On October 14, 2012, 43-year-old Austrian daredevel Felix Baumgartner broke the
+world record for the highest-ever skydive after jumping from a balloon at an
+altitude of 127,852 feet. He reached a top speed of 843.6 mph, spent approximately
+4 minutes 19 seconds in free-fall, and landed safely approximately 11 minutes
+after jumping.
+
+### Simulation
+The simulation only considers the forces of gravity and drag, and only motion in
+the vertical. The acceleration of the skydiver is determined by summing the
+forces of gravity and drag acting on him and then dividing by his mass. His
+velocity is determined by integrating his acceleration over time, and his
+altitude by integrating his velocity over time.
+
+Desired outputs are:
+* Plot of altitude vs. time.
+* Plot of velocity vs. time.
+* Time of touchdown.
+
+##### Gravity
+<math xmlns="http://www.w3.org/1998/Math/MathML">
+<mrow>
+    <msub>
+        <mi>F</mi>
+        <mi>g</mi>
+    </msub>
+    <mo>=</mo>
+    <mi>m</mi>
+    <mo>&InvisibleTimes;</mo>
+    <mi>g</mi>
+</mrow>
+</math>
+
+Where:
+
+* m = mass of the skydiver.
+* g = acceleration of gravity.
+
+##### Drag
+<math xmlns="http://www.w3.org/1998/Math/MathML">
+    <mrow>
+        <msub>
+            <mi>F</mi>
+            <mi>d</mi>
+        </msub>
+        <mo>=</mo>
+        <mfrac>
+            <mn>1</mn>
+            <mn>2</mn>
+        </mfrac>
+        <mo>&InvisibleTimes;</mo>
+        <mi>&#0961;</mi>
+        <mo>&InvisibleTimes;</mo>
+        <msup>
+            <mi>v</mi>
+            <mn>2</mn>
+        </msup>
+        <mo>&InvisibleTimes;</mo>
+        <msub>
+            <mi>C</mi>
+            <mi>d</mi>
+        </msub>
+        <mo>&InvisibleTimes;</mo>
+        <mi>A</mi>
+    </mrow>
+</math>
+
+Where:
+
+* C<sub>d</sub> = Coefficient of drag
+* &#0961; = air density
+* v = instantaneous velocity
+* A = cross-sectional area of the jumper
+
+#### Air Density and Gravity Data
+The table at:
+<http://www.engineeringtoolbox.com/standard-atmosphere-d_604.html>
+provides both gravity (g) and air density (&#0961;) at various altitudes.
+From these data we interpolate, to approximate the air density and gravity at
+specific altitudes.
+
+#### Parachute Deployment
+Parachute deployment is modeled, using a Trick event (in input.py) that simply increasing skydiver's
+
+1) cross-sectional area and
+2) coefficient of drag.
+
+at the specified time.
+
+### Jump Scenario
+dyn.groundAltitude = 1000
+dyn.parachutist.altitude = 38969.6 meters
+dyn.parachutist.velocity = 0.0
+dyn.parachutist.area = 0.75
+dyn.parachutist.Cd = 0.75
+dyn.parachutist.mass = 82.0
+
+# At 4 minutes and 19 seconds, pop the chute.
+dyn.parachutist.Cd = 1.3
+dyn.parachutist.area = 30.0
+
+#### Results
+
+![Plot of Altitude vs Time](images/plot_altitude_vs_time.png "Altitude vs. Time")
+![Plot of Velocity vs Time](images/plot_velocity_vs_time.png "Velocity vs. Time")
+

trick_sims/SIM_parachute/RUN_Baumgartner_jump/input.py

@@ -1,7 +1,27 @@
-#execfile("Modified_data/realtime.py")
+
 execfile("Modified_data/parachutist.dr")
-
 trick.TMM_reduced_checkpoint(0)
-dyn_integloop.getIntegrator(trick.Runge_Kutta_4, 4)
+dyn_integloop.getIntegrator(trick.Runge_Kutta_4, 2)
 
+# The ground is a 1000 meters above sea level.
+dyn.groundAltitude = 1000
+
+# State of the skydiver when he jumps
+dyn.parachutist.altitude = 38969.6 # 127852 feet
+dyn.parachutist.velocity = 0.0
+dyn.parachutist.area = 0.75
+dyn.parachutist.Cd = 0.75
+dyn.parachutist.mass = 82.0
+
+# At 4 minutes and 19 seconds, pop the chute.
+deployChute = trick.new_event("deployChute")
+deployChute.condition(0, "trick.exec_get_sim_time() >= 259.0")
+deployChute.action(0, "dyn.parachutist.Cd = 1.3");
+deployChute.action(1, "dyn.parachutist.area = 30.0");
+deployChute.action(2, "print \"Parachute Deployed.\"");
+deployChute.set_cycle(1.0)
+deployChute.activate()
+trick.add_event(deployChute)
+
+# Run for 800 seconds.
 trick.stop(800)

trick_sims/SIM_parachute/S_define

@@ -13,6 +13,7 @@ class ParachutistSimObject : public Trick::SimObject {
 
     public:
     Parachutist parachutist;
+    double groundAltitude;
 
     ParachutistSimObject() {
 
@@ -20,8 +21,7 @@ class ParachutistSimObject : public Trick::SimObject {
         ("initialization")  parachutist.state_init() ;
         ("derivative")      parachutist.state_deriv() ;
         ("integration") trick_ret = parachutist.state_integ() ;
-        ("dynamic_event")   parachutist.touch_down() ;
+        ("dynamic_event")   parachutist.touch_down(groundAltitude) ;
-        (0.02, "scheduled") parachutist.parachute_control();
     }
 } ;