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   if (debug_lvl & 2) cout << "Instantiated: FGInitialCondition" << endl;
  95 }
  96
  97 //******************************************************************************
  98
  99 FGInitialCondition::~FGInitialCondition()
 100 {
 101   if (debug_lvl & 2) cout << "Destroyed:    FGInitialCondition" << endl;
 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     }
 370     uw=wnorth*ctheta*cpsi +
 371        weast*ctheta*spsi -
 372        wdown*stheta;
 373     vw=wnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 374         weast*( sphi*stheta*spsi + cphi*cpsi ) +
 375        wdown*sphi*ctheta;
 376     ww=wnorth*(cphi*stheta*cpsi + sphi*spsi) +
 377        weast*(cphi*stheta*spsi - sphi*cpsi) +
 378        wdown*cphi*ctheta;
 379
 380
 381     /* cout << "FGInitialCondition::calcWindUVW: wnorth, weast, wdown "
 382          << wnorth << ", " << weast << ", " << wdown << endl;
 383     cout << "FGInitialCondition::calcWindUVW: theta, phi, psi "
 384           << theta << ", " << phi << ", " << psi << endl;
 385     cout << "FGInitialCondition::calcWindUVW: uw, vw, ww "
 386           << uw << ", " << vw << ", " << ww << endl;   */
 387
 388 }
 389
 390 //******************************************************************************
 391
 392 void FGInitialCondition::SetAltitudeFtIC(double tt) {
 393   altitude=tt;
 394   fdmex->GetPosition()->Seth(altitude);
 395   fdmex->GetAtmosphere()->Run();
 396   //lets try to make sure the user gets what they intended
 397
 398   switch(lastSpeedSet) {
 399   case setned:
 400   case setuvw:
 401   case setvt:
 402     SetVtrueKtsIC(vt*fpstokts);
 403     break;
 404   case setvc:
 405     SetVcalibratedKtsIC(vc*fpstokts);
 406     break;
 407   case setve:
 408     SetVequivalentKtsIC(ve*fpstokts);
 409     break;
 410   case setmach:
 411     SetMachIC(mach);
 412     break;
 413   case setvg:
 414     SetVgroundFpsIC(vg);
 415     break;
 416   }
 417 }
 418
 419 //******************************************************************************
 420
 421 void FGInitialCondition::SetAltitudeAGLFtIC(double tt) {
 422   fdmex->GetPosition()->SetDistanceAGL(tt);
 423   altitude=fdmex->GetPosition()->Geth();
 424   SetAltitudeFtIC(altitude);
 425 }
 426
 427 //******************************************************************************
 428
 429 void FGInitialCondition::SetSeaLevelRadiusFtIC(double tt) {
 430   sea_level_radius = tt;
 431 }
 432
 433 //******************************************************************************
 434
 435 void FGInitialCondition::SetTerrainAltitudeFtIC(double tt) {
 436   terrain_altitude=tt;
 437 }
 438
 439 //******************************************************************************
 440
 441 void FGInitialCondition::calcUVWfromNED(void) {
 442   u=vnorth*ctheta*cpsi +
 443      veast*ctheta*spsi -
 444      vdown*stheta;
 445   v=vnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 446      veast*( sphi*stheta*spsi + cphi*cpsi ) +
 447      vdown*sphi*ctheta;
 448   w=vnorth*( cphi*stheta*cpsi + sphi*spsi ) +
 449      veast*( cphi*stheta*spsi - sphi*cpsi ) +
 450      vdown*cphi*ctheta;
 451 }
 452
 453 //******************************************************************************
 454
 455 void FGInitialCondition::SetVnorthFpsIC(double tt) {
 456   vnorth=tt;
 457   calcUVWfromNED();
 458   vt=sqrt(u*u + v*v + w*w);
 459   lastSpeedSet=setned;
 460 }
 461
 462 //******************************************************************************
 463
 464 void FGInitialCondition::SetVeastFpsIC(double tt) {
 465   veast=tt;
 466   calcUVWfromNED();
 467   vt=sqrt(u*u + v*v + w*w);
 468   lastSpeedSet=setned;
 469 }
 470
 471 //******************************************************************************
 472
 473 void FGInitialCondition::SetVdownFpsIC(double tt) {
 474   vdown=tt;
 475   calcUVWfromNED();
 476   vt=sqrt(u*u + v*v + w*w);
 477   SetClimbRateFpsIC(-1*vdown);
 478   lastSpeedSet=setned;
 479 }
 480
 481 //******************************************************************************
 482
 483 bool FGInitialCondition::getMachFromVcas(double *Mach,double vcas) {
 484
 485   bool result=false;
 486   double guess=1.5;
 487   xlo=xhi=0;
 488   xmin=0;xmax=50;
 489   sfunc=&FGInitialCondition::calcVcas;
 490   if(findInterval(vcas,guess)) {
 491     if(solve(&mach,vcas))
 492       result=true;
 493   }
 494   return result;
 495 }
 496
 497 //******************************************************************************
 498
 499 bool FGInitialCondition::getAlpha(void) {
 500   bool result=false;
 501   double guess=theta-gamma;
 502   xlo=xhi=0;
 503   xmin=fdmex->GetAircraft()->GetAlphaCLMin();
 504   xmax=fdmex->GetAircraft()->GetAlphaCLMax();
 505   sfunc=&FGInitialCondition::GammaEqOfAlpha;
 506   if(findInterval(0,guess)){
 507     if(solve(&alpha,0)){
 508       result=true;
 509       salpha=sin(alpha);
 510       calpha=cos(alpha);
 511     }
 512   }
 513   calcWindUVW();
 514   return result;
 515 }
 516
 517 //******************************************************************************
 518
 519 bool FGInitialCondition::getTheta(void) {
 520   bool result=false;
 521   double guess=alpha+gamma;
 522   xlo=xhi=0;
 523   xmin=-89;xmax=89;
 524   sfunc=&FGInitialCondition::GammaEqOfTheta;
 525   if(findInterval(0,guess)){
 526     if(solve(&theta,0)){
 527       result=true;
 528       stheta=sin(theta);
 529       ctheta=cos(theta);
 530     }
 531   }
 532   calcWindUVW();
 533   return result;
 534 }
 535
 536 //******************************************************************************
 537
 538 double FGInitialCondition::GammaEqOfTheta(double Theta) {
 539   double a,b,c,d;
 540   double sTheta,cTheta;
 541
 542   //theta=Theta; stheta=sin(theta); ctheta=cos(theta);
 543   sTheta=sin(Theta); cTheta=cos(Theta);
 544   calcWindUVW();
 545   a=wdown + vt*calpha*cbeta + uw;
 546   b=vt*sphi*sbeta + vw*sphi;
 547   c=vt*cphi*salpha*cbeta + ww*cphi;
 548   return vt*sgamma - ( a*sTheta - (b+c)*cTheta);
 549 }
 550
 551 //******************************************************************************
 552
 553 double FGInitialCondition::GammaEqOfAlpha(double Alpha) {
 554   double a,b,c,d;
 555   double sAlpha,cAlpha;
 556
 557   sAlpha=sin(Alpha); cAlpha=cos(Alpha);
 558   a=wdown + vt*cAlpha*cbeta + uw;
 559   b=vt*sphi*sbeta + vw*sphi;
 560   c=vt*cphi*sAlpha*cbeta + ww*cphi;
 561
 562   return vt*sgamma - ( a*stheta - (b+c)*ctheta );
 563 }
 564
 565 //******************************************************************************
 566
 567 double FGInitialCondition::calcVcas(double Mach) {
 568
 569   double p=fdmex->GetAtmosphere()->GetPressure();
 570   double psl=fdmex->GetAtmosphere()->GetPressureSL();
 571   double rhosl=fdmex->GetAtmosphere()->GetDensitySL();
 572   double pt,A,B,D,vcas;
 573   if(Mach < 0) Mach=0;
 574   if(Mach < 1)    //calculate total pressure assuming isentropic flow
 575     pt=p*pow((1 + 0.2*Mach*Mach),3.5);
 576   else {
 577     // shock in front of pitot tube, we'll assume its normal and use
 578     // the Rayleigh Pitot Tube Formula, i.e. the ratio of total
 579     // pressure behind the shock to the static pressure in front
 580
 581
 582     //the normal shock assumption should not be a bad one -- most supersonic
 583     //aircraft place the pitot probe out front so that it is the forward
 584     //most point on the aircraft.  The real shock would, of course, take
 585     //on something like the shape of a rounded-off cone but, here again,
 586     //the assumption should be good since the opening of the pitot probe
 587     //is very small and, therefore, the effects of the shock curvature
 588     //should be small as well. AFAIK, this approach is fairly well accepted
 589     //within the aerospace community
 590
 591     B = 5.76*Mach*Mach/(5.6*Mach*Mach - 0.8);
 592
 593     // The denominator above is zero for Mach ~ 0.38, for which
 594     // we'll never be here, so we're safe
 595
 596     D = (2.8*Mach*Mach-0.4)*0.4167;
 597     pt = p*pow(B,3.5)*D;
 598   }
 599
 600   A = pow(((pt-p)/psl+1),0.28571);
 601   vcas = sqrt(7*psl/rhosl*(A-1));
 602   //cout << "calcVcas: vcas= " << vcas*fpstokts << " mach= " << Mach << " pressure: " << pt << endl;
 603   return vcas;
 604 }
 605
 606 //******************************************************************************
 607
 608 bool FGInitialCondition::findInterval(double x,double guess) {
 609   //void find_interval(inter_params &ip,eqfunc f,double y,double constant, int &flag){
 610
 611   int i=0;
 612   bool found=false;
 613   double flo,fhi,fguess;
 614   double lo,hi,step;
 615   step=0.1;
 616   fguess=(this->*sfunc)(guess)-x;
 617   lo=hi=guess;
 618   do {
 619     step=2*step;
 620     lo-=step;
 621     hi+=step;
 622     if(lo < xmin) lo=xmin;
 623     if(hi > xmax) hi=xmax;
 624     i++;
 625     flo=(this->*sfunc)(lo)-x;
 626     fhi=(this->*sfunc)(hi)-x;
 627     if(flo*fhi <=0) {  //found interval with root
 628       found=true;
 629       if(flo*fguess <= 0) {  //narrow interval down a bit
 630         hi=lo+step;    //to pass solver interval that is as
 631         //small as possible
 632       }
 633       else if(fhi*fguess <= 0) {
 634         lo=hi-step;
 635       }
 636     }
 637     //cout << "findInterval: i=" << i << " Lo= " << lo << " Hi= " << hi << endl;
 638   }
 639   while((found == 0) && (i <= 100));
 640   xlo=lo;
 641   xhi=hi;
 642   return found;
 643 }
 644
 645 //******************************************************************************
 646
 647 bool FGInitialCondition::solve(double *y,double x)
 648 {
 649   double x1,x2,x3,f1,f2,f3,d,d0;
 650   double eps=1E-5;
 651   double const relax =0.9;
 652   int i;
 653   bool success=false;
 654
 655   //initializations
 656   d=1;
 657
 658   x1=xlo;x3=xhi;
 659   f1=(this->*sfunc)(x1)-x;
 660   f3=(this->*sfunc)(x3)-x;
 661   d0=fabs(x3-x1);
 662
 663   //iterations
 664   i=0;
 665   while ((fabs(d) > eps) && (i < 100)) {
 666     d=(x3-x1)/d0;
 667     if (f3-f1 != 0.0) x2 = x1-d*d0*f1/(f3-f1);
 668     else              x2 = x1-d*d0*f1/0.01; // TONY: What is correct? This is a WAG.
 669
 670     f2=(this->*sfunc)(x2)-x;
 671     //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
 672     //cout << "                " << f1 << "," << f2 << "," << f3 << endl;
 673
 674     if(fabs(f2) <= 0.001) {
 675       x1=x3=x2;
 676     } else if(f1*f2 <= 0.0) {
 677       x3=x2;
 678       f3=f2;
 679       f1=relax*f1;
 680     } else if(f2*f3 <= 0) {
 681       x1=x2;
 682       f1=f2;
 683       f3=relax*f3;
 684     }
 685     //cout << i << endl;
 686     i++;
 687   }//end while
 688   if(i < 100) {
 689     success=true;
 690     *y=x2;
 691   }
 692
 693   //cout << "Success= " << success << " Vcas: " << vcas*fpstokts << " Mach: " << x2 << endl;
 694   return success;
 695 }
 696
 697 //******************************************************************************
 698
 699 double FGInitialCondition::GetWindDirDegIC(void) {
 700   if(weast != 0.0)
 701     return atan2(weast,wnorth)*radtodeg;
 702   else if(wnorth > 0)
 703     return 0.0;
 704   else
 705     return 180.0;
 706 }
 707
 708 //******************************************************************************
 709
 710 bool FGInitialCondition::Load(string acpath, string acname, string rstfile)
 711 {
 712   string resetDef;
 713   string token="";
 714
 715   double temp;
 716
 717 # ifndef macintosh
 718   resetDef = acpath + "/" + acname + "/" + rstfile + ".xml";
 719 # else
 720   resetDef = acpath + ";" + acname + ";" + rstfile + ".xml";
 721 # endif
 722
 723   FGConfigFile resetfile(resetDef);
 724   if (!resetfile.IsOpen()) return false;
 725
 726   resetfile.GetNextConfigLine();
 727   token = resetfile.GetValue();
 728   if (token != string("initialize")) {
 729     cerr << "The reset file " << resetDef
 730          << " does not appear to be a reset file" << endl;
 731     return false;
 732   }
 733
 734   resetfile.GetNextConfigLine();
 735   resetfile >> token;
 736   while (token != string("/initialize") && token != string("EOF")) {
 737     if (token == "UBODY" ) { resetfile >> temp; SetUBodyFpsIC(temp); }
 738     if (token == "VBODY" ) { resetfile >> temp; SetVBodyFpsIC(temp); }
 739     if (token == "WBODY" ) { resetfile >> temp; SetWBodyFpsIC(temp); }
 740     if (token == "LATITUDE" ) { resetfile >> temp; SetLatitudeDegIC(temp); }
 741     if (token == "LONGITUDE" ) { resetfile >> temp; SetLongitudeDegIC(temp); }
 742     if (token == "PHI" ) { resetfile >> temp; SetRollAngleDegIC(temp); }
 743     if (token == "THETA" ) { resetfile >> temp; SetPitchAngleDegIC(temp); }
 744     if (token == "PSI" ) { resetfile >> temp; SetTrueHeadingDegIC(temp); }
 745     if (token == "ALPHA" ) { resetfile >> temp; SetAlphaDegIC(temp); }
 746     if (token == "BETA" ) { resetfile >> temp; SetBetaDegIC(temp); }
 747     if (token == "GAMMA" ) { resetfile >> temp; SetFlightPathAngleDegIC(temp); }
 748     if (token == "ROC" ) { resetfile >> temp; SetClimbRateFpmIC(temp); }
 749     if (token == "ALTITUDE" ) { resetfile >> temp; SetAltitudeFtIC(temp); }
 750     if (token == "WINDDIR" ) { resetfile >> temp; SetWindDirDegIC(temp); }
 751     if (token == "VWIND" ) { resetfile >> temp; SetWindMagKtsIC(temp); }
 752     if (token == "HWIND" ) { resetfile >> temp; SetHeadWindKtsIC(temp); }
 753     if (token == "XWIND" ) { resetfile >> temp; SetCrossWindKtsIC(temp); }
 754     if (token == "VC" ) { resetfile >> temp; SetVcalibratedKtsIC(temp); }
 755     if (token == "MACH" ) { resetfile >> temp; SetMachIC(temp); }
 756     if (token == "VGROUND" ) { resetfile >> temp; SetVgroundKtsIC(temp); }
 757     resetfile >> token;
 758   }
 759
 760   fdmex->RunIC(this);
 761
 762   return true;
 763 }