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
 505   if(vt < 0.01) return 0;
 506
 507   xlo=xhi=0;
 508   xmin=fdmex->GetAircraft()->GetAlphaCLMin();
 509   xmax=fdmex->GetAircraft()->GetAlphaCLMax();
 510   sfunc=&FGInitialCondition::GammaEqOfAlpha;
 511   if(findInterval(0,guess)){
 512     if(solve(&alpha,0)){
 513       result=true;
 514       salpha=sin(alpha);
 515       calpha=cos(alpha);
 516     }
 517   }
 518   calcWindUVW();
 519   return result;
 520 }
 521
 522 //******************************************************************************
 523
 524 bool FGInitialCondition::getTheta(void) {
 525   bool result=false;
 526   double guess=alpha+gamma;
 527
 528   if(vt < 0.01) return 0;
 529
 530   xlo=xhi=0;
 531   xmin=-89;xmax=89;
 532   sfunc=&FGInitialCondition::GammaEqOfTheta;
 533   if(findInterval(0,guess)){
 534     if(solve(&theta,0)){
 535       result=true;
 536       stheta=sin(theta);
 537       ctheta=cos(theta);
 538     }
 539   }
 540   calcWindUVW();
 541   return result;
 542 }
 543
 544 //******************************************************************************
 545
 546 double FGInitialCondition::GammaEqOfTheta(double Theta) {
 547   double a,b,c;
 548   double sTheta,cTheta;
 549
 550   //theta=Theta; stheta=sin(theta); ctheta=cos(theta);
 551   sTheta=sin(Theta); cTheta=cos(Theta);
 552   calcWindUVW();
 553   a=wdown + vt*calpha*cbeta + uw;
 554   b=vt*sphi*sbeta + vw*sphi;
 555   c=vt*cphi*salpha*cbeta + ww*cphi;
 556   return vt*sgamma - ( a*sTheta - (b+c)*cTheta);
 557 }
 558
 559 //******************************************************************************
 560
 561 double FGInitialCondition::GammaEqOfAlpha(double Alpha) {
 562   double a,b,c;
 563   double sAlpha,cAlpha;
 564   sAlpha=sin(Alpha); cAlpha=cos(Alpha);
 565   a=wdown + vt*cAlpha*cbeta + uw;
 566   b=vt*sphi*sbeta + vw*sphi;
 567   c=vt*cphi*sAlpha*cbeta + ww*cphi;
 568
 569   return vt*sgamma - ( a*stheta - (b+c)*ctheta );
 570 }
 571
 572 //******************************************************************************
 573
 574 double FGInitialCondition::calcVcas(double Mach) {
 575
 576   double p=fdmex->GetAtmosphere()->GetPressure();
 577   double psl=fdmex->GetAtmosphere()->GetPressureSL();
 578   double rhosl=fdmex->GetAtmosphere()->GetDensitySL();
 579   double pt,A,B,D,vcas;
 580   if(Mach < 0) Mach=0;
 581   if(Mach < 1)    //calculate total pressure assuming isentropic flow
 582     pt=p*pow((1 + 0.2*Mach*Mach),3.5);
 583   else {
 584     // shock in front of pitot tube, we'll assume its normal and use
 585     // the Rayleigh Pitot Tube Formula, i.e. the ratio of total
 586     // pressure behind the shock to the static pressure in front
 587
 588
 589     //the normal shock assumption should not be a bad one -- most supersonic
 590     //aircraft place the pitot probe out front so that it is the forward
 591     //most point on the aircraft.  The real shock would, of course, take
 592     //on something like the shape of a rounded-off cone but, here again,
 593     //the assumption should be good since the opening of the pitot probe
 594     //is very small and, therefore, the effects of the shock curvature
 595     //should be small as well. AFAIK, this approach is fairly well accepted
 596     //within the aerospace community
 597
 598     B = 5.76*Mach*Mach/(5.6*Mach*Mach - 0.8);
 599
 600     // The denominator above is zero for Mach ~ 0.38, for which
 601     // we'll never be here, so we're safe
 602
 603     D = (2.8*Mach*Mach-0.4)*0.4167;
 604     pt = p*pow(B,3.5)*D;
 605   }
 606
 607   A = pow(((pt-p)/psl+1),0.28571);
 608   vcas = sqrt(7*psl/rhosl*(A-1));
 609   //cout << "calcVcas: vcas= " << vcas*fpstokts << " mach= " << Mach << " pressure: " << pt << endl;
 610   return vcas;
 611 }
 612
 613 //******************************************************************************
 614
 615 bool FGInitialCondition::findInterval(double x,double guess) {
 616   //void find_interval(inter_params &ip,eqfunc f,double y,double constant, int &flag){
 617
 618   int i=0;
 619   bool found=false;
 620   double flo,fhi,fguess;
 621   double lo,hi,step;
 622   step=0.1;
 623   fguess=(this->*sfunc)(guess)-x;
 624   lo=hi=guess;
 625   do {
 626     step=2*step;
 627     lo-=step;
 628     hi+=step;
 629     if(lo < xmin) lo=xmin;
 630     if(hi > xmax) hi=xmax;
 631     i++;
 632     flo=(this->*sfunc)(lo)-x;
 633     fhi=(this->*sfunc)(hi)-x;
 634     if(flo*fhi <=0) {  //found interval with root
 635       found=true;
 636       if(flo*fguess <= 0) {  //narrow interval down a bit
 637         hi=lo+step;    //to pass solver interval that is as
 638         //small as possible
 639       }
 640       else if(fhi*fguess <= 0) {
 641         lo=hi-step;
 642       }
 643     }
 644     //cout << "findInterval: i=" << i << " Lo= " << lo << " Hi= " << hi << endl;
 645   }
 646   while((found == 0) && (i <= 100));
 647   xlo=lo;
 648   xhi=hi;
 649   return found;
 650 }
 651
 652 //******************************************************************************
 653
 654 bool FGInitialCondition::solve(double *y,double x)
 655 {
 656   double x1,x2,x3,f1,f2,f3,d,d0;
 657   double eps=1E-5;
 658   double const relax =0.9;
 659   int i;
 660   bool success=false;
 661
 662   //initializations
 663   d=1;
 664   x2 = 0;
 665   x1=xlo;x3=xhi;
 666   f1=(this->*sfunc)(x1)-x;
 667   f3=(this->*sfunc)(x3)-x;
 668   d0=fabs(x3-x1);
 669
 670   //iterations
 671   i=0;
 672   while ((fabs(d) > eps) && (i < 100)) {
 673     d=(x3-x1)/d0;
 674     x2 = x1-d*d0*f1/(f3-f1);
 675
 676     f2=(this->*sfunc)(x2)-x;
 677     //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
 678     //cout << "                " << f1 << "," << f2 << "," << f3 << endl;
 679
 680     if(fabs(f2) <= 0.001) {
 681       x1=x3=x2;
 682     } else if(f1*f2 <= 0.0) {
 683       x3=x2;
 684       f3=f2;
 685       f1=relax*f1;
 686     } else if(f2*f3 <= 0) {
 687       x1=x2;
 688       f1=f2;
 689       f3=relax*f3;
 690     }
 691     //cout << i << endl;
 692     i++;
 693   }//end while
 694   if(i < 100) {
 695     success=true;
 696     *y=x2;
 697   }
 698
 699   //cout << "Success= " << success << " Vcas: " << vcas*fpstokts << " Mach: " << x2 << endl;
 700   return success;
 701 }
 702
 703 //******************************************************************************
 704
 705 double FGInitialCondition::GetWindDirDegIC(void) {
 706   if(weast != 0.0)
 707     return atan2(weast,wnorth)*radtodeg;
 708   else if(wnorth > 0)
 709     return 0.0;
 710   else
 711     return 180.0;
 712 }
 713
 714 //******************************************************************************
 715
 716 bool FGInitialCondition::Load(string acpath, string acname, string rstfile)
 717 {
 718   string resetDef;
 719   string token="";
 720
 721   double temp;
 722
 723 # ifndef macintosh
 724   resetDef = acpath + "/" + acname + "/" + rstfile + ".xml";
 725 # else
 726   resetDef = acpath + ";" + acname + ";" + rstfile + ".xml";
 727 # endif
 728
 729   FGConfigFile resetfile(resetDef);
 730   if (!resetfile.IsOpen()) {
 731     cerr << "Failed to open reset file: " << resetDef << endl;
 732     return false;
 733   }
 734
 735   resetfile.GetNextConfigLine();
 736   token = resetfile.GetValue();
 737   if (token != string("initialize")) {
 738     cerr << "The reset file " << resetDef
 739          << " does not appear to be a reset file" << endl;
 740     return false;
 741   }
 742
 743   resetfile.GetNextConfigLine();
 744   resetfile >> token;
 745   while (token != string("/initialize") && token != string("EOF")) {
 746     if (token == "UBODY" ) { resetfile >> temp; SetUBodyFpsIC(temp); }
 747     if (token == "VBODY" ) { resetfile >> temp; SetVBodyFpsIC(temp); }
 748     if (token == "WBODY" ) { resetfile >> temp; SetWBodyFpsIC(temp); }
 749     if (token == "LATITUDE" ) { resetfile >> temp; SetLatitudeDegIC(temp); }
 750     if (token == "LONGITUDE" ) { resetfile >> temp; SetLongitudeDegIC(temp); }
 751     if (token == "PHI" ) { resetfile >> temp; SetRollAngleDegIC(temp); }
 752     if (token == "THETA" ) { resetfile >> temp; SetPitchAngleDegIC(temp); }
 753     if (token == "PSI" ) { resetfile >> temp; SetTrueHeadingDegIC(temp); }
 754     if (token == "ALPHA" ) { resetfile >> temp; SetAlphaDegIC(temp); }
 755     if (token == "BETA" ) { resetfile >> temp; SetBetaDegIC(temp); }
 756     if (token == "GAMMA" ) { resetfile >> temp; SetFlightPathAngleDegIC(temp); }
 757     if (token == "ROC" ) { resetfile >> temp; SetClimbRateFpmIC(temp); }
 758     if (token == "ALTITUDE" ) { resetfile >> temp; SetAltitudeFtIC(temp); }
 759     if (token == "WINDDIR" ) { resetfile >> temp; SetWindDirDegIC(temp); }
 760     if (token == "VWIND" ) { resetfile >> temp; SetWindMagKtsIC(temp); }
 761     if (token == "HWIND" ) { resetfile >> temp; SetHeadWindKtsIC(temp); }
 762     if (token == "XWIND" ) { resetfile >> temp; SetCrossWindKtsIC(temp); }
 763     if (token == "VC" ) { resetfile >> temp; SetVcalibratedKtsIC(temp); }
 764     if (token == "MACH" ) { resetfile >> temp; SetMachIC(temp); }
 765     if (token == "VGROUND" ) { resetfile >> temp; SetVgroundKtsIC(temp); }
 766     resetfile >> token;
 767   }
 768
 769   fdmex->RunIC(this);
 770
 771   return true;
 772 }
 773
 774 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 775 //    The bitmasked value choices are as follows:
 776 //    unset: In this case (the default) JSBSim would only print
 777 //       out the normally expected messages, essentially echoing
 778 //       the config files as they are read. If the environment
 779 //       variable is not set, debug_lvl is set to 1 internally
 780 //    0: This requests JSBSim not to output any messages
 781 //       whatsoever.
 782 //    1: This value explicity requests the normal JSBSim
 783 //       startup messages
 784 //    2: This value asks for a message to be printed out when
 785 //       a class is instantiated
 786 //    4: When this value is set, a message is displayed when a
 787 //       FGModel object executes its Run() method
 788 //    8: When this value is set, various runtime state variables
 789 //       are printed out periodically
 790 //    16: When set various parameters are sanity checked and
 791 //       a message is printed out when they go out of bounds
 792
 793 void FGInitialCondition::Debug(int from)
 794 {
 795   if (debug_lvl <= 0) return;
 796
 797   if (debug_lvl & 1) { // Standard console startup message output
 798   }
 799   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
 800     if (from == 0) cout << "Instantiated: FGInitialCondition" << endl;
 801     if (from == 1) cout << "Destroyed:    FGInitialCondition" << endl;
 802   }
 803   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
 804   }
 805   if (debug_lvl & 8 ) { // Runtime state variables
 806   }
 807   if (debug_lvl & 16) { // Sanity checking
 808   }
 809   if (debug_lvl & 64) {
 810     if (from == 0) { // Constructor
 811       cout << IdSrc << endl;
 812       cout << IdHdr << endl;
 813     }
 814   }
 815 }
 816