/** This is an improved, simplified model of Mars's atmosphere by Alan Eliasen
    based on NASA's Mars-GRAM observations model.  (The old version provided by
    NASA has been moved to MarsAtmosphereOld.frink )  This version has the full
    coefficients as found in the curve fit, as opposed to the simpler, truncated
    coefficients found in MarsAtmosphere.frink 

    It is intended to fit the data with a minimum number of parameters but with
    a very good fit to observed data.

    The equations are produced by parsing multiple Mars-GRAM data files by
    MarsAtmosherePlot.frink .  It is calibrated between the latitudes of
    -30 to +30 degrees.

    It is fitted on data from -5 to 170 km altitude so those should be
    considered the limits of validity, but the exponential fits on pressure and
    density seem to be quite linear over the range of the data so you might be
    justified in trying to extrapolate beyond that.
*/
class MarsAtmosphere
{
   /** Calculate the temperature at a given altitude.  This is a piecewise
       polynomial fit of degree 6 in two pieces.  Temperature is a curve that
       goes back and forth which is why this is the most complicated equation.
       It has an r-value of 0.9992498 */
   class temperature[altitude is length] :=
   {
      if altitude < 75000 m
         return 5.282633379488119848e-22 K m^-5 altitude^5 +
              -1.0843333839705542966e-16 K m^-4 altitude^4 +
               8.1990464437373983514e-12 K m^-3 altitude^3 +
              -2.4440144442633200965e-7  K m^-2 altitude^2 +
               0.00045304141784463639336 K m^-1 altitude +
               218.95470785694358513 K
      else
         return 2.46395223850889803e-22 K m^-5 altitude^5 +
             -1.5089093732887027681e-16 K m^-4 altitude^4 +
              3.5614428973920825364e-11 K m^-3 altitude^3 +
           -0.0000040171473690951854995 K m^-2 altitude^2 +
                  0.2152205818448618262 K m^-1 altitude +
                 -4225.6165589569379028 K
   }

   /** Calculate the pressure at a given altitude.  This is an exponential
       piecewise fit in two pieces.  It has an r-value of 0.99997429. */
   class pressure[altitude is length] :=
   {
      if altitude < 15000 m
         return 575.4966557798 e^(-0.000091171019659 m^-1 altitude) kg m^-1 s^-2
      else
         return 887.5439923258 e^(-0.000124149016880 m^-1 altitude) kg m^-1 s^-2
   }

   /** Calculate the atmospheric density at a given altitude.  This is an
       exponential piecewise fit in two pieces.  It has an r-value of
       0.99963887. */
   class density[altitude is length] :=
   {
      if altitude < 40000 m
         return 0.014330001122388 e^(-0.000093523163879 m^-1 altitude) kg m^-3
      else
         return 0.039535054930745 e^(-0.0001274456027426 m^-1 altitude) kg m^-3
   }

   /** Convenience method to return [temperature, pressure, density] for the
       specified altitude.  This is to match the old model. */
   class getTPD[altitude is length] :=
   {
      return [temperature[altitude], pressure[altitude], density[altitude]]
   }
}