/** 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[]