/** This is a simplified model to examine the burning of an Estes rocket.
    It is a program written in the programming language Frink,
    https://frinklang.org/

    This was inspired by a tweet from Dianna Cowern, 
    https://twitter.com/thephysicsgirl/status/1302000329224007682

    It asks,
    "Anyone know why those little model rocket engines have an initial spike
     in thrust, then level out?" and it has a graph of force curves, showing
     an initial spike in thrust.

    This model examines the implications of an initially higher thrust.  Is it
    beneficial?

    Executive Summary:
    Yes, initially higher thrust is beneficial.  Burning harder earlier gives
    the rocket more energy than would be obtained by a constant burn, or by
    burning harder later.

    Some sample results for different scenarios:
*/

// Initial conditions.  Note you can mix units of measure freely
// in Frink and it will do the right thing with them.
v = 0 mph
t = 0 s
h = 0 ft
timestep = .001 s

// Mass of the rocket.  This is a guess.  Plug in your own value.
// Higher masses will show more clearly why early, stronger burns are benefical.
m = 150 g

do
{
   ft = f[t]   // Rocket force at a given instant, supplied by the f[t] function.
   F = ft - m gravity  // Net force is rocket force up minus gravity force down.
   a = F / m           // Calculate acceleration
   v = v + a timestep  // Recalculate velocity from deltaV = a timestep
   h = h + v timestep  // Recalculate distance from deltaD = v timestep
   E = 1/2 m v^2 + m gravity h  // Calculate kinetic energy plus potential energy
   t = t + timestep    // Increment timestep
} while ft > 0 N       // While force function is greater than 0 N


// Print final time, velocity, height, and net energy at burnout time.
println["t="+format[t, "s", 3] + "\tv=" + format[v, "m/s", 3] + "\th=" + format[h, "m", 3] + "\tE=" + format[E, "J", 3]]


// Force function.  This is modifiable to produce any force curve
// (maybe taken from a rocket engine specification.)
//
// It currently is just set up to produce two different force
// strengths:
//
// From time t=0 s to firstStageDuration, it returns firstStageForce.
// from time t=firstStageDuration to firstStageDuration+secondStageDuration
//     (this phase has a duration of secondStageDuration), it returns
//     secondStageForce.
// From time firstStageDuration+secondStageDuration to infinity, it returns
// 0 newtons (this is after burnout.)
//
// To see the effects that different burn profiles have, alter
// the numbers below (but preserve the net sum of duration*force if you
// want to compare one burn profile to another.)
//
// For example, the numbers below very roughly model the
// thrust curve posted by Diana Cowern at:
// https://twitter.com/thephysicsgirl/status/1302000329224007682
//
// About 10 N for .2 seconds, followed by about 5 N for 1.7 seconds.
//
// You can reverse these to have the "harder" burn at the end of the
// model, or a constant burn (both stages with same force), or
// an even more severe burn at the beginning.
f[t] :=
{
   firstStageDuration = .2 s
   firstStageForce = 10 N
   
   secondStageDuration = 1.7 s
   secondStageForce = 5 N

   if t < firstStageDuration
      return firstStageForce

   if t < firstStageDuration + secondStageDuration
      return secondStageForce

   return 0 N
}