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
minheight = -5 km
maxheight = 200 km

g = new graphics

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

dummyx =!= dummyx

AITemp[h is height] :=
{
   h1 = h / km
   return K (242 - 0.087 h1 + 1.15 × 10^-4 h^2)
}

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 and height <= maxheight and height >= minheight
   {
      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 and height <= maxheight and height >= minheight
   {
      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 = exponentialFitFunction[meanpressure]
rpressure = rValue[pressureFunc, toArray[meanpressure]]
//println[inputForm[functionBody[pressureFunc]]]
println[inputForm[pressureFunc]]
println["rvalue = $rpressure"]
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 = exponentialFitFunction[meandensity]
rdensity = rValue[densityFunc, toArray[meandensity]]
//println[inputForm[functionBody[densityFunc]]]
println[inputForm[densityFunc]]
println["rvalue = $rdensity"]
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[]


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 19988 days, 15 hours, 31 minutes ago.