]> git.mxchange.org Git - flightgear.git/blob - src/WeatherCM/FGPhysicalProperties.cpp
Don't start paused for in air starts.
[flightgear.git] / src / WeatherCM / FGPhysicalProperties.cpp
1 /*****************************************************************************
2
3  Module:       FGPhysicalProperties.cpp
4  Author:       Christian Mayer
5  Date started: 28.05.99
6  Called by:    main program
7
8  -------- Copyright (C) 1999 Christian Mayer (fgfs@christianmayer.de) --------
9
10  This program is free software; you can redistribute it and/or modify it under
11  the terms of the GNU General Public License as published by the Free Software
12  Foundation; either version 2 of the License, or (at your option) any later
13  version.
14
15  This program is distributed in the hope that it will be useful, but WITHOUT
16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
18  details.
19
20  You should have received a copy of the GNU General Public License along with
21  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
22  Place - Suite 330, Boston, MA  02111-1307, USA.
23
24  Further information about the GNU General Public License can also be found on
25  the world wide web at http://www.gnu.org.
26
27 FUNCTIONAL DESCRIPTION
28 ------------------------------------------------------------------------------
29 Initialice the FGPhysicalProperties struct to something sensible(?)
30
31 HISTORY
32 ------------------------------------------------------------------------------
33 29.05.1999 Christian Mayer      Created
34 16.06.1999 Durk Talsma          Portability for Linux
35 20.06.1999 Christian Mayer      added lots of consts
36 11.10.1999 Christian Mayer      changed set<> to map<> on Bernie Bright's 
37                                 suggestion
38 19.10.1999 Christian Mayer      change to use PLIB's sg instead of Point[2/3]D
39                                 and lots of wee code cleaning
40 *****************************************************************************/
41
42 /****************************************************************************/
43 /* INCLUDES                                                                 */
44 /****************************************************************************/
45 #include "FGPhysicalProperties.h"
46 #include "FGWeatherDefs.h"
47
48 #include "FGWeatherUtils.h"
49
50 /****************************************************************************/
51 /********************************** CODE ************************************/
52 /****************************************************************************/
53 FGPhysicalProperties::FGPhysicalProperties()
54 {
55     sgVec3 zero;
56     sgZeroVec3( zero );
57     /************************************************************************/
58     /* This standart constructor fills the class with a standard weather    */
59     /************************************************************************/
60
61     Wind[-1000.0] = FGWindItem(zero);               //no Wind by default
62     Wind[10000.0] = FGWindItem(zero);               //no Wind by default
63
64     Turbulence[-1000.0] = FGTurbulenceItem(zero);   //no Turbulence by default
65     Turbulence[10000.0] = FGTurbulenceItem(zero);   //no Turbulence by default
66
67     //Initialice with the CINA atmosphere
68     Temperature[    0.0] = +15.0 + 273.16;  
69     Temperature[11000.0] = -56.5 + 273.16;                          
70     Temperature[20000.0] = -56.5 + 273.16;                          
71
72     AirPressure = FGAirPressureItem(101325.0);  
73
74     VaporPressure[-1000.0] = FG_WEATHER_DEFAULT_VAPORPRESSURE;   //in Pa (I *only* accept SI!)
75     VaporPressure[10000.0] = FG_WEATHER_DEFAULT_VAPORPRESSURE;   //in Pa (I *only* accept SI!)
76
77     //Clouds.insert(FGCloudItem())    => none
78     SnowRainIntensity = 0.0;     
79     snowRainType = Rain;            
80     LightningProbability = 0.0;
81 }
82
83 /*
84 The Methods:
85
86     WeatherPrecision getWind_x( int number ) const;
87     WeatherPrecision getWind_y( int number ) const;
88     WeatherPrecision getWind_z( int number ) const;
89     WeatherPrecision getWind_a( int number ) const;
90     void setWind_x( int number, WeatherPrecision x);
91     void setWind_y( int number, WeatherPrecision y);
92     void setWind_z( int number, WeatherPrecision z);
93     void setWind_a( int number, WeatherPrecision a);
94     WeatherPrecision getTurbulence_x( int number ) const;
95     WeatherPrecision getTurbulence_y( int number ) const;
96     WeatherPrecision getTurbulence_z( int number ) const;
97     WeatherPrecision getTurbulence_a( int number ) const;
98     void setTurbulence_x( int number, WeatherPrecision x);
99     void setTurbulence_y( int number, WeatherPrecision y);
100     void setTurbulence_z( int number, WeatherPrecision z);
101     void setTurbulence_a( int number, WeatherPrecision a);
102     WeatherPrecision getTemperature_x( int number ) const;
103     WeatherPrecision getTemperature_a( int number ) const;
104     void setTemperature_x( int number, WeatherPrecision x);
105     void setTemperature_a( int number, WeatherPrecision a);
106     WeatherPrecision getVaporPressure_x( int number ) const;
107     WeatherPrecision getVaporPressure_a( int number ) const;
108     void setVaporPressure_x( int number, WeatherPrecision x);
109     void setVaporPressure_a( int number, WeatherPrecision a);
110
111 are in the extra file FGPhysicalProperties_bind.cpp
112 */
113
114 unsigned int FGPhysicalProperties::getNumberOfCloudLayers(void) const
115 {
116     return Clouds.size();
117 }
118
119 FGCloudItem FGPhysicalProperties::getCloudLayer(unsigned int nr) const
120 {
121     map<WeatherPrecision,FGCloudItem>::const_iterator CloudsIt = Clouds.begin();
122
123     //set the iterator to the 'nr'th entry
124     for (; nr > 0; nr--)
125         CloudsIt++;
126
127     return CloudsIt->second;
128 }
129
130 ostream& operator<< ( ostream& out, const FGPhysicalProperties2D& p )
131 {
132     typedef map<FGPhysicalProperties::Altitude, FGWindItem      >::const_iterator wind_iterator;
133     typedef map<FGPhysicalProperties::Altitude, FGTurbulenceItem>::const_iterator turbulence_iterator;
134     typedef map<FGPhysicalProperties::Altitude, WeatherPrecision>::const_iterator scalar_iterator;
135
136     out << "Position: (" << p.p[0] << ", " << p.p[1] << ", " << p.p[2] << ")\n";
137     
138     out << "Stored Wind: ";
139     for (wind_iterator       WindIt = p.Wind.begin(); 
140                              WindIt != p.Wind.end(); 
141                              WindIt++)
142         out << "(" << WindIt->second.x() << ", " << WindIt->second.y() << ", " << WindIt->second.z() << ") m/s at (" << WindIt->first << ") m; ";
143     out << "\n";
144
145     out << "Stored Turbulence: ";
146     for (turbulence_iterator TurbulenceIt = p.Turbulence.begin(); 
147                              TurbulenceIt != p.Turbulence.end(); 
148                              TurbulenceIt++)
149         out << "(" << TurbulenceIt->second.x() << ", " << TurbulenceIt->second.y() << ", " << TurbulenceIt->second.z() << ") m/s at (" << TurbulenceIt->first << ") m; ";
150     out << "\n";
151
152     out << "Stored Temperature: ";
153     for (scalar_iterator     TemperatureIt = p.Temperature.begin(); 
154                              TemperatureIt != p.Temperature.end(); 
155                              TemperatureIt++)
156         out << Kelvin2Celsius(TemperatureIt->second) << " degC at " << TemperatureIt->first << "m; ";
157     out << "\n";
158
159     out << "Stored AirPressure: ";
160     out << p.AirPressure.getValue()/100.0 << " hPa at " << 0.0 << "m; ";
161     out << "\n";
162
163     out << "Stored VaporPressure: ";
164     for (scalar_iterator     VaporPressureIt = p.VaporPressure.begin(); 
165                              VaporPressureIt != p.VaporPressure.end(); 
166                              VaporPressureIt++)
167         out << VaporPressureIt->second/100.0 << " hPa at " << VaporPressureIt->first << "m; ";
168     out << "\n";
169
170     return out << "\n";
171 }
172
173
174 inline double F(const WeatherPrecision factor, const WeatherPrecision a, const WeatherPrecision b, const WeatherPrecision r, const WeatherPrecision x)
175 {
176     const double c = 1.0 / (-b + a * r);
177     return factor * c * ( 1.0 / (r + x) + a * c * log(fabs((r + x) * (b + a * x))) );
178 }
179
180 WeatherPrecision FGPhysicalProperties::AirPressureAt(const WeatherPrecision x) const
181 {
182     const double rho0 = (AirPressure.getValue()*FG_WEATHER_DEFAULT_AIRDENSITY*FG_WEATHER_DEFAULT_TEMPERATURE)/(TemperatureAt(0)*FG_WEATHER_DEFAULT_AIRPRESSURE);
183     const double G = 6.673e-11;    //Gravity; in m^3 kg^-1 s^-2
184     const double m = 5.977e24;      //mass of the earth in kg
185     const double r = 6368e3;        //radius of the earth in metres
186     const double factor = -(rho0 * TemperatureAt(0) * G * m) / AirPressure.getValue();
187
188     double a, b, FF = 0.0;
189
190     //ok, integrate from 0 to a now.
191     if (Temperature.size() < 2)
192     {   //take care of the case that there aren't enough points
193         //actually this should be impossible...
194
195         if (Temperature.size() == 0)
196         {
197             cerr << "ERROR in FGPhysicalProperties: Air pressure at " << x << " metres altiude requested,\n";
198             cerr << "      but there isn't enough data stored! No temperature is aviable!\n";
199             return FG_WEATHER_DEFAULT_AIRPRESSURE;
200         }
201
202         //ok, I've got only one point. So I'm assuming that that temperature is
203         //the same for all altitudes. 
204         a = 1;
205         b = TemperatureAt(0);
206         FF += F(factor, a, b, r, x  );
207         FF -= F(factor, a, b, r, 0.0);
208     }
209     else
210     {   //I've got at least two entries now
211         
212         //integrate 'backwards' by integrating the strip ]n,x] first, then ]n-1,n] ... to [0,n-m]
213         
214         if (x>=0.0)
215         {
216             map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
217             map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
218             
219             if (temp1->first == x)
220             {   //ignore that interval
221                 temp1--; temp2--;
222             }
223
224             bool first_pass = true;
225             while(true)
226             {
227                 if (temp2 == Temperature.end())
228                 {
229                     //temp2 doesn't exist. So cheat by assuming that the slope is the
230                     //same as between the two earlier temperatures
231                     temp1--; temp2--;
232                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
233                     b = temp1->second - a * temp1->first;
234                     temp1++; temp2++;
235                 }
236                 else
237                 {
238                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
239                     b = temp1->second - a * temp1->first;
240                 }
241                 
242                 if (first_pass)
243                 {
244                     
245                     FF += F(factor, a, b, r, x);
246                     first_pass = false;
247                 }
248                 else
249                 {
250                     FF += F(factor, a, b, r, temp2->first);
251                 }
252
253                 if (temp1->first>0.0)
254                 {
255                     FF -= F(factor, a, b, r, temp1->first);
256                     temp1--; temp2--;
257                 }
258                 else
259                 {
260                     FF -= F(factor, a, b, r, 0.0);
261                     return AirPressure.getValue() * exp(FF);
262                 }
263             }
264         }
265         else
266         {   //ok x is smaller than 0.0, so do everything in reverse
267             map<WeatherPrecision, WeatherPrecision>::const_iterator temp2 = Temperature.upper_bound(x);
268             map<WeatherPrecision, WeatherPrecision>::const_iterator temp1 = temp2; temp1--;
269             
270             bool first_pass = true;
271             while(true)
272             {
273                 if (temp2 == Temperature.begin())
274                 {
275                     //temp1 doesn't exist. So cheat by assuming that the slope is the
276                     //same as between the two earlier temperatures
277                     temp1 = Temperature.begin(); temp2++;
278                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
279                     b = temp1->second - a * temp1->first;
280                     temp2--;
281                 }
282                 else
283                 {
284                     a = (temp2->second - temp1->second)/(temp2->first - temp1->first);
285                     b = temp1->second - a * temp1->first;
286                 }
287                 
288                 if (first_pass)
289                 {
290                     
291                     FF += F(factor, a, b, r, x);
292                     first_pass = false;
293                 }
294                 else
295                 {
296                     FF += F(factor, a, b, r, temp2->first);
297                 }
298                 
299                 if (temp2->first<0.0)
300                 {
301                     FF -= F(factor, a, b, r, temp1->first);
302                     
303                     if (temp2 == Temperature.begin())
304                     {
305                         temp1 = Temperature.begin(); temp2++;
306                     }
307                     else
308                     {
309                         temp1++; temp2++;
310                     }
311                 }
312                 else
313                 {
314                     FF -= F(factor, a, b, r, 0.0);
315                     return AirPressure.getValue() * exp(FF);
316                 }
317             }
318         }
319     }
320     
321     return AirPressure.getValue() * exp(FF);
322 }