MarsAtmospherePlot.frink

Download or view MarsAtmospherePlot.frink in plain text format


/** This plots atmospheric data for Mars based on the Mars-GRAM data.

    See the Mars-GRAM User's Guide at:
https://ntrs.nasa.gov/api/citations/20210023957/downloads/Mars-GRAM%20User%20Guide.pdf
    (jump to page 62 of the PDF file which is labeled page 49 on the document itself)

   For a map of landing spots, see:
   https://www.planetary.org/space-images/mars_landing_site_map_lakdawalla
*/

use Grid.frink
use MarsAtmosphere.frink
use curveFit.frink
use functionUtils.frink

url = filenameToURL["/home/eliasen/builds/Mars-GRAM 2010/Release1.0_Nov10/txtFiles/tpdloy11.txt"]

maxlat = 30 deg

g = new graphics

dtemp     = new dict
ddensity  = new dict
dpressure = new dict

dummyx =!= dummyx   

for line = rest[lines[url]]  // Skips first line of headers
{
   line = trim[line]
   fields = eval[split[%r/\s+/, line]]
   long =   fields@0 degrees
   height = fields@1 km
   lat =    fields@2 deg
   temp =   fields@3 K
   density = fields@8 kg/m^3
   pressure = fields@13 Pa

   // Change these to focus on a certain latitude.
   // See
   // https://www.planetary.org/space-images/mars_landing_site_map_lakdawalla
   if abs[lat] <= maxlat
   {
      dtemp.addToList[height, temp]
      ddensity.addToList[height, density]
      dpressure.addToList[height, pressure]
   }
}

// Medium altitude file
url = filenameToURL["/home/eliasen/builds/Mars-GRAM 2010/Release1.0_Nov10/txtFiles/tpdms101.txt"]

for line = rest[lines[url]]  // Skips first line of headers
{
   line = trim[line]
   fields = eval[split[%r/\s+/, line]]
   long =   fields@0 degrees
   height = fields@1 km
   lat =    fields@2 deg
   temp =   fields@3 K
   pressure = fields@8 Pa
   density = fields@13 kg/m^3

   // Change these to focus on a certain latitude.
   // See
   // https://www.planetary.org/space-images/mars_landing_site_map_lakdawalla
   if abs[lat] <= maxlat
   {
      dtemp.addToList[height, temp]
      ddensity.addToList[height, density]
      dpressure.addToList[height, pressure]
   }
}

println["\nTemperature"]
meantemp = new dict
ptemp = new polyline
ptemp1 = new polyline
for alt = sort[dtemp.keys[]]
{
   meantemp@alt = meanAndSD[dtemp@alt, false]@0
   println["$alt\t" + meantemp@alt]
   ptemp.addPoint[meantemp@alt/K, -alt/km]
   [t,p,d] = MarsAtmosphere.getTPD[alt]
   ptemp1.addPoint[t/K, -alt/km]
}

gtemp = new graphics
gtemp.font["Monospaced", .05]
gtemp.add[ptemp]
gtemp.color[0,0,1]
gtemp.add[ptemp1]
grid = new Grid
grid.setUnits[1, -1]
grid.auto[gtemp]
gtemp.add[grid.getGrid[]]
gtemp.caption["Temperature, K", "top"]
gtemp.caption["altitude (km)", "right", -90 deg]
gtemp.show[]

println["\nPressure"]
meanpressure = new dict
ppressure = new polyline
ppressure1 = new polyline
for alt = sort[dpressure.keys[]]
{
   meanpressure@alt = meanAndSD[dpressure@alt, false]@0
   println["$alt\t" + meanpressure@alt]
   ppressure.addPoint[log[meanpressure@alt / Pa] dummyx, -alt/km]
   [t,p,d] = MarsAtmosphere.getTPD[alt]
   ppressure1.addPoint[log[p / Pa] dummyx, -alt/km]
}

pressureFunc = exponentialFitFunction1[meanpressure]
println[inputForm[functionBody[pressureFunc]]]
ppressurefit = new polyline
for alt = sort[meanpressure.keys[]]
   ppressurefit.addPoint[log[pressureFunc[alt] / Pa] dummyx, -alt/km]

gpressure = new graphics
gpressure.add[ppressure]
gpressure.color[0,0,1]
gpressure.add[ppressure1]
gpressure.color[1,0,0]
gpressure.add[ppressurefit]
grid = new Grid
grid.setUnits[dummyx, -1]
grid.auto[gpressure]
gpressure.add[grid.getGrid[]]
gpressure.caption["Pressure", "top"]
gpressure.show[]


println["\nDensity"]
meandensity = new dict
pdensity = new polyline
pdensity1 = new polyline
for alt = sort[ddensity.keys[]]
{
   meandensity@alt = meanAndSD[ddensity@alt, false]@0
   println["$alt\t" + meandensity@alt]
   pdensity.addPoint[log[meandensity@alt / (kg/m^3)] dummyx, -alt/km]
   [t,p,d] = MarsAtmosphere.getTPD[alt]
   if d != 0 kg/m^3
      pdensity1.addPoint[log[d / (kg/m^3)] dummyx, -alt/km]
}

densityFunc = exponentialFitFunction1[meandensity]
println[inputForm[functionBody[densityFunc]]]
pdensityfit = new polyline
for alt = sort[meandensity.keys[]]
   pdensityfit.addPoint[log[densityFunc[alt] / (kg/m^3)] dummyx, -alt/km]


gdensity = new graphics
gdensity.add[pdensity]
gdensity.color[0,0,1]
gdensity.add[pdensity1]
gdensity.color[1,0,0]
gdensity.add[pdensityfit]
grid = new Grid
grid.setUnits[dummyx, -1]
grid.auto[gdensity]
gdensity.add[grid.getGrid[]]
gdensity.caption["Density", "top"]
gdensity.show[]


/** Calculates the mean and standard deviation of an array or
    enumerating expression.

    This uses Welford's algorithm as cited in Knuth, The Art of Computer
    Programming, Vol. 2, 3rd edition, page 232.

    Arguments:
     [list, sample]

    where
     list: is an array or enumerating expression of the items to average.
     sample: a boolean flag indicating if we want the standard deviation to be
         the sample standard variation (=true)  or the population standard
         variation (=false).

    Returns:
     [mean, sd, number]
   */

meanAndSD[list, sample] :=
{
   M = 0 list@0     // Make units of measure come out right
   S = 0 (list@0)^2 // Make units of measure come out right
   k = 1

   for v = list
   {
      oldM = M
      diff = v-oldM
      M = M + diff / k
      S = S + diff * (v - M)
      k = k+1
   }

   if sample == true
      sub = 1            // Sample standard deviation, subtract 1 from num
   else
      sub = 0            // Population standard deviation
   
   return [M, sqrt[S/(k-1-sub)], k-1]
}


Download or view MarsAtmospherePlot.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 19979 days, 18 hours, 24 minutes ago.