]> git.mxchange.org Git - flightgear.git/blob - src/FDM/JSBSim/models/atmosphere/FGStandardAtmosphere.cpp
Sync. With JSBSim CVS
[flightgear.git] / src / FDM / JSBSim / models / atmosphere / FGStandardAtmosphere.cpp
1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2
3  Module:       FGStandardAtmosphere.cpp
4  Author:       Jon Berndt, Tony Peden
5  Date started: 5/2011
6  Purpose:      Models the 1976 U.S. Standard Atmosphere
7  Called by:    FGFDMExec
8
9  ------------- Copyright (C) 2011  Jon S. Berndt (jon@jsbsim.org) -------------
10
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
14  version.
15
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
19  details.
20
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.
24
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.
27
28 FUNCTIONAL DESCRIPTION
29 --------------------------------------------------------------------------------
30
31
32 HISTORY
33 --------------------------------------------------------------------------------
34
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
40
41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 INCLUDES
43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
44
45 #include <iostream>
46 #include <iomanip>
47 #include <cstdlib>
48 #include "FGFDMExec.h"
49 #include "FGStandardAtmosphere.h"
50
51 namespace JSBSim {
52
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;
55
56 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 CLASS IMPLEMENTATION
58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
59
60 FGStandardAtmosphere::FGStandardAtmosphere(FGFDMExec* fdmex) : FGAtmosphere(fdmex),
61                                                                TemperatureBias(0.0),
62                                                                TemperatureDeltaGradient(0.0)
63 {
64   Name = "FGStandardAtmosphere";
65
66   StdAtmosTemperatureTable = new FGTable(9);
67
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. 
71   //
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 
84
85   LapseRateVector.resize(StdAtmosTemperatureTable->GetNumRows()-1);
86   PressureBreakpointVector.resize(StdAtmosTemperatureTable->GetNumRows());
87
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);
92
93   bind();
94   Debug(0);
95 }
96
97 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
98
99 FGStandardAtmosphere::~FGStandardAtmosphere()
100 {
101   delete StdAtmosTemperatureTable;
102   LapseRateVector.clear();
103   Debug(1);
104 }
105
106 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107
108 bool FGStandardAtmosphere::InitModel(void)
109 {
110   PressureBreakpointVector[0] = StdSLpressure = 2116.22; // psf
111   TemperatureDeltaGradient = 0.0;
112   TemperatureBias = 0.0;
113   CalculateLapseRates();
114   CalculatePressureBreakpoints();
115   Calculate(0.0);
116   StdSLtemperature = SLtemperature = Temperature;
117   SLpressure = Pressure;
118   StdSLdensity     = SLdensity = Density;
119   StdSLsoundspeed  = SLsoundspeed = Soundspeed;
120
121   rSLtemperature = 1/SLtemperature ;
122   rSLpressure    = 1/SLpressure    ;
123   rSLdensity     = 1/SLdensity     ;
124   rSLsoundspeed  = 1/SLsoundspeed  ;
125
126 //  PrintStandardAtmosphereTable();
127
128   return true;
129 }
130
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.
135
136 double FGStandardAtmosphere::GetPressure(double altitude) const
137 {
138   unsigned int b=0;
139   double pressure = 0.0;
140   double Lmb, Exp, Tmb, deltaH, factor;
141   double numRows = StdAtmosTemperatureTable->GetNumRows();
142
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)) {
150     b++;
151     testAlt = (*StdAtmosTemperatureTable)(b+1,0);
152   }
153   if (b>0) b--;
154
155   double BaseAlt = (*StdAtmosTemperatureTable)(b+1,0);
156   Tmb = GetTemperature(BaseAlt);
157   deltaH = altitude - BaseAlt;
158
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);
164   } else {
165     pressure = PressureBreakpointVector[b]*exp(-Mair*deltaH/(Rstar*Tmb));
166   }
167
168   return pressure;
169 }
170
171 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
172
173 void FGStandardAtmosphere::SetPressureSL(double pressure, ePressure unit)
174 {
175   double press = ConvertToPSF(pressure, unit);
176
177   PressureBreakpointVector[0] = press;
178   CalculatePressureBreakpoints();
179 }
180
181 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182 // Get the modeled temperature at a specified altitude, including any bias or gradient
183 // effects.
184
185 double FGStandardAtmosphere::GetTemperature(double altitude) const
186 {
187   double T = StdAtmosTemperatureTable->GetValue(altitude) + TemperatureBias;
188   if (altitude <= GradientFadeoutAltitude)
189     T += TemperatureDeltaGradient * (GradientFadeoutAltitude - altitude);
190
191   return T;
192 }
193
194 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
195 // Retrieves the standard temperature at a particular altitude.
196
197 double FGStandardAtmosphere::GetStdTemperature(double altitude) const
198 {
199   double Lk9 = 0.00658368; // deg R per foot
200   double Tinf = 1800.0; // Same as 1000 Kelvin
201   double temp = Tinf;
202
203   if (altitude < 298556.4) {                // 91 km - station 8
204
205     temp = StdAtmosTemperatureTable->GetValue(altitude);
206
207   } else if (altitude < 360892.4) {        // 110 km - station 9
208
209     temp = 473.7429 - 137.38176 * sqrt(1.0 - pow((altitude - 298556.4)/65429.462, 2.0));
210
211   } else if (altitude < 393700.8) {        // 120 km - station 10
212
213     temp = 432 + Lk9 * (altitude - 360892.4);
214
215   } else if (altitude < 3280839.9) {        // 1000 km station 12
216
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);
220
221   }
222
223   return temp;
224 }
225
226 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
227
228 double FGStandardAtmosphere::GetStdPressure(double altitude) const
229 {
230   double press=0;
231   if (TemperatureBias == 0.0 && TemperatureDeltaGradient == 0.0 && PressureBreakpointVector[0] == StdSLpressure) {
232     press = GetPressure(altitude);
233   } else if (altitude <= 100000.0) {
234     GetStdPressure100K(altitude);
235   } else {
236     // Cannot currently retrieve the standard pressure
237   }
238   return press;
239 }
240
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.
247
248 double FGStandardAtmosphere::GetStdPressure100K(double altitude) const
249 {
250   // Limit this equation to input altitudes of 100000 ft.
251   if (altitude > 100000.0) altitude = 100000.0;
252
253   double alt[5];
254   const double coef[5] = {  2116.217,
255                           -7.648932746E-2,
256                            1.0925498604E-6,
257                           -7.1135726027E-12,
258                            1.7470331356E-17 };
259
260   alt[0] = 1;
261   for (int pwr=1; pwr<=4; pwr++) alt[pwr] = alt[pwr-1]*altitude;
262
263   double press = 0.0;
264   for (int ctr=0; ctr<=4; ctr++) press += coef[ctr]*alt[ctr];
265   return press;
266 }
267
268 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
269 // Get the standard density at a specified altitude
270
271 double FGStandardAtmosphere::GetStdDensity(double altitude) const
272 {
273   return GetStdPressure(altitude)/(Reng * GetStdTemperature(altitude));
274 }
275
276 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
277
278 void FGStandardAtmosphere::SetTemperature(double t, double h, eTemperature unit)
279 {
280   double targetSLtemp = ConvertToRankine(t, unit);
281
282   TemperatureBias = 0.0;
283   TemperatureBias = targetSLtemp - GetTemperature(h);
284   CalculatePressureBreakpoints();
285 }
286
287 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
288
289 void FGStandardAtmosphere::SetTemperatureBias(eTemperature unit, double t)
290 {
291   if (unit == eCelsius || unit == eKelvin)
292     t *= 1.80; // If temp delta "t" is given in metric, scale up to English
293
294   TemperatureBias = t;
295   CalculatePressureBreakpoints();
296 }
297
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.
303
304 void FGStandardAtmosphere::SetTemperatureSL(double t, eTemperature unit)
305 {
306   SetTemperature(t, 0.0, unit);
307 }
308
309 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310 // Sets a Sea Level temperature delta that is ramped out by 86 km (282,152 ft).
311
312 void FGStandardAtmosphere::SetSLTemperatureGradedDelta(eTemperature unit, double deltemp)
313 {
314   SetTemperatureGradedDelta(deltemp, 0.0, unit);
315 }
316
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.
323
324 void FGStandardAtmosphere::SetTemperatureGradedDelta(double deltemp, double h, eTemperature unit)
325 {
326   if (unit == eCelsius || unit == eKelvin)
327     deltemp *= 1.80; // If temp delta "t" is given in metric, scale up to English
328
329   TemperatureDeltaGradient = deltemp/(GradientFadeoutAltitude - h);
330   CalculateLapseRates();
331   CalculatePressureBreakpoints();
332 }
333
334 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
335
336   void FGStandardAtmosphere::PrintStandardAtmosphereTable()
337 {
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) {
341     Calculate(i);
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
346        << std::endl;
347   }
348
349   // Re-execute the Run() method to reset the calculated values
350   Run(false);
351 }
352
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.
359
360 void FGStandardAtmosphere::CalculateLapseRates()
361 {
362   for (unsigned int bh=0; bh<LapseRateVector.size(); bh++)
363   {
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;
369   }
370 }
371
372 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
373
374 void FGStandardAtmosphere::CalculatePressureBreakpoints()
375 {
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
382                  + TemperatureBias 
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);
389     } else {
390       PressureBreakpointVector[b+1] = PressureBreakpointVector[b]*exp(-Mair*deltaH/(Rstar*Tmb));
391     }
392   }
393 }
394
395 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
396
397 void FGStandardAtmosphere::ResetSLTemperature()
398 {
399   TemperatureBias = TemperatureDeltaGradient = 0.0;
400   CalculateLapseRates();
401   CalculatePressureBreakpoints();
402 }
403
404 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
405
406 void FGStandardAtmosphere::ResetSLPressure()
407 {
408   PressureBreakpointVector[0] = StdSLpressure; // psf
409   CalculatePressureBreakpoints();
410 }
411
412 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
413
414 void FGStandardAtmosphere::bind(void)
415 {
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);
424 }
425
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
433 //       whatsoever.
434 //    1: This value explicity requests the normal JSBSim
435 //       startup messages
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
444
445 void FGStandardAtmosphere::Debug(int from)
446 {
447   if (debug_lvl <= 0) return;
448
449   if (debug_lvl & 1) { // Standard console startup message output
450     if (from == 0) { // Constructor
451     }
452   }
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;
456   }
457   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
458   }
459   if (debug_lvl & 8 ) { // Runtime state variables
460   }
461   if (debug_lvl & 16) { // Sanity checking
462   }
463   if (debug_lvl & 128) { // 
464   }
465   if (debug_lvl & 64) {
466     if (from == 0) { // Constructor
467       std::cout << IdSrc << std::endl;
468       std::cout << IdHdr << std::endl;
469     }
470   }
471 }
472
473 } // namespace JSBSim