/** This plots atmospheric data for Mars based on the Mars-GRAM data. It is used to create the much better MarsAtmosphere.frink model 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) It currently is modeled between the latitudes between -30 deg and 30 deg where most landings occur. For a map of landing spots, see: https://www.planetary.org/space-images/mars_landing_site_map_lakdawalla */ use Grid.frink use MarsAtmosphereOld.frink use curveFit.frink use functionUtils.frink use LeastSquares.frink url = filenameToURL["/home/eliasen/builds/Mars-GRAM 2010/Release1.0_Nov10/txtFiles/tpdloy11.txt"] maxlat = 30 deg //minheight = -5 km //maxheight = 200 km minheight = -5 km maxheight = 170 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] = MarsAtmosphereOld.getTPD[alt] ptemp1.addPoint[t/K, -alt/km] } // Fit the values to a piecewise polynomial of the specified degree [tf1, rtf1] = LeastSquares.fitDegreePiecewise[meantemp, 5] println[inputForm[tf1]] println["rValue: $rtf1"] [tempFunc, rtemp] = quadraticFitFunctionPiecewise3[meantemp] //tempFunc = quadraticFitFunction[meantemp] //rtemp = rValue[tempFunc, toArray[meantemp]] //println[inputForm[functionBody[tempFunc]]] println[inputForm[tempFunc]] println["rvalue = $rtemp"] ptempfit = new polyline ptempfit1 = new polyline for alt = sort[meantemp.keys[]] { ptempfit.addPoint[tempFunc[alt]/K, -alt/km] ptempfit1.addPoint[tf1[alt]/K, -alt/km] } gtemp = new graphics gtemp.font["Monospaced", .05] gtemp.add[ptemp] gtemp.color[0,0,1] gtemp.add[ptemp1] gtemp.color[1,0,0] gtemp.add[ptempfit] gtemp.color[0,1,0] gtemp.add[ptempfit1] 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[] gtemp.write["MarsAtmosphereTemperature.svg", 1000, 800] gtemp.write["MarsAtmosphereTemperature.png", 1000, 800] 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] = MarsAtmosphereOld.getTPD[alt] ppressure1.addPoint[log[p / Pa] dummyx, -alt/km] } [pressureFunc, rpressure] = exponentialFitFunctionPiecewise[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.caption["log10[pressure/pascal]"] gpressure.show[] gpressure.write["MarsAtmospherePressure.svg", 1000, 800] gpressure.write["MarsAtmospherePressure.png", 1000, 800] 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] = MarsAtmosphereOld.getTPD[alt] if d != 0 kg/m^3 pdensity1.addPoint[log[d / (kg/m^3)] dummyx, -alt/km] } [densityFunc, rdensity] = exponentialFitFunctionPiecewise[meandensity] 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.caption["log10[density/(kg/m^3)]"] gdensity.show[] gdensity.write["MarsAtmosphereDensity.svg", 1000, 800] gdensity.write["MarsAtmosphereDensity.png", 1000, 800]