/** This program calculates the drag through the atmosphere of the StarTram
    design:  

    https://en.wikipedia.org/wiki/StarTram
*/

use StandardAtmosphere.frink

getDensity[altitude] :=
{
   [temp, pressure] = StandardAtmosphere.getTemperatureAndPressure[altitude]
   return StandardAtmosphere.getDensity[altitude, temp, pressure]
}

getDensityVelocitySquaredProduct[altitude, velocity] :=
{
   getDensity[altitude] * velocity^2
}


altitude = 6000 m  // Top of evacuated tunnel
v = 8.8 km/s       // Fast enough for orbital velocity (with corrective burn)

mass = 40 tons
Frocket =  0 N     // StarTram is ballistic, no rocket thrust.

target = 350 km
Cd = 0.09
A = pi (1 m)^2    // 33 ft^2

launchAngle = 90 degrees   // We'll call 90 degrees straight up.

timestep = .01 s

t = 0 s
Esum = 0 J
while (altitude < target)
{
   density = getDensity[altitude]
   Fdrag = 1/2 density v^2 A Cd

   weight = G earthmass mass / (altitude + earthradius)^2   // Correct weight for decreasing gravity
   Fup = Frocket - Fdrag - weight

   a = Fup/mass                                             // Actual upward acceleration

   vnew = v + a timestep
   vave = (v + vnew) / 2

   stepdist = vave * timestep
   Esum = Esum + Fdrag * stepdist         // Energy lost this timestep

   println[(t->"s") + "\t" + format[altitude,"km", 3] + "\t" + format[v, "mach", 2] + "\t" + format[a,"gee",3] + "\t" + format[Fdrag, "lbf", 5] + "\t" + format[Esum, "MJ", 2]]
   
   v = vnew

   // TODO:  Refine this to model changing distance around earth's curve.
   altitude = altitude + stepdist sin[launchAngle]

   t = t + timestep
}