/** This program calculates the effect of the Coriolis force on a kicked
football.  (Or any other projectile).  It models the effects of the
Coriolis force in 3 dimensions.

For the equations and coordinate system used, see:
   https://en.wikipedia.org/wiki/Coriolis_effect#Rotating_sphere
*/

var t1
var t2
var e1
var e2
var n1
var n2
var u1
var u2

for useCoriolis = [false, true]
{
   // Initial velocity components
   veast  = 0 mph
   vnorth = 70 mph
   vup    = 10 mph

   // Initial positions; will terminate when up = 0 m
   east  = 0 yards
   north = 0 yards
   up    = 3 inches   // Initial height above ground.

   timestep = 0.0001 s

   omega = 1 revolution/day      // Rotation rate of the earth

   latitude = +40 degrees        // Boulder, Colorado

   t = 0 s
   while up > 0 mm
   {
      t = t + timestep
      
      // Eastward component
      if useCoriolis
         aeast = 2 omega (vnorth sin[latitude] - vup cos[latitude])
      else
         aeast = 0 m/s^2
      
      veast = veast + aeast timestep
      east = east + veast timestep

      
      // Northward component
      if useCoriolis
         anorth = 2 omega (-veast sin[latitude])
      else
         anorth = 0 m/s^2
      
      vnorth = vnorth + anorth timestep
      north = north + vnorth timestep

      
      // Upward component
      aup = -gravity                             // Constant gravity
      if useCoriolis
         aup = aup + 2 omega (veast cos[latitude])  // Add coriolis effect upwards
      vup = vup + aup timestep    // deltaV = a t
      up  = up + vup timestep     // deltaDistance = v t
   }

   println["\nuseCoriolis = $useCoriolis:\nt: "     + format[t, "s", 3]       + "\t" +
           "East: "  + format[east, "yards", 3] + "\t" +
           "North: " + format[north,"yards", 3] + "\t" +
           "Up: "    + format[up,   "yards", 3]]

   if useCoriolis == false
      [t1, e1, n1, u1] =[t, east, north, up]
   else
      [t2, e2, n2, u2] =[t, east, north, up]
}

// Calculate total deflection
dtotal = sqrt[(e2-e1)^2 + (n2-n1)^2 + (u2-u1)^2]
println["\nDeflection difference due to Coriolis effect:"]
println["dt: "     + format[t2-t1, "s", 3] + "\t" +
        "dEast: "  + format[e2-e1,"mm",3]  + "\t" +
        "dNorth: " + format[n2-n1,"mm",3]  + "\t" +
        "dUp: "    + format[u2-u1,"mm",3]  + "\n" +
        "total: " + format[dtotal,"mm",3]]