/** This is a 2-dimensional ballistics / aerospace simulation that accurately
    models air resistance (using the U.S. Standard Atmosphere) and the curvature
    of the earth (in 2 dimensions.)

    This program models the entry of the Meteor Crater, Arizona impactor
    through layers of atmosphere and demonstrates the huge energy radiation
    due to pushing through layers of atmosphere (463 kilotons TNT).

    This program can be altered to model many types of impactors.  Some may
    slow down to terminal velocity.
*/

use StandardAtmosphere.frink

/*for v = mach 0 to mach 5 step mach 1/10
   println[(v->mach) + "\t" + Cd[v]]
exit[]*/

// h is height above ground (initial height)
h = 900 km

radius =  50 m / 2
volume = 4/3 pi radius^3
density = 7.874 g/cm^3

// m1 is the mass of the projectile
m1 = volume density

// Cd is drag coefficient of projectile
Cd = 1.4

// area is area of projectile in direction of travel.
// area of a sphere can be calculated from mass and density as:
// area = (3/4)^(2/3) m1^(2/3) pi^(1/3) density^(-2/3)

density = 7.874 g/cm^3
area = (3/4)^(2/3) m1^(2/3) pi^(1/3) density^(-2/3)
//println["Area is " + (area->"cm^2")]
//area = 0.7 m^2

// Initial velocity in the x direction (horizontal)
vx = 9.19 km/s

// Initial velocity in the y direction (vertical, positive is up)
vy = -9.19 km/s

// mp is mass of planet
mp = earthmass
r = earthradius

// x and y are a cartesian coordinate system with the center of the planet
// at x=0 m, y=0 m.  The projectile typically begins its journey at x=0 and
// at a given height-above-ground.
x = -700 km
y = r + h

initialGeopotentialHeight = (r * h) / (r + h)
//println["Geopotential height = " + (geopotentialHeight -> "ft")]
initialGeopotentialEnergy = m1 gravity initialGeopotentialHeight

initialKineticEnergy = 1/2 m1 (vx^2 + vy^2)

timestep = .01 s

t = 0 s

// Energy lost to drag
Edrag = 0 J

g = new graphics
line = new polyline
height = new polyline

do
{
   // l is distance from center of earth
   l2 = x^2 + y^2
   l = sqrt[l2]
   h = l - r

   // Angle with respect to center of the earth
   alpha = arctan[x,y]

   // Force due to gravity
   fg = - G m1 mp / l2

   // Acceleration due to gravity
   ag =  fg / m1
   agx = ag sin[alpha]
   agy = ag cos[alpha]

   // Atmospheric drag
   v2 = vx^2 + vy^2
   v = sqrt[v2]

   // Angle of travel (0 is in x direction, 90 degrees in y direction)
   theta = arctan[vy, vx]

   [temp, pressure] = StandardAtmosphere.getTemperatureAndPressure[h]
   density = StandardAtmosphere.getDensity[h, temp, pressure]

   // Calculate drag coefficient as a function of velocity
   Cd = Cd[v]

   fdrag = 1/2 density v2 Cd area

   adrag =  -fdrag / m1
   adragx = adrag cos[theta]
   adragy = adrag sin[theta]

   t = t + timestep

   // Total acceleration
   axtotal = agx + adragx
   aytotal = agy + adragy
   atotal = sqrt[axtotal^2 + aytotal^2]
   
   dvx = axtotal timestep
   dvy = aytotal timestep
   
   vx = vx + dvx
   vy = vy + dvy

   dx = vx timestep
   dy = vy timestep

   // Power due to drag
   Pdrag = fdrag v
   
   // Energy lost to drag
   // E = f * d = f * v * t
   Edrag = Edrag + Pdrag timestep
   
   x = x + dx
   y = y + dy

   line.addPoint[x/m, -y/m]
   height.addPoint[t, -h/m]

   geopotentialHeight = (r * h) / (r + h)
   geopotentialEnergy = m1 gravity geopotentialHeight

   kineticEnergy = 1/2 m1 (vx^2 + vy^2)

   totalEnergy = geopotentialEnergy + kineticEnergy
//if t mod (1 s) == 0 s
      println[formatFixed[t,"s",2] + "\t" + format[x,"m",8] + "\t" + format[y,"m",9] + "\t" + format[h,"m",9] + "\t" + format[adragx,"gee",4] + "\t" + format[adragy,"gee",4] + "\t" + format[adrag,"gee",4] + "\t" + format[v,"m/s",7] + "\t" + formatEng[Pdrag, "W", 3]]
} while h >= 5600 ft

initialEnergy = initialGeopotentialEnergy + initialKineticEnergy
println["Initial potential energy = $initialGeopotentialEnergy"]
println["Initial kinetic energy   = $initialKineticEnergy"]
println["Final kinetic energy     = " + (1/2 m1 v2)]
println["Initial energy           = " + initialEnergy]
println["Energy lost to drag      = $Edrag\t(" + formatEng[Edrag, "tons TNT", 3] + ")"]
println["Fraction of energy lost to drag = " + format[Edrag / initialEnergy, "percent", 4]]

g.add[line]
g.show[]

g2=new graphics
g2.add[height]
g2.show[]

// Drag coefficient equation for a sphere in the subsonic-to-hypersonic environment
// See ESTIMATING THE DRAG COEFFICIENTS OF METEORITES FOR ALL MACH NUMBER REGIMES.
// R. T. Carter, P. S. Jandir, and M. E. Kress,
// 40th Lunar and Planetary Science Conference (2009)
// They don't actually display the "quadratic" behavior which occurs below
// about Mach 0.8, so this is a curve-fit.
// For hypersonic velocities, this is about 0.92, while other authors state that
// a constant 0.7 is okay.
Cd[v] :=
{
   M0 = v/mach
   
   if (M0 >= 0.8)
   {
      // Mach 0.8 and above, use eq. from paper.
      M1 = M0 + 0.35
      return 2.1 e^(-1.2 M1) - 8.9 e^(-2.2 M1) + 0.92
   } else
   {
      // Below Mach 0.8
      // Use "quadratic" range picked off chart and fitted
      return 0.424 + 0.472 M0^2
   }
}