// This program draws the effects of a rocket changing an orbit.
//  All of the burns use the same change in speed  but burn in
// 4 different directions.
//

/** This file demonstrates how thrust in a rocket changes an orbit.   

   White:          Typical orbit
     Red:  Stop:   Slows down rocket
   Green:  Go:     Speeds up rocket
    Blue:  Down:   Burns downward toward planet
  Yellow:  Up:     Burns upward away from planet

   Each velocity change is regarded as instantaneous and occurs at the "top"
  of the diagram.

*/

[orbitalHeight, dv] = eval[input["Enter orbital changes",[["Height above ground", "1000 km"], ["Velocity change", "1600 km/hr"]]]]

y = earthradius + orbitalHeight
x = 0 miles

vy = 0 m/s

// Calculate horizontal velocity for a circular orbit.
// This is derived sybolically by orbit.frink.  It's amazing. Check it out.
vx = G^(1/2) earthmass^(1/2) (earthradius + orbitalHeight)^(-1/2)

println["vx = " + (vx -> "mph")]

// Calculate circular orbital period.
// This is derived sybolically by orbit.frink.  It's amazing. Check it out.
period = 2 pi vx^-1 (earthradius + orbitalHeight)
println["Period is " + (period -> HMS)]

g = new graphics
g.backgroundColor[0,0,0]
g.color[0,0,.5]
g.fillEllipseCenter[0, 0, 2 earthradius/mile, 2 earthradius/mile]
g.color[1,1,1]
g.stroke[30 km/mile]

drawOrbit[g, x, y, vx, vy, period]     // White: Nominal orbit

g.color[0,.5,1]                  // Blue: down: Accelerate towards earth.
drawOrbit[g, x, y, vx, vy-dv, period]

g.color[1,1,0]                  // Yellow: up:  Accelerate up
drawOrbit[g, x, y, vx, vy+dv, period]

g.color[1,0,0]
drawOrbit[g, x, y, vx-dv, vy, period]  // Red:  stop: Backwards (slow down)

g.color[0,1,0]
drawOrbit[g, x, y, vx+dv, vy, period*2]  // Green:  go: Forwards (speed up)

g.show[]
g.write["orbitalBurn.svg", 1000, undef]
g.write["orbitalBurn.html", 1000, undef]
g.write["orbitalBurn.png", 1000, undef]


drawOrbit[g, x, y, vx, vy, period=220 min] :=
{
   p = new polyline
   
   timestep = period/1000
   for t = 0 min to period step timestep
   {
      dsqr = x^2 + y^2
      a = G earthmass / dsqr
      theta = arctan[y,x]
      ax = - a cos[theta]
      ay = - a sin[theta]

      vx = vx + timestep ax
      vy = vy + timestep ay
      x = x + timestep vx
      y = y + timestep vy
      p.addPoint[x/mile,-y/mile]
   }

   g.add[p]
}