]> git.mxchange.org Git - flightgear.git/blob - src/FDM/JSBSim/models/atmosphere/FGWinds.cpp
Merge branch 'next' into comm-subsystem
[flightgear.git] / src / FDM / JSBSim / models / atmosphere / FGWinds.cpp
1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2
3  Module:       FGWinds.cpp
4  Author:       Jon Berndt, Tony Peden, Andreas Gaeb
5  Date started: Extracted from FGAtmosphere, which originated in 1998
6                5/2011
7  Purpose:      Models winds, gusts, turbulence, and other atmospheric disturbances
8  Called by:    FGFDMExec
9
10  ------------- Copyright (C) 2011  Jon S. Berndt (jon@jsbsim.org) -------------
11
12  This program is free software; you can redistribute it and/or modify it under
13  the terms of the GNU Lesser General Public License as published by the Free Software
14  Foundation; either version 2 of the License, or (at your option) any later
15  version.
16
17  This program is distributed in the hope that it will be useful, but WITHOUT
18  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
20  details.
21
22  You should have received a copy of the GNU Lesser General Public License along with
23  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
24  Place - Suite 330, Boston, MA  02111-1307, USA.
25
26  Further information about the GNU Lesser General Public License can also be found on
27  the world wide web at http://www.gnu.org.
28
29 FUNCTIONAL DESCRIPTION
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 <cstdlib>
47 #include "FGWinds.h"
48 #include "FGFDMExec.h"
49
50 using namespace std;
51
52 namespace JSBSim {
53
54 static const char *IdSrc = "$Id: FGWinds.cpp,v 1.5 2011/09/11 11:36:04 bcoconni Exp $";
55 static const char *IdHdr = ID_WINDS;
56
57 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 CLASS IMPLEMENTATION
59 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
60
61 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62 // square a value, but preserve the original sign
63
64 static inline double square_signed (double value)
65 {
66     if (value < 0)
67         return value * value * -1;
68     else
69         return value * value;
70 }
71
72 /// simply square a value
73 static inline double sqr(double x) { return x*x; }
74
75 FGWinds::FGWinds(FGFDMExec* fdmex) : FGModel(fdmex)
76 {
77   Name = "FGWinds";
78
79   MagnitudedAccelDt = MagnitudeAccel = Magnitude = 0.0;
80   SetTurbType( ttMilspec );
81   TurbGain = 1.0;
82   TurbRate = 10.0;
83   Rhythmicity = 0.1;
84   spike = target_time = strength = 0.0;
85   wind_from_clockwise = 0.0;
86   psiw = 0.0;
87
88   vGustNED.InitMatrix();
89   vTurbulenceNED.InitMatrix();
90
91   // Milspec turbulence model
92   windspeed_at_20ft = 0.;
93   probability_of_exceedence_index = 0;
94   POE_Table = new FGTable(7,12);
95   // this is Figure 7 from p. 49 of MIL-F-8785C
96   // rows: probability of exceedance curve index, cols: altitude in ft
97   *POE_Table
98            << 500.0 << 1750.0 << 3750.0 << 7500.0 << 15000.0 << 25000.0 << 35000.0 << 45000.0 << 55000.0 << 65000.0 << 75000.0 << 80000.0
99     << 1   <<   3.2 <<    2.2 <<    1.5 <<    0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0
100     << 2   <<   4.2 <<    3.6 <<    3.3 <<    1.6 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0
101     << 3   <<   6.6 <<    6.9 <<    7.4 <<    6.7 <<     4.6 <<     2.7 <<     0.4 <<     0.0 <<     0.0 <<     0.0 <<     0.0 <<     0.0
102     << 4   <<   8.6 <<    9.6 <<   10.6 <<   10.1 <<     8.0 <<     6.6 <<     5.0 <<     4.2 <<     2.7 <<     0.0 <<     0.0 <<     0.0
103     << 5   <<  11.8 <<   13.0 <<   16.0 <<   15.1 <<    11.6 <<     9.7 <<     8.1 <<     8.2 <<     7.9 <<     4.9 <<     3.2 <<     2.1
104     << 6   <<  15.6 <<   17.6 <<   23.0 <<   23.6 <<    22.1 <<    20.0 <<    16.0 <<    15.1 <<    12.1 <<     7.9 <<     6.2 <<     5.1
105     << 7   <<  18.7 <<   21.5 <<   28.4 <<   30.2 <<    30.7 <<    31.0 <<    25.2 <<    23.1 <<    17.5 <<    10.7 <<     8.4 <<     7.2;
106
107   bind();
108   Debug(0);
109 }
110
111 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
112
113 FGWinds::~FGWinds()
114 {
115   delete(POE_Table);
116   Debug(1);
117 }
118
119 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
120
121 bool FGWinds::InitModel(void)
122 {
123   return true;
124 }
125
126 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127
128 bool FGWinds::Run(bool Holding)
129 {
130   if (FGModel::Run(Holding)) return true;
131   if (Holding) return false;
132
133   RunPreFunctions();
134
135   if (turbType != ttNone) Turbulence(in.AltitudeASL);
136   if (oneMinusCosineGust.gustProfile.Running) CosineGust();
137
138   vTotalWindNED = vWindNED + vGustNED + vCosineGust + vTurbulenceNED;
139
140    // psiw (Wind heading) is the direction the wind is blowing towards
141   if (vWindNED(eX) != 0.0) psiw = atan2( vWindNED(eY), vWindNED(eX) );
142   if (psiw < 0) psiw += 2*M_PI;
143
144   RunPostFunctions();
145
146   Debug(2);
147   return false;
148 }
149
150 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151 //
152 // psi is the angle that the wind is blowing *towards*
153
154 void FGWinds::SetWindspeed(double speed)
155 {
156   if (vWindNED.Magnitude() == 0.0) {
157     psiw = 0.0;
158     vWindNED(eNorth) = speed;
159   } else {
160     vWindNED(eNorth) = speed * cos(psiw);
161     vWindNED(eEast) = speed * sin(psiw);
162     vWindNED(eDown) = 0.0;
163   }
164 }
165
166 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
167
168 double FGWinds::GetWindspeed(void) const
169 {
170   return vWindNED.Magnitude();
171 }
172
173 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
174 //
175 // psi is the angle that the wind is blowing *towards*
176
177 void FGWinds::SetWindPsi(double dir)
178 {
179   double mag = GetWindspeed();
180   psiw = dir;
181   SetWindspeed(mag);  
182 }
183
184 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
185
186 void FGWinds::Turbulence(double h)
187 {
188   switch (turbType) {
189
190   case ttCulp: { 
191
192     vTurbPQR(eP) = wind_from_clockwise;
193     if (TurbGain == 0.0) return;
194   
195     // keep the inputs within allowable limts for this model
196     if (TurbGain < 0.0) TurbGain = 0.0;
197     if (TurbGain > 1.0) TurbGain = 1.0;
198     if (TurbRate < 0.0) TurbRate = 0.0;
199     if (TurbRate > 30.0) TurbRate = 30.0;
200     if (Rhythmicity < 0.0) Rhythmicity = 0.0;
201     if (Rhythmicity > 1.0) Rhythmicity = 1.0;
202
203     // generate a sine wave corresponding to turbulence rate in hertz
204     double time = FDMExec->GetSimTime();
205     double sinewave = sin( time * TurbRate * 6.283185307 );
206
207     double random = 0.0;
208     if (target_time == 0.0) {
209       strength = random = 1 - 2.0*(double(rand())/double(RAND_MAX));
210       target_time = time + 0.71 + (random * 0.5);
211     }
212     if (time > target_time) {
213       spike = 1.0;
214       target_time = 0.0;
215     }    
216
217     // max vertical wind speed in fps, corresponds to TurbGain = 1.0
218     double max_vs = 40;
219
220     vTurbulenceNED(1) = vTurbulenceNED(2) = vTurbulenceNED(3) = 0.0;
221     double delta = strength * max_vs * TurbGain * (1-Rhythmicity) * spike;
222
223     // Vertical component of turbulence.
224     vTurbulenceNED(3) = sinewave * max_vs * TurbGain * Rhythmicity;
225     vTurbulenceNED(3)+= delta;
226     if (in.DistanceAGL/in.wingspan < 3.0)
227         vTurbulenceNED(3) *= in.DistanceAGL/in.wingspan * 0.3333;
228  
229     // Yaw component of turbulence.
230     vTurbulenceNED(1) = sin( delta * 3.0 );
231     vTurbulenceNED(2) = cos( delta * 3.0 );
232
233     // Roll component of turbulence. Clockwise vortex causes left roll.
234     vTurbPQR(eP) += delta * 0.04;
235
236     spike = spike * 0.9;
237     break;
238   }
239   case ttMilspec:
240   case ttTustin: {
241
242     // an index of zero means turbulence is disabled
243     // airspeed occurs as divisor in the code below
244     if (probability_of_exceedence_index == 0 || in.V == 0) {
245       vTurbulenceNED(1) = vTurbulenceNED(2) = vTurbulenceNED(3) = 0.0;
246       vTurbPQR(1) = vTurbPQR(2) = vTurbPQR(3) = 0.0;
247       return;
248     }
249
250     // Turbulence model according to MIL-F-8785C (Flying Qualities of Piloted Aircraft)
251     double b_w = in.wingspan, L_u, L_w, sig_u, sig_w;
252
253       if (b_w == 0.) b_w = 30.;
254
255     // clip height functions at 10 ft
256     if (h <= 10.) h = 10;
257
258     // Scale lengths L and amplitudes sigma as function of height
259     if (h <= 1000) {
260       L_u = h/pow(0.177 + 0.000823*h, 1.2); // MIL-F-8785c, Fig. 10, p. 55
261       L_w = h;
262       sig_w = 0.1*windspeed_at_20ft;
263       sig_u = sig_w/pow(0.177 + 0.000823*h, 0.4); // MIL-F-8785c, Fig. 11, p. 56
264     } else if (h <= 2000) {
265       // linear interpolation between low altitude and high altitude models
266       L_u = L_w = 1000 + (h-1000.)/1000.*750.;
267       sig_u = sig_w = 0.1*windspeed_at_20ft
268                     + (h-1000.)/1000.*(POE_Table->GetValue(probability_of_exceedence_index, h) - 0.1*windspeed_at_20ft);
269     } else {
270       L_u = L_w = 1750.; //  MIL-F-8785c, Sec. 3.7.2.1, p. 48
271       sig_u = sig_w = POE_Table->GetValue(probability_of_exceedence_index, h);
272     }
273
274     // keep values from last timesteps
275     // TODO maybe use deque?
276     static double
277       xi_u_km1 = 0, nu_u_km1 = 0,
278       xi_v_km1 = 0, xi_v_km2 = 0, nu_v_km1 = 0, nu_v_km2 = 0,
279       xi_w_km1 = 0, xi_w_km2 = 0, nu_w_km1 = 0, nu_w_km2 = 0,
280       xi_p_km1 = 0, nu_p_km1 = 0,
281       xi_q_km1 = 0, xi_r_km1 = 0;
282
283
284     double
285       T_V = in.totalDeltaT, // for compatibility of nomenclature
286       sig_p = 1.9/sqrt(L_w*b_w)*sig_w, // Yeager1998, eq. (8)
287       //sig_q = sqrt(M_PI/2/L_w/b_w), // eq. (14)
288       //sig_r = sqrt(2*M_PI/3/L_w/b_w), // eq. (17)
289       L_p = sqrt(L_w*b_w)/2.6, // eq. (10)
290       tau_u = L_u/in.V, // eq. (6)
291       tau_w = L_w/in.V, // eq. (3)
292       tau_p = L_p/in.V, // eq. (9)
293       tau_q = 4*b_w/M_PI/in.V, // eq. (13)
294       tau_r =3*b_w/M_PI/in.V, // eq. (17)
295       nu_u = GaussianRandomNumber(),
296       nu_v = GaussianRandomNumber(),
297       nu_w = GaussianRandomNumber(),
298       nu_p = GaussianRandomNumber(),
299       xi_u=0, xi_v=0, xi_w=0, xi_p=0, xi_q=0, xi_r=0;
300
301     // values of turbulence NED velocities
302
303     if (turbType == ttTustin) {
304       // the following is the Tustin formulation of Yeager's report
305       double
306         omega_w = in.V/L_w, // hidden in nomenclature p. 3
307         omega_v = in.V/L_u, // this is defined nowhere
308         C_BL  = 1/tau_u/tan(T_V/2/tau_u), // eq. (19)
309         C_BLp = 1/tau_p/tan(T_V/2/tau_p), // eq. (22)
310         C_BLq = 1/tau_q/tan(T_V/2/tau_q), // eq. (24)
311         C_BLr = 1/tau_r/tan(T_V/2/tau_r); // eq. (26)
312
313       // all values calculated so far are strictly positive, except for
314       // the random numbers nu_*. This means that in the code below, all
315       // divisors are strictly positive, too, and no floating point
316       // exception should occur.
317       xi_u = -(1 - C_BL*tau_u)/(1 + C_BL*tau_u)*xi_u_km1
318            + sig_u*sqrt(2*tau_u/T_V)/(1 + C_BL*tau_u)*(nu_u + nu_u_km1); // eq. (18)
319       xi_v = -2*(sqr(omega_v) - sqr(C_BL))/sqr(omega_v + C_BL)*xi_v_km1
320            - sqr(omega_v - C_BL)/sqr(omega_v + C_BL) * xi_v_km2
321            + sig_u*sqrt(3*omega_v/T_V)/sqr(omega_v + C_BL)*(
322                  (C_BL + omega_v/sqrt(3.))*nu_v
323                + 2/sqrt(3.)*omega_v*nu_v_km1
324                + (omega_v/sqrt(3.) - C_BL)*nu_v_km2); // eq. (20) for v
325       xi_w = -2*(sqr(omega_w) - sqr(C_BL))/sqr(omega_w + C_BL)*xi_w_km1
326            - sqr(omega_w - C_BL)/sqr(omega_w + C_BL) * xi_w_km2
327            + sig_w*sqrt(3*omega_w/T_V)/sqr(omega_w + C_BL)*(
328                  (C_BL + omega_w/sqrt(3.))*nu_w
329                + 2/sqrt(3.)*omega_w*nu_w_km1
330                + (omega_w/sqrt(3.) - C_BL)*nu_w_km2); // eq. (20) for w
331       xi_p = -(1 - C_BLp*tau_p)/(1 + C_BLp*tau_p)*xi_p_km1
332            + sig_p*sqrt(2*tau_p/T_V)/(1 + C_BLp*tau_p) * (nu_p + nu_p_km1); // eq. (21)
333       xi_q = -(1 - 4*b_w*C_BLq/M_PI/in.V)/(1 + 4*b_w*C_BLq/M_PI/in.V) * xi_q_km1
334            + C_BLq/in.V/(1 + 4*b_w*C_BLq/M_PI/in.V) * (xi_w - xi_w_km1); // eq. (23)
335       xi_r = - (1 - 3*b_w*C_BLr/M_PI/in.V)/(1 + 3*b_w*C_BLr/M_PI/in.V) * xi_r_km1
336            + C_BLr/in.V/(1 + 3*b_w*C_BLr/M_PI/in.V) * (xi_v - xi_v_km1); // eq. (25)
337
338     } else if (turbType == ttMilspec) {
339       // the following is the MIL-STD-1797A formulation
340       // as cited in Yeager's report
341       xi_u = (1 - T_V/tau_u)  *xi_u_km1 + sig_u*sqrt(2*T_V/tau_u)*nu_u;  // eq. (30)
342       xi_v = (1 - 2*T_V/tau_u)*xi_v_km1 + sig_u*sqrt(4*T_V/tau_u)*nu_v;  // eq. (31)
343       xi_w = (1 - 2*T_V/tau_w)*xi_w_km1 + sig_w*sqrt(4*T_V/tau_w)*nu_w;  // eq. (32)
344       xi_p = (1 - T_V/tau_p)  *xi_p_km1 + sig_p*sqrt(2*T_V/tau_p)*nu_p;  // eq. (33)
345       xi_q = (1 - T_V/tau_q)  *xi_q_km1 + M_PI/4/b_w*(xi_w - xi_w_km1);  // eq. (34)
346       xi_r = (1 - T_V/tau_r)  *xi_r_km1 + M_PI/3/b_w*(xi_v - xi_v_km1);  // eq. (35)
347     }
348
349     // rotate by wind azimuth and assign the velocities
350     double cospsi = cos(psiw), sinpsi = sin(psiw);
351     vTurbulenceNED(1) =  cospsi*xi_u + sinpsi*xi_v;
352     vTurbulenceNED(2) = -sinpsi*xi_u + cospsi*xi_v;
353     vTurbulenceNED(3) = xi_w;
354
355     vTurbPQR(1) =  cospsi*xi_p + sinpsi*xi_q;
356     vTurbPQR(2) = -sinpsi*xi_p + cospsi*xi_q;
357     vTurbPQR(3) = xi_r;
358
359     // vTurbPQR is in the body fixed frame, not NED
360     vTurbPQR = in.Tl2b*vTurbPQR;
361
362     // hand on the values for the next timestep
363     xi_u_km1 = xi_u; nu_u_km1 = nu_u;
364     xi_v_km2 = xi_v_km1; xi_v_km1 = xi_v; nu_v_km2 = nu_v_km1; nu_v_km1 = nu_v;
365     xi_w_km2 = xi_w_km1; xi_w_km1 = xi_w; nu_w_km2 = nu_w_km1; nu_w_km1 = nu_w;
366     xi_p_km1 = xi_p; nu_p_km1 = nu_p;
367     xi_q_km1 = xi_q;
368     xi_r_km1 = xi_r;
369
370   }
371   default:
372     break;
373   }
374 }
375
376 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
377
378 double FGWinds::CosineGustProfile(double startDuration, double steadyDuration, double endDuration, double elapsedTime)
379 {
380   double factor = 0.0;
381   if (elapsedTime >= 0 && elapsedTime <= startDuration) {
382     factor = (1.0 - cos(M_PI*elapsedTime/startDuration))/2.0;
383   } else if (elapsedTime > startDuration && (elapsedTime <= (startDuration + steadyDuration))) {
384     factor = 1.0;
385   } else if (elapsedTime > (startDuration + steadyDuration) && elapsedTime <= (startDuration + steadyDuration + endDuration)) {
386     factor = (1-cos(M_PI*(1-(elapsedTime-(startDuration + steadyDuration))/endDuration)))/2.0;
387   } else {
388     factor = 0.0;
389   }
390
391   return factor;
392 }
393
394 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
395
396 void FGWinds::CosineGust()
397 {
398   struct OneMinusCosineProfile& profile = oneMinusCosineGust.gustProfile;
399
400   double factor = CosineGustProfile( profile.startupDuration,
401                                      profile.steadyDuration,
402                                      profile.endDuration,
403                                      profile.elapsedTime);
404   // Normalize the gust wind vector
405   oneMinusCosineGust.vWind.Normalize();
406
407   if (oneMinusCosineGust.vWindTransformed.Magnitude() == 0.0) {
408     switch (oneMinusCosineGust.gustFrame) {
409     case gfBody:
410       oneMinusCosineGust.vWindTransformed = in.Tl2b.Inverse() * oneMinusCosineGust.vWind;
411       break;
412     case gfWind:
413       oneMinusCosineGust.vWindTransformed = in.Tl2b.Inverse() * in.Tw2b * oneMinusCosineGust.vWind;
414       break;
415     case gfLocal:
416       // this is the native frame - and the default.
417       oneMinusCosineGust.vWindTransformed = oneMinusCosineGust.vWind;
418       break;
419     default:
420       break;
421     }
422   }
423
424   vCosineGust = factor * oneMinusCosineGust.vWindTransformed * oneMinusCosineGust.magnitude;
425
426   profile.elapsedTime += in.totalDeltaT;
427
428   if (profile.elapsedTime > (profile.startupDuration + profile.steadyDuration + profile.endDuration)) {
429     profile.Running = false;
430     profile.elapsedTime = 0.0;
431     oneMinusCosineGust.vWindTransformed.InitMatrix(0.0);
432     vCosineGust.InitMatrix(0);
433   }
434 }
435
436 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
437
438 void FGWinds::bind(void)
439 {
440   typedef double (FGWinds::*PMF)(int) const;
441   typedef int (FGWinds::*PMFt)(void) const;
442   typedef void   (FGWinds::*PMFd)(int,double);
443   typedef void   (FGWinds::*PMFi)(int);
444   typedef double (FGWinds::*Ptr)(void) const;
445
446   // User-specified steady, constant, wind properties (local navigational/geographic frame: N-E-D)
447   PropertyManager->Tie("atmosphere/psiw-rad", this, &FGWinds::GetWindPsi, &FGWinds::SetWindPsi);
448   PropertyManager->Tie("atmosphere/wind-north-fps", this, eNorth, (PMF)&FGWinds::GetWindNED,
449                                                           (PMFd)&FGWinds::SetWindNED);
450   PropertyManager->Tie("atmosphere/wind-east-fps",  this, eEast, (PMF)&FGWinds::GetWindNED,
451                                                           (PMFd)&FGWinds::SetWindNED);
452   PropertyManager->Tie("atmosphere/wind-down-fps",  this, eDown, (PMF)&FGWinds::GetWindNED,
453                                                           (PMFd)&FGWinds::SetWindNED);
454   PropertyManager->Tie("atmosphere/wind-mag-fps", this, &FGWinds::GetWindspeed,
455                                                         &FGWinds::SetWindspeed);
456
457   // User-specifieded gust (local navigational/geographic frame: N-E-D)
458   PropertyManager->Tie("atmosphere/gust-north-fps", this, eNorth, (PMF)&FGWinds::GetGustNED,
459                                                           (PMFd)&FGWinds::SetGustNED);
460   PropertyManager->Tie("atmosphere/gust-east-fps",  this, eEast, (PMF)&FGWinds::GetGustNED,
461                                                           (PMFd)&FGWinds::SetGustNED);
462   PropertyManager->Tie("atmosphere/gust-down-fps",  this, eDown, (PMF)&FGWinds::GetGustNED,
463                                                           (PMFd)&FGWinds::SetGustNED);
464   
465   // User-specified 1 - cosine gust parameters (in specified frame)
466   PropertyManager->Tie("atmosphere/cosine-gust/startup-duration-sec", this, (Ptr)0L, &FGWinds::StartupGustDuration);
467   PropertyManager->Tie("atmosphere/cosine-gust/steady-duration-sec", this, (Ptr)0L, &FGWinds::SteadyGustDuration);
468   PropertyManager->Tie("atmosphere/cosine-gust/end-duration-sec", this, (Ptr)0L, &FGWinds::EndGustDuration);
469   PropertyManager->Tie("atmosphere/cosine-gust/magnitude-ft_sec", this, (Ptr)0L, &FGWinds::GustMagnitude);
470   PropertyManager->Tie("atmosphere/cosine-gust/frame", this, (PMFt)0L, (PMFi)&FGWinds::GustFrame);
471   PropertyManager->Tie("atmosphere/cosine-gust/X-velocity-ft_sec", this, (Ptr)0L, &FGWinds::GustXComponent);
472   PropertyManager->Tie("atmosphere/cosine-gust/Y-velocity-ft_sec", this, (Ptr)0L, &FGWinds::GustYComponent);
473   PropertyManager->Tie("atmosphere/cosine-gust/Z-velocity-ft_sec", this, (Ptr)0L, &FGWinds::GustZComponent);
474   PropertyManager->Tie("atmosphere/cosine-gust/start", this, (PMFt)0L, (PMFi)&FGWinds::StartGust);
475
476   // User-specified turbulence (local navigational/geographic frame: N-E-D)
477   PropertyManager->Tie("atmosphere/turb-north-fps", this, eNorth, (PMF)&FGWinds::GetTurbNED,
478                                                           (PMFd)&FGWinds::SetTurbNED);
479   PropertyManager->Tie("atmosphere/turb-east-fps",  this, eEast, (PMF)&FGWinds::GetTurbNED,
480                                                           (PMFd)&FGWinds::SetTurbNED);
481   PropertyManager->Tie("atmosphere/turb-down-fps",  this, eDown, (PMF)&FGWinds::GetTurbNED,
482                                                           (PMFd)&FGWinds::SetTurbNED);
483   // Experimental turbulence parameters
484   PropertyManager->Tie("atmosphere/p-turb-rad_sec", this,1, (PMF)&FGWinds::GetTurbPQR);
485   PropertyManager->Tie("atmosphere/q-turb-rad_sec", this,2, (PMF)&FGWinds::GetTurbPQR);
486   PropertyManager->Tie("atmosphere/r-turb-rad_sec", this,3, (PMF)&FGWinds::GetTurbPQR);
487   PropertyManager->Tie("atmosphere/turb-type", this, (PMFt)&FGWinds::GetTurbType, (PMFi)&FGWinds::SetTurbType);
488   PropertyManager->Tie("atmosphere/turb-rate", this, &FGWinds::GetTurbRate, &FGWinds::SetTurbRate);
489   PropertyManager->Tie("atmosphere/turb-gain", this, &FGWinds::GetTurbGain, &FGWinds::SetTurbGain);
490   PropertyManager->Tie("atmosphere/turb-rhythmicity", this, &FGWinds::GetRhythmicity,
491                                                             &FGWinds::SetRhythmicity);
492
493   // Parameters for milspec turbulence
494   PropertyManager->Tie("atmosphere/turbulence/milspec/windspeed_at_20ft_AGL-fps",
495                        this, &FGWinds::GetWindspeed20ft,
496                              &FGWinds::SetWindspeed20ft);
497   PropertyManager->Tie("atmosphere/turbulence/milspec/severity",
498                        this, &FGWinds::GetProbabilityOfExceedence,
499                              &FGWinds::SetProbabilityOfExceedence);
500
501   // Total, calculated winds (local navigational/geographic frame: N-E-D). Read only.
502   PropertyManager->Tie("atmosphere/total-wind-north-fps", this, eNorth, (PMF)&FGWinds::GetTotalWindNED);
503   PropertyManager->Tie("atmosphere/total-wind-east-fps",  this, eEast,  (PMF)&FGWinds::GetTotalWindNED);
504   PropertyManager->Tie("atmosphere/total-wind-down-fps",  this, eDown,  (PMF)&FGWinds::GetTotalWindNED);
505
506 }
507
508 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
509 //    The bitmasked value choices are as follows:
510 //    unset: In this case (the default) JSBSim would only print
511 //       out the normally expected messages, essentially echoing
512 //       the config files as they are read. If the environment
513 //       variable is not set, debug_lvl is set to 1 internally
514 //    0: This requests JSBSim not to output any messages
515 //       whatsoever.
516 //    1: This value explicity requests the normal JSBSim
517 //       startup messages
518 //    2: This value asks for a message to be printed out when
519 //       a class is instantiated
520 //    4: When this value is set, a message is displayed when a
521 //       FGModel object executes its Run() method
522 //    8: When this value is set, various runtime state variables
523 //       are printed out periodically
524 //    16: When set various parameters are sanity checked and
525 //       a message is printed out when they go out of bounds
526
527 void FGWinds::Debug(int from)
528 {
529   if (debug_lvl <= 0) return;
530
531   if (debug_lvl & 1) { // Standard console startup message output
532     if (from == 0) { // Constructor
533     }
534   }
535   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
536     if (from == 0) cout << "Instantiated: FGWinds" << endl;
537     if (from == 1) cout << "Destroyed:    FGWinds" << endl;
538   }
539   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
540   }
541   if (debug_lvl & 8 ) { // Runtime state variables
542   }
543   if (debug_lvl & 16) { // Sanity checking
544   }
545   if (debug_lvl & 128) { // 
546   }
547   if (debug_lvl & 64) {
548     if (from == 0) { // Constructor
549       cout << IdSrc << endl;
550       cout << IdHdr << endl;
551     }
552   }
553 }
554
555 } // namespace JSBSim