1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 Module: FGStandardAtmosphere.cpp
4 Author: Jon Berndt, Tony Peden
6 Purpose: Models the 1976 U.S. Standard Atmosphere
9 ------------- Copyright (C) 2011 Jon S. Berndt (jon@jsbsim.org) -------------
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU Lesser General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
21 You should have received a copy of the GNU Lesser General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 59 Temple
23 Place - Suite 330, Boston, MA 02111-1307, USA.
25 Further information about the GNU Lesser General Public License can also be found on
26 the world wide web at http://www.gnu.org.
28 FUNCTIONAL DESCRIPTION
29 --------------------------------------------------------------------------------
33 --------------------------------------------------------------------------------
35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 COMMENTS, REFERENCES, and NOTES
37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38 [1] Anderson, John D. "Introduction to Flight, Third Edition", McGraw-Hill,
39 1989, ISBN 0-07-001641-0
41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
48 #include "FGFDMExec.h"
49 #include "FGStandardAtmosphere.h"
53 static const char *IdSrc = "$Id: FGStandardAtmosphere.cpp,v 1.20 2011/09/18 12:06:21 bcoconni Exp $";
54 static const char *IdHdr = ID_STANDARDATMOSPHERE;
56 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
60 FGStandardAtmosphere::FGStandardAtmosphere(FGFDMExec* fdmex) : FGAtmosphere(fdmex),
62 TemperatureDeltaGradient(0.0)
64 Name = "FGStandardAtmosphere";
66 StdAtmosTemperatureTable = new FGTable(9);
68 // This is the U.S. Standard Atmosphere table for temperature in degrees
69 // Rankine, based on geometric altitude. The table values are often given
70 // in literature relative to geopotential altitude.
72 // GeoMet Alt Temp GeoPot Alt GeoMet Alt
73 // (ft) (deg R) (km) (km)
74 // -------- -------- ---------- ----------
75 *StdAtmosTemperatureTable << 0.0 << 518.67 // 0.000 0.000
76 << 36151.6 << 390.0 // 11.000 11.019
77 << 65823.5 << 390.0 // 20.000 20.063
78 << 105518.4 << 411.6 // 32.000 32.162
79 << 155347.8 << 487.2 // 47.000 47.350
80 << 168677.8 << 487.2 // 51.000 51.413
81 << 235570.9 << 386.4 // 71.000 71.802
82 << 282152.2 << 336.5 // 84.852 86.000
83 << 298556.4 << 336.5; // 91.000 - First layer in high altitude regime
85 LapseRateVector.resize(StdAtmosTemperatureTable->GetNumRows()-1);
86 PressureBreakpointVector.resize(StdAtmosTemperatureTable->GetNumRows());
88 // Assume the altitude to fade out the gradient at is at the highest
89 // altitude in the table. Above that, other functions are used to
90 // calculate temperature.
91 GradientFadeoutAltitude = (*StdAtmosTemperatureTable)(StdAtmosTemperatureTable->GetNumRows(),0);
97 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99 FGStandardAtmosphere::~FGStandardAtmosphere()
101 delete StdAtmosTemperatureTable;
102 LapseRateVector.clear();
106 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108 bool FGStandardAtmosphere::InitModel(void)
110 PressureBreakpointVector[0] = StdSLpressure = 2116.22; // psf
111 TemperatureDeltaGradient = 0.0;
112 TemperatureBias = 0.0;
113 CalculateLapseRates();
114 CalculatePressureBreakpoints();
116 StdSLtemperature = SLtemperature = Temperature;
117 SLpressure = Pressure;
118 StdSLdensity = SLdensity = Density;
119 StdSLsoundspeed = SLsoundspeed = Soundspeed;
121 rSLtemperature = 1/SLtemperature ;
122 rSLpressure = 1/SLpressure ;
123 rSLdensity = 1/SLdensity ;
124 rSLsoundspeed = 1/SLsoundspeed ;
126 // PrintStandardAtmosphereTable();
131 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132 // Get the actual pressure as modeled at a specified altitude
133 // These calculations are from equations 33a and 33b in the U.S. Standard Atmosphere
134 // document referenced in the documentation for this code.
136 double FGStandardAtmosphere::GetPressure(double altitude) const
139 double pressure = 0.0;
140 double Lmb, Exp, Tmb, deltaH, factor;
141 double numRows = StdAtmosTemperatureTable->GetNumRows();
143 // Iterate through the altitudes to find the current Base Altitude
144 // in the table. That is, if the current altitude (the argument passed in)
145 // is 20000 ft, then the base altitude from the table is 0.0. If the
146 // passed-in altitude is 40000 ft, the base altitude is 36151.6 ft (and
147 // the index "b" is 2 - the second entry in the table).
148 double testAlt = (*StdAtmosTemperatureTable)(b+1,0);
149 while ((altitude >= testAlt) && (b <= numRows-2)) {
151 testAlt = (*StdAtmosTemperatureTable)(b+1,0);
155 double BaseAlt = (*StdAtmosTemperatureTable)(b+1,0);
156 Tmb = GetTemperature(BaseAlt);
157 deltaH = altitude - BaseAlt;
159 if (LapseRateVector[b] != 0.00) {
160 Lmb = LapseRateVector[b];
161 Exp = Mair/(Rstar*Lmb);
162 factor = Tmb/(Tmb + Lmb*deltaH);
163 pressure = PressureBreakpointVector[b]*pow(factor, Exp);
165 pressure = PressureBreakpointVector[b]*exp(-Mair*deltaH/(Rstar*Tmb));
171 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
173 void FGStandardAtmosphere::SetPressureSL(double pressure, ePressure unit)
175 double press = ConvertToPSF(pressure, unit);
177 PressureBreakpointVector[0] = press;
178 CalculatePressureBreakpoints();
181 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182 // Get the modeled temperature at a specified altitude, including any bias or gradient
185 double FGStandardAtmosphere::GetTemperature(double altitude) const
187 double T = StdAtmosTemperatureTable->GetValue(altitude) + TemperatureBias;
188 if (altitude <= GradientFadeoutAltitude)
189 T += TemperatureDeltaGradient * (GradientFadeoutAltitude - altitude);
194 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
195 // Retrieves the standard temperature at a particular altitude.
197 double FGStandardAtmosphere::GetStdTemperature(double altitude) const
199 double Lk9 = 0.00658368; // deg R per foot
200 double Tinf = 1800.0; // Same as 1000 Kelvin
203 if (altitude < 298556.4) { // 91 km - station 8
205 temp = StdAtmosTemperatureTable->GetValue(altitude);
207 } else if (altitude < 360892.4) { // 110 km - station 9
209 temp = 473.7429 - 137.38176 * sqrt(1.0 - pow((altitude - 298556.4)/65429.462, 2.0));
211 } else if (altitude < 393700.8) { // 120 km - station 10
213 temp = 432 + Lk9 * (altitude - 360892.4);
215 } else if (altitude < 3280839.9) { // 1000 km station 12
217 double lambda = 0.00001870364;
218 double eps = (altitude - 393700.8) * (20855531.5 + 393700.8) / (20855531.5 + altitude);
219 temp = Tinf - (Tinf - 648.0) * exp(-lambda*eps);
226 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
228 double FGStandardAtmosphere::GetStdPressure(double altitude) const
231 if (TemperatureBias == 0.0 && TemperatureDeltaGradient == 0.0 && PressureBreakpointVector[0] == StdSLpressure) {
232 press = GetPressure(altitude);
233 } else if (altitude <= 100000.0) {
234 GetStdPressure100K(altitude);
236 // Cannot currently retrieve the standard pressure
241 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
242 // This function calculates an approximation of the standard atmospheric pressure
243 // up to an altitude of about 100,000 ft. If the temperature and pressure are not
244 // altered for local conditions, the GetPressure(h) function should be used,
245 // as that is valid to a much higher altitude. This function is accurate to within
246 // a couple of psf up to 100K ft. This polynomial fit was determined using Excel.
248 double FGStandardAtmosphere::GetStdPressure100K(double altitude) const
250 // Limit this equation to input altitudes of 100000 ft.
251 if (altitude > 100000.0) altitude = 100000.0;
254 const double coef[5] = { 2116.217,
261 for (int pwr=1; pwr<=4; pwr++) alt[pwr] = alt[pwr-1]*altitude;
264 for (int ctr=0; ctr<=4; ctr++) press += coef[ctr]*alt[ctr];
268 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
269 // Get the standard density at a specified altitude
271 double FGStandardAtmosphere::GetStdDensity(double altitude) const
273 return GetStdPressure(altitude)/(Reng * GetStdTemperature(altitude));
276 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
278 void FGStandardAtmosphere::SetTemperature(double t, double h, eTemperature unit)
280 double targetSLtemp = ConvertToRankine(t, unit);
282 TemperatureBias = 0.0;
283 TemperatureBias = targetSLtemp - GetTemperature(h);
284 CalculatePressureBreakpoints();
287 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
289 void FGStandardAtmosphere::SetTemperatureBias(eTemperature unit, double t)
291 if (unit == eCelsius || unit == eKelvin)
292 t *= 1.80; // If temp delta "t" is given in metric, scale up to English
295 CalculatePressureBreakpoints();
298 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
299 // This function calculates a bias based on the supplied temperature for sea
300 // level. The bias is applied to the entire temperature profile at all altitudes.
301 // Internally, the Rankine scale is used for calculations, so any temperature
302 // supplied must be converted to that unit.
304 void FGStandardAtmosphere::SetTemperatureSL(double t, eTemperature unit)
306 SetTemperature(t, 0.0, unit);
309 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310 // Sets a Sea Level temperature delta that is ramped out by 86 km (282,152 ft).
312 void FGStandardAtmosphere::SetSLTemperatureGradedDelta(eTemperature unit, double deltemp)
314 SetTemperatureGradedDelta(deltemp, 0.0, unit);
317 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
318 // Sets a temperature delta at the supplied altitude that is ramped out by 86 km.
319 // After this calculation is performed, the lapse rates and pressure breakpoints
320 // must be recalculated. Since we are calculating a delta here and not an actual
321 // temperature, we only need to be concerned about a scale factor and not
322 // the actual temperature itself.
324 void FGStandardAtmosphere::SetTemperatureGradedDelta(double deltemp, double h, eTemperature unit)
326 if (unit == eCelsius || unit == eKelvin)
327 deltemp *= 1.80; // If temp delta "t" is given in metric, scale up to English
329 TemperatureDeltaGradient = deltemp/(GradientFadeoutAltitude - h);
330 CalculateLapseRates();
331 CalculatePressureBreakpoints();
334 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
336 void FGStandardAtmosphere::PrintStandardAtmosphereTable()
338 std::cout << "Altitude (ft) Temp (F) Pressure (psf) Density (sl/ft3)" << std::endl;
339 std::cout << "------------- -------- -------------- ----------------" << std::endl;
340 for (int i=0; i<280000; i+=1000) {
342 std::cout << std::setw(12) << std::setprecision(2) << i
343 << " " << std::setw(9) << std::setprecision(2) << Temperature - 459.67
344 << " " << std::setw(13) << std::setprecision(4) << Pressure
345 << " " << std::setw(18) << std::setprecision(8) << Density
349 // Re-execute the Run() method to reset the calculated values
353 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354 // This function calculates (or recalculates) the lapse rate over an altitude range
355 // where the "bh" in this case refers to the index of the base height in the
356 // StdAtmosTemperatureTable table. This function should be called anytime the
357 // temperature table is altered, such as when a gradient is applied across the
358 // temperature table for a range of altitudes.
360 void FGStandardAtmosphere::CalculateLapseRates()
362 for (unsigned int bh=0; bh<LapseRateVector.size(); bh++)
364 double t0 = (*StdAtmosTemperatureTable)(bh+1,1);
365 double t1 = (*StdAtmosTemperatureTable)(bh+2,1);
366 double h0 = (*StdAtmosTemperatureTable)(bh+1,0);
367 double h1 = (*StdAtmosTemperatureTable)(bh+2,0);
368 LapseRateVector[bh] = (t1 - t0) / (h1 - h0) + TemperatureDeltaGradient;
372 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
374 void FGStandardAtmosphere::CalculatePressureBreakpoints()
376 for (unsigned int b=0; b<PressureBreakpointVector.size()-1; b++) {
377 double BaseTemp = (*StdAtmosTemperatureTable)(b+1,1);
378 double BaseAlt = (*StdAtmosTemperatureTable)(b+1,0);
379 double UpperAlt = (*StdAtmosTemperatureTable)(b+2,0);
380 double deltaH = UpperAlt - BaseAlt;
381 double Tmb = BaseTemp
383 + (GradientFadeoutAltitude - BaseAlt)*TemperatureDeltaGradient;
384 if (LapseRateVector[b] != 0.00) {
385 double Lmb = LapseRateVector[b];
386 double Exp = Mair/(Rstar*Lmb);
387 double factor = Tmb/(Tmb + Lmb*deltaH);
388 PressureBreakpointVector[b+1] = PressureBreakpointVector[b]*pow(factor, Exp);
390 PressureBreakpointVector[b+1] = PressureBreakpointVector[b]*exp(-Mair*deltaH/(Rstar*Tmb));
395 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
397 void FGStandardAtmosphere::ResetSLTemperature()
399 TemperatureBias = TemperatureDeltaGradient = 0.0;
400 CalculateLapseRates();
401 CalculatePressureBreakpoints();
404 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
406 void FGStandardAtmosphere::ResetSLPressure()
408 PressureBreakpointVector[0] = StdSLpressure; // psf
409 CalculatePressureBreakpoints();
412 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
414 void FGStandardAtmosphere::bind(void)
416 typedef double (FGStandardAtmosphere::*PMFi)(int) const;
417 typedef void (FGStandardAtmosphere::*PMF)(int, double);
418 PropertyManager->Tie("atmosphere/delta-T", this, eRankine,
419 (PMFi)&FGStandardAtmosphere::GetTemperatureBias,
420 (PMF)&FGStandardAtmosphere::SetTemperatureBias);
421 PropertyManager->Tie("atmosphere/SL-graded-delta-T", this, eRankine,
422 (PMFi)&FGStandardAtmosphere::GetTemperatureDeltaGradient,
423 (PMF)&FGStandardAtmosphere::SetSLTemperatureGradedDelta);
426 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
427 // The bitmasked value choices are as follows:
428 // unset: In this case (the default) JSBSim would only print
429 // out the normally expected messages, essentially echoing
430 // the config files as they are read. If the environment
431 // variable is not set, debug_lvl is set to 1 internally
432 // 0: This requests JSBSim not to output any messages
434 // 1: This value explicity requests the normal JSBSim
436 // 2: This value asks for a message to be printed out when
437 // a class is instantiated
438 // 4: When this value is set, a message is displayed when a
439 // FGModel object executes its Run() method
440 // 8: When this value is set, various runtime state variables
441 // are printed out periodically
442 // 16: When set various parameters are sanity checked and
443 // a message is printed out when they go out of bounds
445 void FGStandardAtmosphere::Debug(int from)
447 if (debug_lvl <= 0) return;
449 if (debug_lvl & 1) { // Standard console startup message output
450 if (from == 0) { // Constructor
453 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
454 if (from == 0) std::cout << "Instantiated: FGStandardAtmosphere" << std::endl;
455 if (from == 1) std::cout << "Destroyed: FGStandardAtmosphere" << std::endl;
457 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
459 if (debug_lvl & 8 ) { // Runtime state variables
461 if (debug_lvl & 16) { // Sanity checking
463 if (debug_lvl & 128) { //
465 if (debug_lvl & 64) {
466 if (from == 0) { // Constructor
467 std::cout << IdSrc << std::endl;
468 std::cout << IdHdr << std::endl;
473 } // namespace JSBSim