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.7 2011/12/11 17:03:05 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   if (turbType != ttNone) Turbulence(in.AltitudeASL);
 134   if (oneMinusCosineGust.gustProfile.Running) CosineGust();
 135
 136   vTotalWindNED = vWindNED + vGustNED + vCosineGust + vTurbulenceNED;
 137
 138    // psiw (Wind heading) is the direction the wind is blowing towards
 139   if (vWindNED(eX) != 0.0) psiw = atan2( vWindNED(eY), vWindNED(eX) );
 140   if (psiw < 0) psiw += 2*M_PI;
 141
 142   Debug(2);
 143   return false;
 144 }
 145
 146 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 147 //
 148 // psi is the angle that the wind is blowing *towards*
 149
 150 void FGWinds::SetWindspeed(double speed)
 151 {
 152   if (vWindNED.Magnitude() == 0.0) {
 153     psiw = 0.0;
 154     vWindNED(eNorth) = speed;
 155   } else {
 156     vWindNED(eNorth) = speed * cos(psiw);
 157     vWindNED(eEast) = speed * sin(psiw);
 158     vWindNED(eDown) = 0.0;
 159   }
 160 }
 161
 162 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 163
 164 double FGWinds::GetWindspeed(void) const
 165 {
 166   return vWindNED.Magnitude();
 167 }
 168
 169 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 170 //
 171 // psi is the angle that the wind is blowing *towards*
 172
 173 void FGWinds::SetWindPsi(double dir)
 174 {
 175   double mag = GetWindspeed();
 176   psiw = dir;
 177   SetWindspeed(mag);
 178 }
 179
 180 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 181
 182 void FGWinds::Turbulence(double h)
 183 {
 184   switch (turbType) {
 185
 186   case ttCulp: {
 187
 188     vTurbPQR(eP) = wind_from_clockwise;
 189     if (TurbGain == 0.0) return;
 190
 191     // keep the inputs within allowable limts for this model
 192     if (TurbGain < 0.0) TurbGain = 0.0;
 193     if (TurbGain > 1.0) TurbGain = 1.0;
 194     if (TurbRate < 0.0) TurbRate = 0.0;
 195     if (TurbRate > 30.0) TurbRate = 30.0;
 196     if (Rhythmicity < 0.0) Rhythmicity = 0.0;
 197     if (Rhythmicity > 1.0) Rhythmicity = 1.0;
 198
 199     // generate a sine wave corresponding to turbulence rate in hertz
 200     double time = FDMExec->GetSimTime();
 201     double sinewave = sin( time * TurbRate * 6.283185307 );
 202
 203     double random = 0.0;
 204     if (target_time == 0.0) {
 205       strength = random = 1 - 2.0*(double(rand())/double(RAND_MAX));
 206       target_time = time + 0.71 + (random * 0.5);
 207     }
 208     if (time > target_time) {
 209       spike = 1.0;
 210       target_time = 0.0;
 211     }
 212
 213     // max vertical wind speed in fps, corresponds to TurbGain = 1.0
 214     double max_vs = 40;
 215
 216     vTurbulenceNED(1) = vTurbulenceNED(2) = vTurbulenceNED(3) = 0.0;
 217     double delta = strength * max_vs * TurbGain * (1-Rhythmicity) * spike;
 218
 219     // Vertical component of turbulence.
 220     vTurbulenceNED(3) = sinewave * max_vs * TurbGain * Rhythmicity;
 221     vTurbulenceNED(3)+= delta;
 222     if (in.DistanceAGL/in.wingspan < 3.0)
 223         vTurbulenceNED(3) *= in.DistanceAGL/in.wingspan * 0.3333;
 224
 225     // Yaw component of turbulence.
 226     vTurbulenceNED(1) = sin( delta * 3.0 );
 227     vTurbulenceNED(2) = cos( delta * 3.0 );
 228
 229     // Roll component of turbulence. Clockwise vortex causes left roll.
 230     vTurbPQR(eP) += delta * 0.04;
 231
 232     spike = spike * 0.9;
 233     break;
 234   }
 235   case ttMilspec:
 236   case ttTustin: {
 237
 238     // an index of zero means turbulence is disabled
 239     // airspeed occurs as divisor in the code below
 240     if (probability_of_exceedence_index == 0 || in.V == 0) {
 241       vTurbulenceNED(1) = vTurbulenceNED(2) = vTurbulenceNED(3) = 0.0;
 242       vTurbPQR(1) = vTurbPQR(2) = vTurbPQR(3) = 0.0;
 243       return;
 244     }
 245
 246     // Turbulence model according to MIL-F-8785C (Flying Qualities of Piloted Aircraft)
 247     double b_w = in.wingspan, L_u, L_w, sig_u, sig_w;
 248
 249       if (b_w == 0.) b_w = 30.;
 250
 251     // clip height functions at 10 ft
 252     if (h <= 10.) h = 10;
 253
 254     // Scale lengths L and amplitudes sigma as function of height
 255     if (h <= 1000) {
 256       L_u = h/pow(0.177 + 0.000823*h, 1.2); // MIL-F-8785c, Fig. 10, p. 55
 257       L_w = h;
 258       sig_w = 0.1*windspeed_at_20ft;
 259       sig_u = sig_w/pow(0.177 + 0.000823*h, 0.4); // MIL-F-8785c, Fig. 11, p. 56
 260     } else if (h <= 2000) {
 261       // linear interpolation between low altitude and high altitude models
 262       L_u = L_w = 1000 + (h-1000.)/1000.*750.;
 263       sig_u = sig_w = 0.1*windspeed_at_20ft
 264                     + (h-1000.)/1000.*(POE_Table->GetValue(probability_of_exceedence_index, h) - 0.1*windspeed_at_20ft);
 265     } else {
 266       L_u = L_w = 1750.; //  MIL-F-8785c, Sec. 3.7.2.1, p. 48
 267       sig_u = sig_w = POE_Table->GetValue(probability_of_exceedence_index, h);
 268     }
 269
 270     // keep values from last timesteps
 271     // TODO maybe use deque?
 272     static double
 273       xi_u_km1 = 0, nu_u_km1 = 0,
 274       xi_v_km1 = 0, xi_v_km2 = 0, nu_v_km1 = 0, nu_v_km2 = 0,
 275       xi_w_km1 = 0, xi_w_km2 = 0, nu_w_km1 = 0, nu_w_km2 = 0,
 276       xi_p_km1 = 0, nu_p_km1 = 0,
 277       xi_q_km1 = 0, xi_r_km1 = 0;
 278
 279
 280     double
 281       T_V = in.totalDeltaT, // for compatibility of nomenclature
 282       sig_p = 1.9/sqrt(L_w*b_w)*sig_w, // Yeager1998, eq. (8)
 283       //sig_q = sqrt(M_PI/2/L_w/b_w), // eq. (14)
 284       //sig_r = sqrt(2*M_PI/3/L_w/b_w), // eq. (17)
 285       L_p = sqrt(L_w*b_w)/2.6, // eq. (10)
 286       tau_u = L_u/in.V, // eq. (6)
 287       tau_w = L_w/in.V, // eq. (3)
 288       tau_p = L_p/in.V, // eq. (9)
 289       tau_q = 4*b_w/M_PI/in.V, // eq. (13)
 290       tau_r =3*b_w/M_PI/in.V, // eq. (17)
 291       nu_u = GaussianRandomNumber(),
 292       nu_v = GaussianRandomNumber(),
 293       nu_w = GaussianRandomNumber(),
 294       nu_p = GaussianRandomNumber(),
 295       xi_u=0, xi_v=0, xi_w=0, xi_p=0, xi_q=0, xi_r=0;
 296
 297     // values of turbulence NED velocities
 298
 299     if (turbType == ttTustin) {
 300       // the following is the Tustin formulation of Yeager's report
 301       double
 302         omega_w = in.V/L_w, // hidden in nomenclature p. 3
 303         omega_v = in.V/L_u, // this is defined nowhere
 304         C_BL  = 1/tau_u/tan(T_V/2/tau_u), // eq. (19)
 305         C_BLp = 1/tau_p/tan(T_V/2/tau_p), // eq. (22)
 306         C_BLq = 1/tau_q/tan(T_V/2/tau_q), // eq. (24)
 307         C_BLr = 1/tau_r/tan(T_V/2/tau_r); // eq. (26)
 308
 309       // all values calculated so far are strictly positive, except for
 310       // the random numbers nu_*. This means that in the code below, all
 311       // divisors are strictly positive, too, and no floating point
 312       // exception should occur.
 313       xi_u = -(1 - C_BL*tau_u)/(1 + C_BL*tau_u)*xi_u_km1
 314            + sig_u*sqrt(2*tau_u/T_V)/(1 + C_BL*tau_u)*(nu_u + nu_u_km1); // eq. (18)
 315       xi_v = -2*(sqr(omega_v) - sqr(C_BL))/sqr(omega_v + C_BL)*xi_v_km1
 316            - sqr(omega_v - C_BL)/sqr(omega_v + C_BL) * xi_v_km2
 317            + sig_u*sqrt(3*omega_v/T_V)/sqr(omega_v + C_BL)*(
 318                  (C_BL + omega_v/sqrt(3.))*nu_v
 319                + 2/sqrt(3.)*omega_v*nu_v_km1
 320                + (omega_v/sqrt(3.) - C_BL)*nu_v_km2); // eq. (20) for v
 321       xi_w = -2*(sqr(omega_w) - sqr(C_BL))/sqr(omega_w + C_BL)*xi_w_km1
 322            - sqr(omega_w - C_BL)/sqr(omega_w + C_BL) * xi_w_km2
 323            + sig_w*sqrt(3*omega_w/T_V)/sqr(omega_w + C_BL)*(
 324                  (C_BL + omega_w/sqrt(3.))*nu_w
 325                + 2/sqrt(3.)*omega_w*nu_w_km1
 326                + (omega_w/sqrt(3.) - C_BL)*nu_w_km2); // eq. (20) for w
 327       xi_p = -(1 - C_BLp*tau_p)/(1 + C_BLp*tau_p)*xi_p_km1
 328            + sig_p*sqrt(2*tau_p/T_V)/(1 + C_BLp*tau_p) * (nu_p + nu_p_km1); // eq. (21)
 329       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
 330            + C_BLq/in.V/(1 + 4*b_w*C_BLq/M_PI/in.V) * (xi_w - xi_w_km1); // eq. (23)
 331       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
 332            + C_BLr/in.V/(1 + 3*b_w*C_BLr/M_PI/in.V) * (xi_v - xi_v_km1); // eq. (25)
 333
 334     } else if (turbType == ttMilspec) {
 335       // the following is the MIL-STD-1797A formulation
 336       // as cited in Yeager's report
 337       xi_u = (1 - T_V/tau_u)  *xi_u_km1 + sig_u*sqrt(2*T_V/tau_u)*nu_u;  // eq. (30)
 338       xi_v = (1 - 2*T_V/tau_u)*xi_v_km1 + sig_u*sqrt(4*T_V/tau_u)*nu_v;  // eq. (31)
 339       xi_w = (1 - 2*T_V/tau_w)*xi_w_km1 + sig_w*sqrt(4*T_V/tau_w)*nu_w;  // eq. (32)
 340       xi_p = (1 - T_V/tau_p)  *xi_p_km1 + sig_p*sqrt(2*T_V/tau_p)*nu_p;  // eq. (33)
 341       xi_q = (1 - T_V/tau_q)  *xi_q_km1 + M_PI/4/b_w*(xi_w - xi_w_km1);  // eq. (34)
 342       xi_r = (1 - T_V/tau_r)  *xi_r_km1 + M_PI/3/b_w*(xi_v - xi_v_km1);  // eq. (35)
 343     }
 344
 345     // rotate by wind azimuth and assign the velocities
 346     double cospsi = cos(psiw), sinpsi = sin(psiw);
 347     vTurbulenceNED(1) =  cospsi*xi_u + sinpsi*xi_v;
 348     vTurbulenceNED(2) = -sinpsi*xi_u + cospsi*xi_v;
 349     vTurbulenceNED(3) = xi_w;
 350
 351     vTurbPQR(1) =  cospsi*xi_p + sinpsi*xi_q;
 352     vTurbPQR(2) = -sinpsi*xi_p + cospsi*xi_q;
 353     vTurbPQR(3) = xi_r;
 354
 355     // vTurbPQR is in the body fixed frame, not NED
 356     vTurbPQR = in.Tl2b*vTurbPQR;
 357
 358     // hand on the values for the next timestep
 359     xi_u_km1 = xi_u; nu_u_km1 = nu_u;
 360     xi_v_km2 = xi_v_km1; xi_v_km1 = xi_v; nu_v_km2 = nu_v_km1; nu_v_km1 = nu_v;
 361     xi_w_km2 = xi_w_km1; xi_w_km1 = xi_w; nu_w_km2 = nu_w_km1; nu_w_km1 = nu_w;
 362     xi_p_km1 = xi_p; nu_p_km1 = nu_p;
 363     xi_q_km1 = xi_q;
 364     xi_r_km1 = xi_r;
 365
 366   }
 367   default:
 368     break;
 369   }
 370 }
 371
 372 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 373
 374 double FGWinds::CosineGustProfile(double startDuration, double steadyDuration, double endDuration, double elapsedTime)
 375 {
 376   double factor = 0.0;
 377   if (elapsedTime >= 0 && elapsedTime <= startDuration) {
 378     factor = (1.0 - cos(M_PI*elapsedTime/startDuration))/2.0;
 379   } else if (elapsedTime > startDuration && (elapsedTime <= (startDuration + steadyDuration))) {
 380     factor = 1.0;
 381   } else if (elapsedTime > (startDuration + steadyDuration) && elapsedTime <= (startDuration + steadyDuration + endDuration)) {
 382     factor = (1-cos(M_PI*(1-(elapsedTime-(startDuration + steadyDuration))/endDuration)))/2.0;
 383   } else {
 384     factor = 0.0;
 385   }
 386
 387   return factor;
 388 }
 389
 390 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 391
 392 void FGWinds::CosineGust()
 393 {
 394   struct OneMinusCosineProfile& profile = oneMinusCosineGust.gustProfile;
 395
 396   double factor = CosineGustProfile( profile.startupDuration,
 397                                      profile.steadyDuration,
 398                                      profile.endDuration,
 399                                      profile.elapsedTime);
 400   // Normalize the gust wind vector
 401   oneMinusCosineGust.vWind.Normalize();
 402
 403   if (oneMinusCosineGust.vWindTransformed.Magnitude() == 0.0) {
 404     switch (oneMinusCosineGust.gustFrame) {
 405     case gfBody:
 406       oneMinusCosineGust.vWindTransformed = in.Tl2b.Inverse() * oneMinusCosineGust.vWind;
 407       break;
 408     case gfWind:
 409       oneMinusCosineGust.vWindTransformed = in.Tl2b.Inverse() * in.Tw2b * oneMinusCosineGust.vWind;
 410       break;
 411     case gfLocal:
 412       // this is the native frame - and the default.
 413       oneMinusCosineGust.vWindTransformed = oneMinusCosineGust.vWind;
 414       break;
 415     default:
 416       break;
 417     }
 418   }
 419
 420   vCosineGust = factor * oneMinusCosineGust.vWindTransformed * oneMinusCosineGust.magnitude;
 421
 422   profile.elapsedTime += in.totalDeltaT;
 423
 424   if (profile.elapsedTime > (profile.startupDuration + profile.steadyDuration + profile.endDuration)) {
 425     profile.Running = false;
 426     profile.elapsedTime = 0.0;
 427     oneMinusCosineGust.vWindTransformed.InitMatrix(0.0);
 428     vCosineGust.InitMatrix(0);
 429   }
 430 }
 431
 432 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 433
 434 void FGWinds::NumberOfUpDownburstCells(int num)
 435 {
 436   for (unsigned int i=0; i<UpDownBurstCells.size();i++) delete UpDownBurstCells[i];
 437   UpDownBurstCells.clear();
 438   if (num >= 0) {
 439     for (unsigned int i=0; i<num; i++) UpDownBurstCells.push_back(new struct UpDownBurst);
 440   }
 441 }
 442
 443 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 444 // Calculates the distance between a specified point (where presumably the
 445 // Up/Downburst is centered) and the current vehicle location. The distance
 446 // here is calculated from the Haversine formula.
 447
 448 double FGWinds::DistanceFromRingCenter(double lat, double lon)
 449 {
 450   double deltaLat = in.latitude - lat;
 451   double deltaLong = in.longitude - lon;
 452   double dLat2 = deltaLat/2.0;
 453   double dLong2 = deltaLong/2.0;
 454   double a = sin(dLat2)*sin(dLat2)
 455              + cos(lat)*cos(in.latitude)*sin(dLong2)*sin(dLong2);
 456   double c = 2.0*atan2(sqrt(a), sqrt(1.0 - a));
 457   double d = in.planetRadius*c;
 458   return d;
 459 }
 460
 461 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 462
 463 void FGWinds::UpDownBurst()
 464 {
 465
 466   for (unsigned int i=0; i<UpDownBurstCells.size(); i++) {
 467     double d = DistanceFromRingCenter(UpDownBurstCells[i]->ringLatitude, UpDownBurstCells[i]->ringLongitude);
 468
 469   }
 470 }
 471
 472 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 473
 474 void FGWinds::bind(void)
 475 {
 476   typedef double (FGWinds::*PMF)(int) const;
 477   typedef int (FGWinds::*PMFt)(void) const;
 478   typedef void   (FGWinds::*PMFd)(int,double);
 479   typedef void   (FGWinds::*PMFi)(int);
 480   typedef double (FGWinds::*Ptr)(void) const;
 481
 482   // User-specified steady, constant, wind properties (local navigational/geographic frame: N-E-D)
 483   PropertyManager->Tie("atmosphere/psiw-rad", this, &FGWinds::GetWindPsi, &FGWinds::SetWindPsi);
 484   PropertyManager->Tie("atmosphere/wind-north-fps", this, eNorth, (PMF)&FGWinds::GetWindNED,
 485                                                           (PMFd)&FGWinds::SetWindNED);
 486   PropertyManager->Tie("atmosphere/wind-east-fps",  this, eEast, (PMF)&FGWinds::GetWindNED,
 487                                                           (PMFd)&FGWinds::SetWindNED);
 488   PropertyManager->Tie("atmosphere/wind-down-fps",  this, eDown, (PMF)&FGWinds::GetWindNED,
 489                                                           (PMFd)&FGWinds::SetWindNED);
 490   PropertyManager->Tie("atmosphere/wind-mag-fps", this, &FGWinds::GetWindspeed,
 491                                                         &FGWinds::SetWindspeed);
 492
 493   // User-specifieded gust (local navigational/geographic frame: N-E-D)
 494   PropertyManager->Tie("atmosphere/gust-north-fps", this, eNorth, (PMF)&FGWinds::GetGustNED,
 495                                                           (PMFd)&FGWinds::SetGustNED);
 496   PropertyManager->Tie("atmosphere/gust-east-fps",  this, eEast, (PMF)&FGWinds::GetGustNED,
 497                                                           (PMFd)&FGWinds::SetGustNED);
 498   PropertyManager->Tie("atmosphere/gust-down-fps",  this, eDown, (PMF)&FGWinds::GetGustNED,
 499                                                           (PMFd)&FGWinds::SetGustNED);
 500
 501   // User-specified 1 - cosine gust parameters (in specified frame)
 502   PropertyManager->Tie("atmosphere/cosine-gust/startup-duration-sec", this, (Ptr)0L, &FGWinds::StartupGustDuration);
 503   PropertyManager->Tie("atmosphere/cosine-gust/steady-duration-sec", this, (Ptr)0L, &FGWinds::SteadyGustDuration);
 504   PropertyManager->Tie("atmosphere/cosine-gust/end-duration-sec", this, (Ptr)0L, &FGWinds::EndGustDuration);
 505   PropertyManager->Tie("atmosphere/cosine-gust/magnitude-ft_sec", this, (Ptr)0L, &FGWinds::GustMagnitude);
 506   PropertyManager->Tie("atmosphere/cosine-gust/frame", this, (PMFt)0L, (PMFi)&FGWinds::GustFrame);
 507   PropertyManager->Tie("atmosphere/cosine-gust/X-velocity-ft_sec", this, (Ptr)0L, &FGWinds::GustXComponent);
 508   PropertyManager->Tie("atmosphere/cosine-gust/Y-velocity-ft_sec", this, (Ptr)0L, &FGWinds::GustYComponent);
 509   PropertyManager->Tie("atmosphere/cosine-gust/Z-velocity-ft_sec", this, (Ptr)0L, &FGWinds::GustZComponent);
 510   PropertyManager->Tie("atmosphere/cosine-gust/start", this, (PMFt)0L, (PMFi)&FGWinds::StartGust);
 511
 512   // User-specified Up- Down-burst parameters
 513   PropertyManager->Tie("atmosphere/updownburst/number-of-cells", this, (PMFt)0L, &FGWinds::NumberOfUpDownburstCells);
 514 //  PropertyManager->Tie("atmosphere/updownburst/", this, (Ptr)0L, &FGWinds::);
 515 //  PropertyManager->Tie("atmosphere/updownburst/", this, (Ptr)0L, &FGWinds::);
 516 //  PropertyManager->Tie("atmosphere/updownburst/", this, (Ptr)0L, &FGWinds::);
 517 //  PropertyManager->Tie("atmosphere/updownburst/", this, (Ptr)0L, &FGWinds::);
 518 //  PropertyManager->Tie("atmosphere/updownburst/", this, (Ptr)0L, &FGWinds::);
 519 //  PropertyManager->Tie("atmosphere/updownburst/", this, (Ptr)0L, &FGWinds::);
 520 //  PropertyManager->Tie("atmosphere/updownburst/", this, (Ptr)0L, &FGWinds::);
 521
 522   // User-specified turbulence (local navigational/geographic frame: N-E-D)
 523   PropertyManager->Tie("atmosphere/turb-north-fps", this, eNorth, (PMF)&FGWinds::GetTurbNED,
 524                                                           (PMFd)&FGWinds::SetTurbNED);
 525   PropertyManager->Tie("atmosphere/turb-east-fps",  this, eEast, (PMF)&FGWinds::GetTurbNED,
 526                                                           (PMFd)&FGWinds::SetTurbNED);
 527   PropertyManager->Tie("atmosphere/turb-down-fps",  this, eDown, (PMF)&FGWinds::GetTurbNED,
 528                                                           (PMFd)&FGWinds::SetTurbNED);
 529   // Experimental turbulence parameters
 530   PropertyManager->Tie("atmosphere/p-turb-rad_sec", this,1, (PMF)&FGWinds::GetTurbPQR);
 531   PropertyManager->Tie("atmosphere/q-turb-rad_sec", this,2, (PMF)&FGWinds::GetTurbPQR);
 532   PropertyManager->Tie("atmosphere/r-turb-rad_sec", this,3, (PMF)&FGWinds::GetTurbPQR);
 533   PropertyManager->Tie("atmosphere/turb-type", this, (PMFt)&FGWinds::GetTurbType, (PMFi)&FGWinds::SetTurbType);
 534   PropertyManager->Tie("atmosphere/turb-rate", this, &FGWinds::GetTurbRate, &FGWinds::SetTurbRate);
 535   PropertyManager->Tie("atmosphere/turb-gain", this, &FGWinds::GetTurbGain, &FGWinds::SetTurbGain);
 536   PropertyManager->Tie("atmosphere/turb-rhythmicity", this, &FGWinds::GetRhythmicity,
 537                                                             &FGWinds::SetRhythmicity);
 538
 539   // Parameters for milspec turbulence
 540   PropertyManager->Tie("atmosphere/turbulence/milspec/windspeed_at_20ft_AGL-fps",
 541                        this, &FGWinds::GetWindspeed20ft,
 542                              &FGWinds::SetWindspeed20ft);
 543   PropertyManager->Tie("atmosphere/turbulence/milspec/severity",
 544                        this, &FGWinds::GetProbabilityOfExceedence,
 545                              &FGWinds::SetProbabilityOfExceedence);
 546
 547   // Total, calculated winds (local navigational/geographic frame: N-E-D). Read only.
 548   PropertyManager->Tie("atmosphere/total-wind-north-fps", this, eNorth, (PMF)&FGWinds::GetTotalWindNED);
 549   PropertyManager->Tie("atmosphere/total-wind-east-fps",  this, eEast,  (PMF)&FGWinds::GetTotalWindNED);
 550   PropertyManager->Tie("atmosphere/total-wind-down-fps",  this, eDown,  (PMF)&FGWinds::GetTotalWindNED);
 551
 552 }
 553
 554 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 555 //    The bitmasked value choices are as follows:
 556 //    unset: In this case (the default) JSBSim would only print
 557 //       out the normally expected messages, essentially echoing
 558 //       the config files as they are read. If the environment
 559 //       variable is not set, debug_lvl is set to 1 internally
 560 //    0: This requests JSBSim not to output any messages
 561 //       whatsoever.
 562 //    1: This value explicity requests the normal JSBSim
 563 //       startup messages
 564 //    2: This value asks for a message to be printed out when
 565 //       a class is instantiated
 566 //    4: When this value is set, a message is displayed when a
 567 //       FGModel object executes its Run() method
 568 //    8: When this value is set, various runtime state variables
 569 //       are printed out periodically
 570 //    16: When set various parameters are sanity checked and
 571 //       a message is printed out when they go out of bounds
 572
 573 void FGWinds::Debug(int from)
 574 {
 575   if (debug_lvl <= 0) return;
 576
 577   if (debug_lvl & 1) { // Standard console startup message output
 578     if (from == 0) { // Constructor
 579     }
 580   }
 581   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 582     if (from == 0) cout << "Instantiated: FGWinds" << endl;
 583     if (from == 1) cout << "Destroyed:    FGWinds" << endl;
 584   }
 585   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 586   }
 587   if (debug_lvl & 8 ) { // Runtime state variables
 588   }
 589   if (debug_lvl & 16) { // Sanity checking
 590   }
 591   if (debug_lvl & 128) { //
 592   }
 593   if (debug_lvl & 64) {
 594     if (from == 0) { // Constructor
 595       cout << IdSrc << endl;
 596       cout << IdHdr << endl;
 597     }
 598   }
 599 }
 600
 601 } // namespace JSBSim