src/FDM/JSBSim/FGInitialCondition.cpp

   1 /*******************************************************************************
   2
   3  Header:       FGInitialCondition.cpp
   4  Author:       Tony Peden
   5  Date started: 7/1/99
   6
   7  ------------- Copyright (C) 1999  Anthony K. Peden (apeden@earthlink.net) -------------
   8
   9  This program is free software; you can redistribute it and/or modify it under
  10  the terms of the GNU General Public License as published by the Free Software
  11  Foundation; either version 2 of the License, or (at your option) any later
  12  version.
  13
  14  This program is distributed in the hope that it will be useful, but WITHOUT
  15  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  16  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  17  details.
  18
  19  You should have received a copy of the GNU General Public License along with
  20  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  21  Place - Suite 330, Boston, MA  02111-1307, USA.
  22
  23  Further information about the GNU General Public License can also be found on
  24  the world wide web at http://www.gnu.org.
  25
  26
  27  HISTORY
  28 --------------------------------------------------------------------------------
  29 7/1/99   TP   Created
  30
  31
  32 FUNCTIONAL DESCRIPTION
  33 --------------------------------------------------------------------------------
  34
  35 The purpose of this class is to take a set of initial conditions and provide
  36 a kinematically consistent set of body axis velocity components, euler
  37 angles, and altitude.  This class does not attempt to trim the model i.e.
  38 the sim will most likely start in a very dynamic state (unless, of course,
  39 you have chosen your IC's wisely) even after setting it up with this class.
  40
  41 ********************************************************************************
  42 INCLUDES
  43 *******************************************************************************/
  44
  45 #include "FGInitialCondition.h"
  46 #include "FGFDMExec.h"
  47 #include "FGState.h"
  48 #include "FGAtmosphere.h"
  49 #include "FGFCS.h"
  50 #include "FGAircraft.h"
  51 #include "FGTranslation.h"
  52 #include "FGRotation.h"
  53 #include "FGPosition.h"
  54 #include "FGAuxiliary.h"
  55 #include "FGOutput.h"
  56 #include "FGConfigFile.h"
  57
  58 static const char *IdSrc = "$Id$";
  59 static const char *IdHdr = ID_INITIALCONDITION;
  60
  61 //******************************************************************************
  62
  63 FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec)
  64 {
  65   vt=vc=ve=vg=0;
  66   mach=0;
  67   alpha=beta=gamma=0;
  68   theta=phi=psi=0;
  69   altitude=hdot=0;
  70   latitude=longitude=0;
  71   u=v=w=0;
  72   uw=vw=ww=0;
  73   vnorth=veast=vdown=0;
  74   wnorth=weast=wdown=0;
  75   whead=wcross=0;
  76   wdir=wmag=0;
  77   lastSpeedSet=setvt;
  78   lastWindSet=setwned;
  79   sea_level_radius = FDMExec->GetInertial()->RefRadius();
  80   radius_to_vehicle = FDMExec->GetInertial()->RefRadius();
  81   terrain_altitude = 0;
  82
  83   salpha=sbeta=stheta=sphi=spsi=sgamma=0;
  84   calpha=cbeta=ctheta=cphi=cpsi=cgamma=1;
  85
  86   if(FDMExec != NULL ) {
  87     fdmex=FDMExec;
  88     fdmex->GetPosition()->Seth(altitude);
  89     fdmex->GetAtmosphere()->Run();
  90   } else {
  91     cout << "FGInitialCondition: This class requires a pointer to a valid FGFDMExec object" << endl;
  92   }
  93
  94   Debug(0);
  95 }
  96
  97 //******************************************************************************
  98
  99 FGInitialCondition::~FGInitialCondition()
 100 {
 101   Debug(1);
 102 }
 103
 104 //******************************************************************************
 105
 106 void FGInitialCondition::SetVcalibratedKtsIC(double tt) {
 107
 108   if(getMachFromVcas(&mach,tt*ktstofps)) {
 109     //cout << "Mach: " << mach << endl;
 110     lastSpeedSet=setvc;
 111     vc=tt*ktstofps;
 112     vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
 113     ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 114     //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
 115   }
 116   else {
 117     cout << "Failed to get Mach number for given Vc and altitude, Vc unchanged." << endl;
 118     cout << "Please mail the set of initial conditions used to apeden@earthlink.net" << endl;
 119   }
 120 }
 121
 122 //******************************************************************************
 123
 124 void FGInitialCondition::SetVequivalentKtsIC(double tt) {
 125   ve=tt*ktstofps;
 126   lastSpeedSet=setve;
 127   vt=ve*1/sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 128   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 129   vc=calcVcas(mach);
 130 }
 131
 132 //******************************************************************************
 133
 134 void FGInitialCondition::SetVgroundFpsIC(double tt) {
 135   double ua,va,wa;
 136   double vxz;
 137
 138   vg=tt;
 139   lastSpeedSet=setvg;
 140   vnorth = vg*cos(psi); veast = vg*sin(psi); vdown = 0;
 141   calcUVWfromNED();
 142   ua = u + uw; va = v + vw; wa = w + ww;
 143   vt = sqrt( ua*ua + va*va + wa*wa );
 144   alpha = beta = 0;
 145   vxz = sqrt( u*u + w*w );
 146   if( w != 0 ) alpha = atan2( w, u );
 147   if( vxz != 0 ) beta = atan2( v, vxz );
 148   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 149   vc=calcVcas(mach);
 150   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 151 }
 152
 153 //******************************************************************************
 154
 155 void FGInitialCondition::SetVtrueFpsIC(double tt) {
 156   vt=tt;
 157   lastSpeedSet=setvt;
 158   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 159   vc=calcVcas(mach);
 160   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 161 }
 162
 163 //******************************************************************************
 164
 165 void FGInitialCondition::SetMachIC(double tt) {
 166   mach=tt;
 167   lastSpeedSet=setmach;
 168   vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
 169   vc=calcVcas(mach);
 170   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 171   //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
 172 }
 173
 174 //******************************************************************************
 175
 176 void FGInitialCondition::SetClimbRateFpmIC(double tt) {
 177   SetClimbRateFpsIC(tt/60.0);
 178 }
 179
 180 //******************************************************************************
 181
 182 void FGInitialCondition::SetClimbRateFpsIC(double tt) {
 183
 184   if(vt > 0.1) {
 185     hdot=tt;
 186     gamma=asin(hdot/vt);
 187     sgamma=sin(gamma); cgamma=cos(gamma);
 188   }
 189 }
 190
 191 //******************************************************************************
 192
 193 void FGInitialCondition::SetFlightPathAngleRadIC(double tt) {
 194   gamma=tt;
 195   sgamma=sin(gamma); cgamma=cos(gamma);
 196   getTheta();
 197   hdot=vt*sgamma;
 198 }
 199
 200 //******************************************************************************
 201
 202 void FGInitialCondition::SetAlphaRadIC(double tt) {
 203   alpha=tt;
 204   salpha=sin(alpha); calpha=cos(alpha);
 205   getTheta();
 206 }
 207
 208 //******************************************************************************
 209
 210 void FGInitialCondition::SetPitchAngleRadIC(double tt) {
 211   theta=tt;
 212   stheta=sin(theta); ctheta=cos(theta);
 213   getAlpha();
 214 }
 215
 216 //******************************************************************************
 217
 218 void FGInitialCondition::SetBetaRadIC(double tt) {
 219   beta=tt;
 220   sbeta=sin(beta); cbeta=cos(beta);
 221   getTheta();
 222
 223 }
 224
 225 //******************************************************************************
 226
 227 void FGInitialCondition::SetRollAngleRadIC(double tt) {
 228   phi=tt;
 229   sphi=sin(phi); cphi=cos(phi);
 230   getTheta();
 231 }
 232
 233 //******************************************************************************
 234
 235 void FGInitialCondition::SetTrueHeadingRadIC(double tt) {
 236     psi=tt;
 237     spsi=sin(psi); cpsi=cos(psi);
 238     calcWindUVW();
 239 }
 240
 241 //******************************************************************************
 242
 243 void FGInitialCondition::SetUBodyFpsIC(double tt) {
 244   u=tt;
 245   vt=sqrt(u*u + v*v + w*w);
 246   lastSpeedSet=setuvw;
 247 }
 248
 249 //******************************************************************************
 250
 251 void FGInitialCondition::SetVBodyFpsIC(double tt) {
 252   v=tt;
 253   vt=sqrt(u*u + v*v + w*w);
 254   lastSpeedSet=setuvw;
 255 }
 256
 257 //******************************************************************************
 258
 259 void FGInitialCondition::SetWBodyFpsIC(double tt) {
 260   w=tt;
 261   vt=sqrt( u*u + v*v + w*w );
 262   lastSpeedSet=setuvw;
 263 }
 264
 265 //******************************************************************************
 266
 267 double FGInitialCondition::GetUBodyFpsIC(void) {
 268     if(lastSpeedSet == setvg )
 269       return u;
 270     else
 271       return vt*calpha*cbeta - uw;
 272 }
 273
 274 //******************************************************************************
 275
 276 double FGInitialCondition::GetVBodyFpsIC(void) {
 277     if( lastSpeedSet == setvg )
 278       return v;
 279     else {
 280       return vt*sbeta - vw;
 281     }
 282 }
 283
 284 //******************************************************************************
 285
 286 double FGInitialCondition::GetWBodyFpsIC(void) {
 287     if( lastSpeedSet == setvg )
 288       return w;
 289     else
 290       return vt*salpha*cbeta -ww;
 291 }
 292
 293 //******************************************************************************
 294
 295 void FGInitialCondition::SetWindNEDFpsIC(double wN, double wE, double wD ) {
 296   wnorth = wN; weast = wE; wdown = wD;
 297   lastWindSet = setwned;
 298   calcWindUVW();
 299   if(lastSpeedSet == setvg)
 300     SetVgroundFpsIC(vg);
 301 }
 302
 303 //******************************************************************************
 304
 305 // positive from left
 306 void FGInitialCondition::SetHeadWindKtsIC(double head){
 307     whead=head*ktstofps;
 308     lastWindSet=setwhc;
 309     calcWindUVW();
 310     if(lastSpeedSet == setvg)
 311       SetVgroundFpsIC(vg);
 312
 313 }
 314
 315 //******************************************************************************
 316
 317 void FGInitialCondition::SetCrossWindKtsIC(double cross){
 318     wcross=cross*ktstofps;
 319     lastWindSet=setwhc;
 320     calcWindUVW();
 321     if(lastSpeedSet == setvg)
 322       SetVgroundFpsIC(vg);
 323
 324 }
 325
 326 //******************************************************************************
 327
 328 void FGInitialCondition::SetWindDownKtsIC(double wD) {
 329     wdown=wD;
 330     calcWindUVW();
 331     if(lastSpeedSet == setvg)
 332       SetVgroundFpsIC(vg);
 333 }
 334
 335 //******************************************************************************
 336
 337 void FGInitialCondition::SetWindMagKtsIC(double mag) {
 338   wmag=mag*ktstofps;
 339   lastWindSet=setwmd;
 340   calcWindUVW();
 341   if(lastSpeedSet == setvg)
 342       SetVgroundFpsIC(vg);
 343 }
 344
 345 //******************************************************************************
 346
 347 void FGInitialCondition::SetWindDirDegIC(double dir) {
 348   wdir=dir*degtorad;
 349   lastWindSet=setwmd;
 350   calcWindUVW();
 351   if(lastSpeedSet == setvg)
 352       SetVgroundFpsIC(vg);
 353 }
 354
 355
 356 //******************************************************************************
 357
 358 void FGInitialCondition::calcWindUVW(void) {
 359
 360     switch(lastWindSet) {
 361       case setwmd:
 362         wnorth=wmag*cos(wdir);
 363         weast=wmag*sin(wdir);
 364       break;
 365       case setwhc:
 366         wnorth=whead*cos(psi) + wcross*cos(psi+M_PI/2);
 367         weast=whead*sin(psi) + wcross*sin(psi+M_PI/2);
 368       break;
 369       case setwned:
 370       break;
 371     }
 372     uw=wnorth*ctheta*cpsi +
 373        weast*ctheta*spsi -
 374        wdown*stheta;
 375     vw=wnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 376         weast*( sphi*stheta*spsi + cphi*cpsi ) +
 377        wdown*sphi*ctheta;
 378     ww=wnorth*(cphi*stheta*cpsi + sphi*spsi) +
 379        weast*(cphi*stheta*spsi - sphi*cpsi) +
 380        wdown*cphi*ctheta;
 381
 382
 383     /* cout << "FGInitialCondition::calcWindUVW: wnorth, weast, wdown "
 384          << wnorth << ", " << weast << ", " << wdown << endl;
 385     cout << "FGInitialCondition::calcWindUVW: theta, phi, psi "
 386           << theta << ", " << phi << ", " << psi << endl;
 387     cout << "FGInitialCondition::calcWindUVW: uw, vw, ww "
 388           << uw << ", " << vw << ", " << ww << endl;   */
 389
 390 }
 391
 392 //******************************************************************************
 393
 394 void FGInitialCondition::SetAltitudeFtIC(double tt) {
 395   altitude=tt;
 396   fdmex->GetPosition()->Seth(altitude);
 397   fdmex->GetAtmosphere()->Run();
 398   //lets try to make sure the user gets what they intended
 399
 400   switch(lastSpeedSet) {
 401   case setned:
 402   case setuvw:
 403   case setvt:
 404     SetVtrueKtsIC(vt*fpstokts);
 405     break;
 406   case setvc:
 407     SetVcalibratedKtsIC(vc*fpstokts);
 408     break;
 409   case setve:
 410     SetVequivalentKtsIC(ve*fpstokts);
 411     break;
 412   case setmach:
 413     SetMachIC(mach);
 414     break;
 415   case setvg:
 416     SetVgroundFpsIC(vg);
 417     break;
 418   }
 419 }
 420
 421 //******************************************************************************
 422
 423 void FGInitialCondition::SetAltitudeAGLFtIC(double tt) {
 424   fdmex->GetPosition()->SetDistanceAGL(tt);
 425   altitude=fdmex->GetPosition()->Geth();
 426   SetAltitudeFtIC(altitude);
 427 }
 428
 429 //******************************************************************************
 430
 431 void FGInitialCondition::SetSeaLevelRadiusFtIC(double tt) {
 432   sea_level_radius = tt;
 433 }
 434
 435 //******************************************************************************
 436
 437 void FGInitialCondition::SetTerrainAltitudeFtIC(double tt) {
 438   terrain_altitude=tt;
 439 }
 440
 441 //******************************************************************************
 442
 443 void FGInitialCondition::calcUVWfromNED(void) {
 444   u=vnorth*ctheta*cpsi +
 445      veast*ctheta*spsi -
 446      vdown*stheta;
 447   v=vnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 448      veast*( sphi*stheta*spsi + cphi*cpsi ) +
 449      vdown*sphi*ctheta;
 450   w=vnorth*( cphi*stheta*cpsi + sphi*spsi ) +
 451      veast*( cphi*stheta*spsi - sphi*cpsi ) +
 452      vdown*cphi*ctheta;
 453 }
 454
 455 //******************************************************************************
 456
 457 void FGInitialCondition::SetVnorthFpsIC(double tt) {
 458   vnorth=tt;
 459   calcUVWfromNED();
 460   vt=sqrt(u*u + v*v + w*w);
 461   lastSpeedSet=setned;
 462 }
 463
 464 //******************************************************************************
 465
 466 void FGInitialCondition::SetVeastFpsIC(double tt) {
 467   veast=tt;
 468   calcUVWfromNED();
 469   vt=sqrt(u*u + v*v + w*w);
 470   lastSpeedSet=setned;
 471 }
 472
 473 //******************************************************************************
 474
 475 void FGInitialCondition::SetVdownFpsIC(double tt) {
 476   vdown=tt;
 477   calcUVWfromNED();
 478   vt=sqrt(u*u + v*v + w*w);
 479   SetClimbRateFpsIC(-1*vdown);
 480   lastSpeedSet=setned;
 481 }
 482
 483 //******************************************************************************
 484
 485 bool FGInitialCondition::getMachFromVcas(double *Mach,double vcas) {
 486
 487   bool result=false;
 488   double guess=1.5;
 489   xlo=xhi=0;
 490   xmin=0;xmax=50;
 491   sfunc=&FGInitialCondition::calcVcas;
 492   if(findInterval(vcas,guess)) {
 493     if(solve(&mach,vcas))
 494       result=true;
 495   }
 496   return result;
 497 }
 498
 499 //******************************************************************************
 500
 501 bool FGInitialCondition::getAlpha(void) {
 502   bool result=false;
 503   double guess=theta-gamma;
 504   xlo=xhi=0;
 505   xmin=fdmex->GetAircraft()->GetAlphaCLMin();
 506   xmax=fdmex->GetAircraft()->GetAlphaCLMax();
 507   sfunc=&FGInitialCondition::GammaEqOfAlpha;
 508   if(findInterval(0,guess)){
 509     if(solve(&alpha,0)){
 510       result=true;
 511       salpha=sin(alpha);
 512       calpha=cos(alpha);
 513     }
 514   }
 515   calcWindUVW();
 516   return result;
 517 }
 518
 519 //******************************************************************************
 520
 521 bool FGInitialCondition::getTheta(void) {
 522   bool result=false;
 523   double guess=alpha+gamma;
 524   xlo=xhi=0;
 525   xmin=-89;xmax=89;
 526   sfunc=&FGInitialCondition::GammaEqOfTheta;
 527   if(findInterval(0,guess)){
 528     if(solve(&theta,0)){
 529       result=true;
 530       stheta=sin(theta);
 531       ctheta=cos(theta);
 532     }
 533   }
 534   calcWindUVW();
 535   return result;
 536 }
 537
 538 //******************************************************************************
 539
 540 double FGInitialCondition::GammaEqOfTheta(double Theta) {
 541   double a,b,c;
 542   double sTheta,cTheta;
 543
 544   //theta=Theta; stheta=sin(theta); ctheta=cos(theta);
 545   sTheta=sin(Theta); cTheta=cos(Theta);
 546   calcWindUVW();
 547   a=wdown + vt*calpha*cbeta + uw;
 548   b=vt*sphi*sbeta + vw*sphi;
 549   c=vt*cphi*salpha*cbeta + ww*cphi;
 550   return vt*sgamma - ( a*sTheta - (b+c)*cTheta);
 551 }
 552
 553 //******************************************************************************
 554
 555 double FGInitialCondition::GammaEqOfAlpha(double Alpha) {
 556   double a,b,c;
 557   double sAlpha,cAlpha;
 558
 559   sAlpha=sin(Alpha); cAlpha=cos(Alpha);
 560   a=wdown + vt*cAlpha*cbeta + uw;
 561   b=vt*sphi*sbeta + vw*sphi;
 562   c=vt*cphi*sAlpha*cbeta + ww*cphi;
 563
 564   return vt*sgamma - ( a*stheta - (b+c)*ctheta );
 565 }
 566
 567 //******************************************************************************
 568
 569 double FGInitialCondition::calcVcas(double Mach) {
 570
 571   double p=fdmex->GetAtmosphere()->GetPressure();
 572   double psl=fdmex->GetAtmosphere()->GetPressureSL();
 573   double rhosl=fdmex->GetAtmosphere()->GetDensitySL();
 574   double pt,A,B,D,vcas;
 575   if(Mach < 0) Mach=0;
 576   if(Mach < 1)    //calculate total pressure assuming isentropic flow
 577     pt=p*pow((1 + 0.2*Mach*Mach),3.5);
 578   else {
 579     // shock in front of pitot tube, we'll assume its normal and use
 580     // the Rayleigh Pitot Tube Formula, i.e. the ratio of total
 581     // pressure behind the shock to the static pressure in front
 582
 583
 584     //the normal shock assumption should not be a bad one -- most supersonic
 585     //aircraft place the pitot probe out front so that it is the forward
 586     //most point on the aircraft.  The real shock would, of course, take
 587     //on something like the shape of a rounded-off cone but, here again,
 588     //the assumption should be good since the opening of the pitot probe
 589     //is very small and, therefore, the effects of the shock curvature
 590     //should be small as well. AFAIK, this approach is fairly well accepted
 591     //within the aerospace community
 592
 593     B = 5.76*Mach*Mach/(5.6*Mach*Mach - 0.8);
 594
 595     // The denominator above is zero for Mach ~ 0.38, for which
 596     // we'll never be here, so we're safe
 597
 598     D = (2.8*Mach*Mach-0.4)*0.4167;
 599     pt = p*pow(B,3.5)*D;
 600   }
 601
 602   A = pow(((pt-p)/psl+1),0.28571);
 603   vcas = sqrt(7*psl/rhosl*(A-1));
 604   //cout << "calcVcas: vcas= " << vcas*fpstokts << " mach= " << Mach << " pressure: " << pt << endl;
 605   return vcas;
 606 }
 607
 608 //******************************************************************************
 609
 610 bool FGInitialCondition::findInterval(double x,double guess) {
 611   //void find_interval(inter_params &ip,eqfunc f,double y,double constant, int &flag){
 612
 613   int i=0;
 614   bool found=false;
 615   double flo,fhi,fguess;
 616   double lo,hi,step;
 617   step=0.1;
 618   fguess=(this->*sfunc)(guess)-x;
 619   lo=hi=guess;
 620   do {
 621     step=2*step;
 622     lo-=step;
 623     hi+=step;
 624     if(lo < xmin) lo=xmin;
 625     if(hi > xmax) hi=xmax;
 626     i++;
 627     flo=(this->*sfunc)(lo)-x;
 628     fhi=(this->*sfunc)(hi)-x;
 629     if(flo*fhi <=0) {  //found interval with root
 630       found=true;
 631       if(flo*fguess <= 0) {  //narrow interval down a bit
 632         hi=lo+step;    //to pass solver interval that is as
 633         //small as possible
 634       }
 635       else if(fhi*fguess <= 0) {
 636         lo=hi-step;
 637       }
 638     }
 639     //cout << "findInterval: i=" << i << " Lo= " << lo << " Hi= " << hi << endl;
 640   }
 641   while((found == 0) && (i <= 100));
 642   xlo=lo;
 643   xhi=hi;
 644   return found;
 645 }
 646
 647 //******************************************************************************
 648
 649 bool FGInitialCondition::solve(double *y,double x)
 650 {
 651   double x1,x2,x3,f1,f2,f3,d,d0;
 652   double eps=1E-5;
 653   double const relax =0.9;
 654   int i;
 655   bool success=false;
 656
 657   //initializations
 658   d=1;
 659   x2 = 0;
 660   x1=xlo;x3=xhi;
 661   f1=(this->*sfunc)(x1)-x;
 662   f3=(this->*sfunc)(x3)-x;
 663   d0=fabs(x3-x1);
 664
 665   //iterations
 666   i=0;
 667   while ((fabs(d) > eps) && (i < 100)) {
 668     d=(x3-x1)/d0;
 669     x2 = x1-d*d0*f1/(f3-f1);
 670
 671     f2=(this->*sfunc)(x2)-x;
 672     //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
 673     //cout << "                " << f1 << "," << f2 << "," << f3 << endl;
 674
 675     if(fabs(f2) <= 0.001) {
 676       x1=x3=x2;
 677     } else if(f1*f2 <= 0.0) {
 678       x3=x2;
 679       f3=f2;
 680       f1=relax*f1;
 681     } else if(f2*f3 <= 0) {
 682       x1=x2;
 683       f1=f2;
 684       f3=relax*f3;
 685     }
 686     //cout << i << endl;
 687     i++;
 688   }//end while
 689   if(i < 100) {
 690     success=true;
 691     *y=x2;
 692   }
 693
 694   //cout << "Success= " << success << " Vcas: " << vcas*fpstokts << " Mach: " << x2 << endl;
 695   return success;
 696 }
 697
 698 //******************************************************************************
 699
 700 double FGInitialCondition::GetWindDirDegIC(void) {
 701   if(weast != 0.0)
 702     return atan2(weast,wnorth)*radtodeg;
 703   else if(wnorth > 0)
 704     return 0.0;
 705   else
 706     return 180.0;
 707 }
 708
 709 //******************************************************************************
 710
 711 bool FGInitialCondition::Load(string acpath, string acname, string rstfile)
 712 {
 713   string resetDef;
 714   string token="";
 715
 716   double temp;
 717
 718 # ifndef macintosh
 719   resetDef = acpath + "/" + acname + "/" + rstfile + ".xml";
 720 # else
 721   resetDef = acpath + ";" + acname + ";" + rstfile + ".xml";
 722 # endif
 723
 724   FGConfigFile resetfile(resetDef);
 725   if (!resetfile.IsOpen()) {
 726     cerr << "Failed to open reset file: " << resetDef << endl;
 727     return false;
 728   }
 729
 730   resetfile.GetNextConfigLine();
 731   token = resetfile.GetValue();
 732   if (token != string("initialize")) {
 733     cerr << "The reset file " << resetDef
 734          << " does not appear to be a reset file" << endl;
 735     return false;
 736   }
 737
 738   resetfile.GetNextConfigLine();
 739   resetfile >> token;
 740   while (token != string("/initialize") && token != string("EOF")) {
 741     if (token == "UBODY" ) { resetfile >> temp; SetUBodyFpsIC(temp); }
 742     if (token == "VBODY" ) { resetfile >> temp; SetVBodyFpsIC(temp); }
 743     if (token == "WBODY" ) { resetfile >> temp; SetWBodyFpsIC(temp); }
 744     if (token == "LATITUDE" ) { resetfile >> temp; SetLatitudeDegIC(temp); }
 745     if (token == "LONGITUDE" ) { resetfile >> temp; SetLongitudeDegIC(temp); }
 746     if (token == "PHI" ) { resetfile >> temp; SetRollAngleDegIC(temp); }
 747     if (token == "THETA" ) { resetfile >> temp; SetPitchAngleDegIC(temp); }
 748     if (token == "PSI" ) { resetfile >> temp; SetTrueHeadingDegIC(temp); }
 749     if (token == "ALPHA" ) { resetfile >> temp; SetAlphaDegIC(temp); }
 750     if (token == "BETA" ) { resetfile >> temp; SetBetaDegIC(temp); }
 751     if (token == "GAMMA" ) { resetfile >> temp; SetFlightPathAngleDegIC(temp); }
 752     if (token == "ROC" ) { resetfile >> temp; SetClimbRateFpmIC(temp); }
 753     if (token == "ALTITUDE" ) { resetfile >> temp; SetAltitudeFtIC(temp); }
 754     if (token == "WINDDIR" ) { resetfile >> temp; SetWindDirDegIC(temp); }
 755     if (token == "VWIND" ) { resetfile >> temp; SetWindMagKtsIC(temp); }
 756     if (token == "HWIND" ) { resetfile >> temp; SetHeadWindKtsIC(temp); }
 757     if (token == "XWIND" ) { resetfile >> temp; SetCrossWindKtsIC(temp); }
 758     if (token == "VC" ) { resetfile >> temp; SetVcalibratedKtsIC(temp); }
 759     if (token == "MACH" ) { resetfile >> temp; SetMachIC(temp); }
 760     if (token == "VGROUND" ) { resetfile >> temp; SetVgroundKtsIC(temp); }
 761     resetfile >> token;
 762   }
 763
 764   fdmex->RunIC(this);
 765
 766   return true;
 767 }
 768
 769 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 770 //    The bitmasked value choices are as follows:
 771 //    unset: In this case (the default) JSBSim would only print
 772 //       out the normally expected messages, essentially echoing
 773 //       the config files as they are read. If the environment
 774 //       variable is not set, debug_lvl is set to 1 internally
 775 //    0: This requests JSBSim not to output any messages
 776 //       whatsoever.
 777 //    1: This value explicity requests the normal JSBSim
 778 //       startup messages
 779 //    2: This value asks for a message to be printed out when
 780 //       a class is instantiated
 781 //    4: When this value is set, a message is displayed when a
 782 //       FGModel object executes its Run() method
 783 //    8: When this value is set, various runtime state variables
 784 //       are printed out periodically
 785 //    16: When set various parameters are sanity checked and
 786 //       a message is printed out when they go out of bounds
 787
 788 void FGInitialCondition::Debug(int from)
 789 {
 790   if (debug_lvl <= 0) return;
 791
 792   if (debug_lvl & 1) { // Standard console startup message output
 793   }
 794   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 795     if (from == 0) cout << "Instantiated: FGInitialCondition" << endl;
 796     if (from == 1) cout << "Destroyed:    FGInitialCondition" << endl;
 797   }
 798   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 799   }
 800   if (debug_lvl & 8 ) { // Runtime state variables
 801   }
 802   if (debug_lvl & 16) { // Sanity checking
 803   }
 804   if (debug_lvl & 64) {
 805     if (from == 0) { // Constructor
 806       cout << IdSrc << endl;
 807       cout << IdHdr << endl;
 808     }
 809   }
 810 }
 811