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 "FGAerodynamics.h"
  50 #include "FGFCS.h"
  51 #include "FGAircraft.h"
  52 #include "FGTranslation.h"
  53 #include "FGRotation.h"
  54 #include "FGPosition.h"
  55 #include "FGAuxiliary.h"
  56 #include "FGOutput.h"
  57 #include "FGConfigFile.h"
  58 #include "FGPropertyManager.h"
  59
  60 namespace JSBSim {
  61
  62 static const char *IdSrc = "$Id$";
  63 static const char *IdHdr = ID_INITIALCONDITION;
  64
  65 //******************************************************************************
  66
  67 FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec)
  68 {
  69   vt=vc=ve=vg=0;
  70   mach=0;
  71   alpha=beta=gamma=0;
  72   theta=phi=psi=0;
  73   altitude=hdot=0;
  74   latitude=longitude=0;
  75   u=v=w=0;
  76   uw=vw=ww=0;
  77   vnorth=veast=vdown=0;
  78   wnorth=weast=wdown=0;
  79   whead=wcross=0;
  80   wdir=wmag=0;
  81   lastSpeedSet=setvt;
  82   lastWindSet=setwned;
  83   sea_level_radius = FDMExec->GetInertial()->RefRadius();
  84   radius_to_vehicle = FDMExec->GetInertial()->RefRadius();
  85   terrain_altitude = 0;
  86
  87   salpha=sbeta=stheta=sphi=spsi=sgamma=0;
  88   calpha=cbeta=ctheta=cphi=cpsi=cgamma=1;
  89
  90   if(FDMExec != NULL ) {
  91     fdmex=FDMExec;
  92     fdmex->GetPosition()->Seth(altitude);
  93     fdmex->GetAtmosphere()->Run();
  94     PropertyManager=fdmex->GetPropertyManager();
  95     bind();
  96   } else {
  97     cout << "FGInitialCondition: This class requires a pointer to a valid FGFDMExec object" << endl;
  98   }
  99
 100   Debug(0);
 101 }
 102
 103 //******************************************************************************
 104
 105 FGInitialCondition::~FGInitialCondition()
 106 {
 107   unbind();
 108   Debug(1);
 109 }
 110
 111 //******************************************************************************
 112
 113 void FGInitialCondition::SetVcalibratedKtsIC(double tt) {
 114
 115   if(getMachFromVcas(&mach,tt*ktstofps)) {
 116     //cout << "Mach: " << mach << endl;
 117     lastSpeedSet=setvc;
 118     vc=tt*ktstofps;
 119     vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
 120     ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 121     //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
 122   }
 123   else {
 124     cout << "Failed to get Mach number for given Vc and altitude, Vc unchanged." << endl;
 125     cout << "Please mail the set of initial conditions used to apeden@earthlink.net" << endl;
 126   }
 127 }
 128
 129 //******************************************************************************
 130
 131 void FGInitialCondition::SetVequivalentKtsIC(double tt) {
 132   ve=tt*ktstofps;
 133   lastSpeedSet=setve;
 134   vt=ve*1/sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 135   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 136   vc=calcVcas(mach);
 137 }
 138
 139 //******************************************************************************
 140
 141 void FGInitialCondition::SetVgroundFpsIC(double tt) {
 142   double ua,va,wa;
 143   double vxz;
 144
 145   vg=tt;
 146   lastSpeedSet=setvg;
 147   vnorth = vg*cos(psi); veast = vg*sin(psi); vdown = 0;
 148   calcUVWfromNED();
 149   ua = u + uw; va = v + vw; wa = w + ww;
 150   vt = sqrt( ua*ua + va*va + wa*wa );
 151   alpha = beta = 0;
 152   vxz = sqrt( u*u + w*w );
 153   if( w != 0 ) alpha = atan2( w, u );
 154   if( vxz != 0 ) beta = atan2( v, vxz );
 155   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 156   vc=calcVcas(mach);
 157   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 158 }
 159
 160 //******************************************************************************
 161
 162 void FGInitialCondition::SetVtrueFpsIC(double tt) {
 163   vt=tt;
 164   lastSpeedSet=setvt;
 165   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
 166   vc=calcVcas(mach);
 167   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 168 }
 169
 170 //******************************************************************************
 171
 172 void FGInitialCondition::SetMachIC(double tt) {
 173   mach=tt;
 174   lastSpeedSet=setmach;
 175   vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
 176   vc=calcVcas(mach);
 177   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 178   //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
 179 }
 180
 181 //******************************************************************************
 182
 183 void FGInitialCondition::SetClimbRateFpmIC(double tt) {
 184   SetClimbRateFpsIC(tt/60.0);
 185 }
 186
 187 //******************************************************************************
 188
 189 void FGInitialCondition::SetClimbRateFpsIC(double tt) {
 190
 191   if(vt > 0.1) {
 192     hdot=tt;
 193     gamma=asin(hdot/vt);
 194     sgamma=sin(gamma); cgamma=cos(gamma);
 195   }
 196 }
 197
 198 //******************************************************************************
 199
 200 void FGInitialCondition::SetFlightPathAngleRadIC(double tt) {
 201   gamma=tt;
 202   sgamma=sin(gamma); cgamma=cos(gamma);
 203   getTheta();
 204   hdot=vt*sgamma;
 205 }
 206
 207 //******************************************************************************
 208
 209 void FGInitialCondition::SetAlphaRadIC(double tt) {
 210   alpha=tt;
 211   salpha=sin(alpha); calpha=cos(alpha);
 212   getTheta();
 213 }
 214
 215 //******************************************************************************
 216
 217 void FGInitialCondition::SetPitchAngleRadIC(double tt) {
 218   theta=tt;
 219   stheta=sin(theta); ctheta=cos(theta);
 220   getAlpha();
 221 }
 222
 223 //******************************************************************************
 224
 225 void FGInitialCondition::SetBetaRadIC(double tt) {
 226   beta=tt;
 227   sbeta=sin(beta); cbeta=cos(beta);
 228   getTheta();
 229
 230 }
 231
 232 //******************************************************************************
 233
 234 void FGInitialCondition::SetRollAngleRadIC(double tt) {
 235   phi=tt;
 236   sphi=sin(phi); cphi=cos(phi);
 237   getTheta();
 238 }
 239
 240 //******************************************************************************
 241
 242 void FGInitialCondition::SetTrueHeadingRadIC(double tt) {
 243     psi=tt;
 244     spsi=sin(psi); cpsi=cos(psi);
 245     calcWindUVW();
 246 }
 247
 248 //******************************************************************************
 249
 250 void FGInitialCondition::SetUBodyFpsIC(double tt) {
 251   u=tt;
 252   vt=sqrt(u*u + v*v + w*w);
 253   lastSpeedSet=setuvw;
 254 }
 255
 256 //******************************************************************************
 257
 258 void FGInitialCondition::SetVBodyFpsIC(double tt) {
 259   v=tt;
 260   vt=sqrt(u*u + v*v + w*w);
 261   lastSpeedSet=setuvw;
 262 }
 263
 264 //******************************************************************************
 265
 266 void FGInitialCondition::SetWBodyFpsIC(double tt) {
 267   w=tt;
 268   vt=sqrt( u*u + v*v + w*w );
 269   lastSpeedSet=setuvw;
 270 }
 271
 272 //******************************************************************************
 273
 274 double FGInitialCondition::GetUBodyFpsIC(void) {
 275     if(lastSpeedSet == setvg )
 276       return u;
 277     else
 278       return vt*calpha*cbeta - uw;
 279 }
 280
 281 //******************************************************************************
 282
 283 double FGInitialCondition::GetVBodyFpsIC(void) {
 284     if( lastSpeedSet == setvg )
 285       return v;
 286     else {
 287       return vt*sbeta - vw;
 288     }
 289 }
 290
 291 //******************************************************************************
 292
 293 double FGInitialCondition::GetWBodyFpsIC(void) {
 294     if( lastSpeedSet == setvg )
 295       return w;
 296     else
 297       return vt*salpha*cbeta -ww;
 298 }
 299
 300 //******************************************************************************
 301
 302 void FGInitialCondition::SetWindNEDFpsIC(double wN, double wE, double wD ) {
 303   wnorth = wN; weast = wE; wdown = wD;
 304   lastWindSet = setwned;
 305   calcWindUVW();
 306   if(lastSpeedSet == setvg)
 307     SetVgroundFpsIC(vg);
 308 }
 309
 310 //******************************************************************************
 311
 312 // positive from left
 313 void FGInitialCondition::SetHeadWindKtsIC(double head){
 314     whead=head*ktstofps;
 315     lastWindSet=setwhc;
 316     calcWindUVW();
 317     if(lastSpeedSet == setvg)
 318       SetVgroundFpsIC(vg);
 319
 320 }
 321
 322 //******************************************************************************
 323
 324 void FGInitialCondition::SetCrossWindKtsIC(double cross){
 325     wcross=cross*ktstofps;
 326     lastWindSet=setwhc;
 327     calcWindUVW();
 328     if(lastSpeedSet == setvg)
 329       SetVgroundFpsIC(vg);
 330
 331 }
 332
 333 //******************************************************************************
 334
 335 void FGInitialCondition::SetWindDownKtsIC(double wD) {
 336     wdown=wD;
 337     calcWindUVW();
 338     if(lastSpeedSet == setvg)
 339       SetVgroundFpsIC(vg);
 340 }
 341
 342 //******************************************************************************
 343
 344 void FGInitialCondition::SetWindMagKtsIC(double mag) {
 345   wmag=mag*ktstofps;
 346   lastWindSet=setwmd;
 347   calcWindUVW();
 348   if(lastSpeedSet == setvg)
 349       SetVgroundFpsIC(vg);
 350 }
 351
 352 //******************************************************************************
 353
 354 void FGInitialCondition::SetWindDirDegIC(double dir) {
 355   wdir=dir*degtorad;
 356   lastWindSet=setwmd;
 357   calcWindUVW();
 358   if(lastSpeedSet == setvg)
 359       SetVgroundFpsIC(vg);
 360 }
 361
 362
 363 //******************************************************************************
 364
 365 void FGInitialCondition::calcWindUVW(void) {
 366
 367     switch(lastWindSet) {
 368       case setwmd:
 369         wnorth=wmag*cos(wdir);
 370         weast=wmag*sin(wdir);
 371       break;
 372       case setwhc:
 373         wnorth=whead*cos(psi) + wcross*cos(psi+M_PI/2);
 374         weast=whead*sin(psi) + wcross*sin(psi+M_PI/2);
 375       break;
 376       case setwned:
 377       break;
 378     }
 379     uw=wnorth*ctheta*cpsi +
 380        weast*ctheta*spsi -
 381        wdown*stheta;
 382     vw=wnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 383         weast*( sphi*stheta*spsi + cphi*cpsi ) +
 384        wdown*sphi*ctheta;
 385     ww=wnorth*(cphi*stheta*cpsi + sphi*spsi) +
 386        weast*(cphi*stheta*spsi - sphi*cpsi) +
 387        wdown*cphi*ctheta;
 388
 389
 390     /* cout << "FGInitialCondition::calcWindUVW: wnorth, weast, wdown "
 391          << wnorth << ", " << weast << ", " << wdown << endl;
 392     cout << "FGInitialCondition::calcWindUVW: theta, phi, psi "
 393           << theta << ", " << phi << ", " << psi << endl;
 394     cout << "FGInitialCondition::calcWindUVW: uw, vw, ww "
 395           << uw << ", " << vw << ", " << ww << endl;   */
 396
 397 }
 398
 399 //******************************************************************************
 400
 401 void FGInitialCondition::SetAltitudeFtIC(double tt) {
 402   altitude=tt;
 403   fdmex->GetPosition()->Seth(altitude);
 404   fdmex->GetAtmosphere()->Run();
 405   //lets try to make sure the user gets what they intended
 406
 407   switch(lastSpeedSet) {
 408   case setned:
 409   case setuvw:
 410   case setvt:
 411     SetVtrueKtsIC(vt*fpstokts);
 412     break;
 413   case setvc:
 414     SetVcalibratedKtsIC(vc*fpstokts);
 415     break;
 416   case setve:
 417     SetVequivalentKtsIC(ve*fpstokts);
 418     break;
 419   case setmach:
 420     SetMachIC(mach);
 421     break;
 422   case setvg:
 423     SetVgroundFpsIC(vg);
 424     break;
 425   }
 426 }
 427
 428 //******************************************************************************
 429
 430 void FGInitialCondition::SetAltitudeAGLFtIC(double tt) {
 431   fdmex->GetPosition()->SetDistanceAGL(tt);
 432   altitude=fdmex->GetPosition()->Geth();
 433   SetAltitudeFtIC(altitude);
 434 }
 435
 436 //******************************************************************************
 437
 438 void FGInitialCondition::SetSeaLevelRadiusFtIC(double tt) {
 439   sea_level_radius = tt;
 440 }
 441
 442 //******************************************************************************
 443
 444 void FGInitialCondition::SetTerrainAltitudeFtIC(double tt) {
 445   terrain_altitude=tt;
 446 }
 447
 448 //******************************************************************************
 449
 450 void FGInitialCondition::calcUVWfromNED(void) {
 451   u=vnorth*ctheta*cpsi +
 452      veast*ctheta*spsi -
 453      vdown*stheta;
 454   v=vnorth*( sphi*stheta*cpsi - cphi*spsi ) +
 455      veast*( sphi*stheta*spsi + cphi*cpsi ) +
 456      vdown*sphi*ctheta;
 457   w=vnorth*( cphi*stheta*cpsi + sphi*spsi ) +
 458      veast*( cphi*stheta*spsi - sphi*cpsi ) +
 459      vdown*cphi*ctheta;
 460 }
 461
 462 //******************************************************************************
 463
 464 void FGInitialCondition::SetVnorthFpsIC(double tt) {
 465   vnorth=tt;
 466   calcUVWfromNED();
 467   vt=sqrt(u*u + v*v + w*w);
 468   lastSpeedSet=setned;
 469 }
 470
 471 //******************************************************************************
 472
 473 void FGInitialCondition::SetVeastFpsIC(double tt) {
 474   veast=tt;
 475   calcUVWfromNED();
 476   vt=sqrt(u*u + v*v + w*w);
 477   lastSpeedSet=setned;
 478 }
 479
 480 //******************************************************************************
 481
 482 void FGInitialCondition::SetVdownFpsIC(double tt) {
 483   vdown=tt;
 484   calcUVWfromNED();
 485   vt=sqrt(u*u + v*v + w*w);
 486   SetClimbRateFpsIC(-1*vdown);
 487   lastSpeedSet=setned;
 488 }
 489
 490 //******************************************************************************
 491
 492 bool FGInitialCondition::getMachFromVcas(double *Mach,double vcas) {
 493
 494   bool result=false;
 495   double guess=1.5;
 496   xlo=xhi=0;
 497   xmin=0;xmax=50;
 498   sfunc=&FGInitialCondition::calcVcas;
 499   if(findInterval(vcas,guess)) {
 500     if(solve(&mach,vcas))
 501       result=true;
 502   }
 503   return result;
 504 }
 505
 506 //******************************************************************************
 507
 508 bool FGInitialCondition::getAlpha(void) {
 509   bool result=false;
 510   double guess=theta-gamma;
 511
 512   if(vt < 0.01) return 0;
 513
 514   xlo=xhi=0;
 515   xmin=fdmex->GetAerodynamics()->GetAlphaCLMin();
 516   xmax=fdmex->GetAerodynamics()->GetAlphaCLMax();
 517   sfunc=&FGInitialCondition::GammaEqOfAlpha;
 518   if(findInterval(0,guess)){
 519     if(solve(&alpha,0)){
 520       result=true;
 521       salpha=sin(alpha);
 522       calpha=cos(alpha);
 523     }
 524   }
 525   calcWindUVW();
 526   return result;
 527 }
 528
 529 //******************************************************************************
 530
 531 bool FGInitialCondition::getTheta(void) {
 532   bool result=false;
 533   double guess=alpha+gamma;
 534
 535   if(vt < 0.01) return 0;
 536
 537   xlo=xhi=0;
 538   xmin=-89;xmax=89;
 539   sfunc=&FGInitialCondition::GammaEqOfTheta;
 540   if(findInterval(0,guess)){
 541     if(solve(&theta,0)){
 542       result=true;
 543       stheta=sin(theta);
 544       ctheta=cos(theta);
 545     }
 546   }
 547   calcWindUVW();
 548   return result;
 549 }
 550
 551 //******************************************************************************
 552
 553 double FGInitialCondition::GammaEqOfTheta(double Theta) {
 554   double a,b,c;
 555   double sTheta,cTheta;
 556
 557   //theta=Theta; stheta=sin(theta); ctheta=cos(theta);
 558   sTheta=sin(Theta); cTheta=cos(Theta);
 559   calcWindUVW();
 560   a=wdown + vt*calpha*cbeta + uw;
 561   b=vt*sphi*sbeta + vw*sphi;
 562   c=vt*cphi*salpha*cbeta + ww*cphi;
 563   return vt*sgamma - ( a*sTheta - (b+c)*cTheta);
 564 }
 565
 566 //******************************************************************************
 567
 568 double FGInitialCondition::GammaEqOfAlpha(double Alpha) {
 569   double a,b,c;
 570   double sAlpha,cAlpha;
 571   sAlpha=sin(Alpha); cAlpha=cos(Alpha);
 572   a=wdown + vt*cAlpha*cbeta + uw;
 573   b=vt*sphi*sbeta + vw*sphi;
 574   c=vt*cphi*sAlpha*cbeta + ww*cphi;
 575
 576   return vt*sgamma - ( a*stheta - (b+c)*ctheta );
 577 }
 578
 579 //******************************************************************************
 580
 581 double FGInitialCondition::calcVcas(double Mach) {
 582
 583   double p=fdmex->GetAtmosphere()->GetPressure();
 584   double psl=fdmex->GetAtmosphere()->GetPressureSL();
 585   double rhosl=fdmex->GetAtmosphere()->GetDensitySL();
 586   double pt,A,B,D,vcas;
 587   if(Mach < 0) Mach=0;
 588   if(Mach < 1)    //calculate total pressure assuming isentropic flow
 589     pt=p*pow((1 + 0.2*Mach*Mach),3.5);
 590   else {
 591     // shock in front of pitot tube, we'll assume its normal and use
 592     // the Rayleigh Pitot Tube Formula, i.e. the ratio of total
 593     // pressure behind the shock to the static pressure in front
 594
 595
 596     //the normal shock assumption should not be a bad one -- most supersonic
 597     //aircraft place the pitot probe out front so that it is the forward
 598     //most point on the aircraft.  The real shock would, of course, take
 599     //on something like the shape of a rounded-off cone but, here again,
 600     //the assumption should be good since the opening of the pitot probe
 601     //is very small and, therefore, the effects of the shock curvature
 602     //should be small as well. AFAIK, this approach is fairly well accepted
 603     //within the aerospace community
 604
 605     B = 5.76*Mach*Mach/(5.6*Mach*Mach - 0.8);
 606
 607     // The denominator above is zero for Mach ~ 0.38, for which
 608     // we'll never be here, so we're safe
 609
 610     D = (2.8*Mach*Mach-0.4)*0.4167;
 611     pt = p*pow(B,3.5)*D;
 612   }
 613
 614   A = pow(((pt-p)/psl+1),0.28571);
 615   vcas = sqrt(7*psl/rhosl*(A-1));
 616   //cout << "calcVcas: vcas= " << vcas*fpstokts << " mach= " << Mach << " pressure: " << pt << endl;
 617   return vcas;
 618 }
 619
 620 //******************************************************************************
 621
 622 bool FGInitialCondition::findInterval(double x,double guess) {
 623   //void find_interval(inter_params &ip,eqfunc f,double y,double constant, int &flag){
 624
 625   int i=0;
 626   bool found=false;
 627   double flo,fhi,fguess;
 628   double lo,hi,step;
 629   step=0.1;
 630   fguess=(this->*sfunc)(guess)-x;
 631   lo=hi=guess;
 632   do {
 633     step=2*step;
 634     lo-=step;
 635     hi+=step;
 636     if(lo < xmin) lo=xmin;
 637     if(hi > xmax) hi=xmax;
 638     i++;
 639     flo=(this->*sfunc)(lo)-x;
 640     fhi=(this->*sfunc)(hi)-x;
 641     if(flo*fhi <=0) {  //found interval with root
 642       found=true;
 643       if(flo*fguess <= 0) {  //narrow interval down a bit
 644         hi=lo+step;    //to pass solver interval that is as
 645         //small as possible
 646       }
 647       else if(fhi*fguess <= 0) {
 648         lo=hi-step;
 649       }
 650     }
 651     //cout << "findInterval: i=" << i << " Lo= " << lo << " Hi= " << hi << endl;
 652   }
 653   while((found == 0) && (i <= 100));
 654   xlo=lo;
 655   xhi=hi;
 656   return found;
 657 }
 658
 659 //******************************************************************************
 660
 661 bool FGInitialCondition::solve(double *y,double x)
 662 {
 663   double x1,x2,x3,f1,f2,f3,d,d0;
 664   double eps=1E-5;
 665   double const relax =0.9;
 666   int i;
 667   bool success=false;
 668
 669   //initializations
 670   d=1;
 671   x2 = 0;
 672   x1=xlo;x3=xhi;
 673   f1=(this->*sfunc)(x1)-x;
 674   f3=(this->*sfunc)(x3)-x;
 675   d0=fabs(x3-x1);
 676
 677   //iterations
 678   i=0;
 679   while ((fabs(d) > eps) && (i < 100)) {
 680     d=(x3-x1)/d0;
 681     x2 = x1-d*d0*f1/(f3-f1);
 682
 683     f2=(this->*sfunc)(x2)-x;
 684     //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
 685     //cout << "                " << f1 << "," << f2 << "," << f3 << endl;
 686
 687     if(fabs(f2) <= 0.001) {
 688       x1=x3=x2;
 689     } else if(f1*f2 <= 0.0) {
 690       x3=x2;
 691       f3=f2;
 692       f1=relax*f1;
 693     } else if(f2*f3 <= 0) {
 694       x1=x2;
 695       f1=f2;
 696       f3=relax*f3;
 697     }
 698     //cout << i << endl;
 699     i++;
 700   }//end while
 701   if(i < 100) {
 702     success=true;
 703     *y=x2;
 704   }
 705
 706   //cout << "Success= " << success << " Vcas: " << vcas*fpstokts << " Mach: " << x2 << endl;
 707   return success;
 708 }
 709
 710 //******************************************************************************
 711
 712 double FGInitialCondition::GetWindDirDegIC(void) {
 713   if(weast != 0.0)
 714     return atan2(weast,wnorth)*radtodeg;
 715   else if(wnorth > 0)
 716     return 0.0;
 717   else
 718     return 180.0;
 719 }
 720
 721 //******************************************************************************
 722
 723 bool FGInitialCondition::Load(string rstfile, bool useStoredPath)
 724 {
 725   string resetDef;
 726   string token="";
 727
 728   double temp;
 729   # ifndef macintosh
 730     string sep = "/";
 731   # else
 732     string sep = ";";
 733   #endif
 734
 735   if( useStoredPath ) {
 736     string acpath = fdmex->GetAircraftPath() + sep + fdmex->GetModelName();
 737     resetDef = acpath + sep + rstfile + ".xml";
 738   } else {
 739     resetDef = rstfile;
 740   }
 741
 742   FGConfigFile resetfile(resetDef);
 743   if (!resetfile.IsOpen()) {
 744     cerr << "Failed to open reset file: " << resetDef << endl;
 745     return false;
 746   }
 747
 748   resetfile.GetNextConfigLine();
 749   token = resetfile.GetValue();
 750   if (token != string("initialize")) {
 751     cerr << "The reset file " << resetDef
 752          << " does not appear to be a reset file" << endl;
 753     return false;
 754   }
 755
 756   resetfile.GetNextConfigLine();
 757   resetfile >> token;
 758   while (token != string("/initialize") && token != string("EOF")) {
 759     if (token == "UBODY" ) { resetfile >> temp; SetUBodyFpsIC(temp); }
 760     if (token == "VBODY" ) { resetfile >> temp; SetVBodyFpsIC(temp); }
 761     if (token == "WBODY" ) { resetfile >> temp; SetWBodyFpsIC(temp); }
 762     if (token == "LATITUDE" ) { resetfile >> temp; SetLatitudeDegIC(temp); }
 763     if (token == "LONGITUDE" ) { resetfile >> temp; SetLongitudeDegIC(temp); }
 764     if (token == "PHI" ) { resetfile >> temp; SetRollAngleDegIC(temp); }
 765     if (token == "THETA" ) { resetfile >> temp; SetPitchAngleDegIC(temp); }
 766     if (token == "PSI" ) { resetfile >> temp; SetTrueHeadingDegIC(temp); }
 767     if (token == "ALPHA" ) { resetfile >> temp; SetAlphaDegIC(temp); }
 768     if (token == "BETA" ) { resetfile >> temp; SetBetaDegIC(temp); }
 769     if (token == "GAMMA" ) { resetfile >> temp; SetFlightPathAngleDegIC(temp); }
 770     if (token == "ROC" ) { resetfile >> temp; SetClimbRateFpmIC(temp); }
 771     if (token == "ALTITUDE" ) { resetfile >> temp; SetAltitudeFtIC(temp); }
 772     if (token == "WINDDIR" ) { resetfile >> temp; SetWindDirDegIC(temp); }
 773     if (token == "VWIND" ) { resetfile >> temp; SetWindMagKtsIC(temp); }
 774     if (token == "HWIND" ) { resetfile >> temp; SetHeadWindKtsIC(temp); }
 775     if (token == "XWIND" ) { resetfile >> temp; SetCrossWindKtsIC(temp); }
 776     if (token == "VC" ) { resetfile >> temp; SetVcalibratedKtsIC(temp); }
 777     if (token == "MACH" ) { resetfile >> temp; SetMachIC(temp); }
 778     if (token == "VGROUND" ) { resetfile >> temp; SetVgroundKtsIC(temp); }
 779     resetfile >> token;
 780   }
 781
 782   fdmex->RunIC();
 783
 784   return true;
 785 }
 786
 787 //******************************************************************************
 788
 789 void FGInitialCondition::bind(void){
 790   PropertyManager->Tie("ic/vc-kts", this,
 791                        &FGInitialCondition::GetVcalibratedKtsIC,
 792                        &FGInitialCondition::SetVcalibratedKtsIC,
 793                        true);
 794   PropertyManager->Tie("ic/ve-kts", this,
 795                        &FGInitialCondition::GetVequivalentKtsIC,
 796                        &FGInitialCondition::SetVequivalentKtsIC,
 797                        true);
 798   PropertyManager->Tie("ic/vg-kts", this,
 799                        &FGInitialCondition::GetVgroundKtsIC,
 800                        &FGInitialCondition::SetVgroundKtsIC,
 801                        true);
 802   PropertyManager->Tie("ic/vt-kts", this,
 803                        &FGInitialCondition::GetVtrueKtsIC,
 804                        &FGInitialCondition::SetVtrueKtsIC,
 805                        true);
 806   PropertyManager->Tie("ic/mach-norm", this,
 807                        &FGInitialCondition::GetMachIC,
 808                        &FGInitialCondition::SetMachIC,
 809                        true);
 810   PropertyManager->Tie("ic/roc-fpm", this,
 811                        &FGInitialCondition::GetClimbRateFpmIC,
 812                        &FGInitialCondition::SetClimbRateFpmIC,
 813                        true);
 814   PropertyManager->Tie("ic/gamma-deg", this,
 815                        &FGInitialCondition::GetFlightPathAngleDegIC,
 816                        &FGInitialCondition::SetFlightPathAngleDegIC,
 817                        true);
 818   PropertyManager->Tie("ic/alpha-deg", this,
 819                        &FGInitialCondition::GetAlphaDegIC,
 820                        &FGInitialCondition::SetAlphaDegIC,
 821                        true);
 822   PropertyManager->Tie("ic/beta-deg", this,
 823                        &FGInitialCondition::GetBetaDegIC,
 824                        &FGInitialCondition::SetBetaDegIC,
 825                        true);
 826   PropertyManager->Tie("ic/theta-deg", this,
 827                        &FGInitialCondition::GetPitchAngleDegIC,
 828                        &FGInitialCondition::SetPitchAngleDegIC,
 829                        true);
 830   PropertyManager->Tie("ic/phi-deg", this,
 831                        &FGInitialCondition::GetRollAngleDegIC,
 832                        &FGInitialCondition::SetRollAngleDegIC,
 833                        true);
 834   PropertyManager->Tie("ic/psi-true-deg", this,
 835                        &FGInitialCondition::GetHeadingDegIC );
 836   PropertyManager->Tie("ic/lat-gc-deg", this,
 837                        &FGInitialCondition::GetLatitudeDegIC,
 838                        &FGInitialCondition::SetLatitudeDegIC,
 839                        true);
 840   PropertyManager->Tie("ic/long-gc-deg", this,
 841                        &FGInitialCondition::GetLongitudeDegIC,
 842                        &FGInitialCondition::SetLongitudeDegIC,
 843                        true);
 844   PropertyManager->Tie("ic/h-sl-ft", this,
 845                        &FGInitialCondition::GetAltitudeFtIC,
 846                        &FGInitialCondition::SetAltitudeFtIC,
 847                        true);
 848   PropertyManager->Tie("ic/h-agl-ft", this,
 849                        &FGInitialCondition::GetAltitudeAGLFtIC,
 850                        &FGInitialCondition::SetAltitudeAGLFtIC,
 851                        true);
 852   PropertyManager->Tie("ic/sea-level-radius-ft", this,
 853                        &FGInitialCondition::GetSeaLevelRadiusFtIC,
 854                        &FGInitialCondition::SetSeaLevelRadiusFtIC,
 855                        true);
 856   PropertyManager->Tie("ic/terrain-altitude-ft", this,
 857                        &FGInitialCondition::GetTerrainAltitudeFtIC,
 858                        &FGInitialCondition::SetTerrainAltitudeFtIC,
 859                        true);
 860   PropertyManager->Tie("ic/vg-fps", this,
 861                        &FGInitialCondition::GetVgroundFpsIC,
 862                        &FGInitialCondition::SetVgroundFpsIC,
 863                        true);
 864   PropertyManager->Tie("ic/vt-fps", this,
 865                        &FGInitialCondition::GetVtrueFpsIC,
 866                        &FGInitialCondition::SetVtrueFpsIC,
 867                        true);
 868   PropertyManager->Tie("ic/vw-bx-fps", this,
 869                        &FGInitialCondition::GetWindUFpsIC);
 870   PropertyManager->Tie("ic/vw-by-fps", this,
 871                        &FGInitialCondition::GetWindVFpsIC);
 872   PropertyManager->Tie("ic/vw-bz-fps", this,
 873                        &FGInitialCondition::GetWindWFpsIC);
 874   PropertyManager->Tie("ic/vw-north-fps", this,
 875                        &FGInitialCondition::GetWindNFpsIC);
 876   PropertyManager->Tie("ic/vw-east-fps", this,
 877                        &FGInitialCondition::GetWindEFpsIC);
 878   PropertyManager->Tie("ic/vw-down-fps", this,
 879                        &FGInitialCondition::GetWindDFpsIC);
 880   PropertyManager->Tie("ic/vw-mag-fps", this,
 881                        &FGInitialCondition::GetWindFpsIC);
 882  /*  PropertyManager->Tie("ic/vw-dir-deg", this,
 883                        &FGInitialCondition::GetWindDirDegIC,
 884                        &FGInitialCondition::SetWindDirDegIC,
 885                        true); */
 886
 887   PropertyManager->Tie("ic/roc-fps", this,
 888                        &FGInitialCondition::GetClimbRateFpsIC,
 889                        &FGInitialCondition::SetClimbRateFpsIC,
 890                        true);
 891   /* PropertyManager->Tie("ic/u-fps", this,
 892                        &FGInitialCondition::GetUBodyFpsIC,
 893                        &FGInitialCondition::SetUBodyFpsIC,
 894                        true);
 895   PropertyManager->Tie("ic/v-fps", this,
 896                        &FGInitialCondition::GetVBodyFpsIC,
 897                        &FGInitialCondition::SetVBodyFpsIC,
 898                        true);
 899   PropertyManager->Tie("ic/w-fps", this,
 900                        &FGInitialCondition::GetWBodyFpsIC,
 901                        &FGInitialCondition::SetWBodyFpsIC,
 902                        true); */
 903
 904   PropertyManager->Tie("ic/gamma-rad", this,
 905                        &FGInitialCondition::GetFlightPathAngleRadIC,
 906                        &FGInitialCondition::SetFlightPathAngleRadIC,
 907                        true);
 908   PropertyManager->Tie("ic/alpha-rad", this,
 909                        &FGInitialCondition::GetAlphaRadIC,
 910                        &FGInitialCondition::SetAlphaRadIC,
 911                        true);
 912   PropertyManager->Tie("ic/theta-rad", this,
 913                        &FGInitialCondition::GetPitchAngleRadIC,
 914                        &FGInitialCondition::SetPitchAngleRadIC,
 915                        true);
 916   PropertyManager->Tie("ic/beta-rad", this,
 917                        &FGInitialCondition::GetBetaRadIC,
 918                        &FGInitialCondition::SetBetaRadIC,
 919                        true);
 920   PropertyManager->Tie("ic/phi-rad", this,
 921                        &FGInitialCondition::GetRollAngleRadIC,
 922                        &FGInitialCondition::SetRollAngleRadIC,
 923                        true);
 924   PropertyManager->Tie("ic/psi-true-rad", this,
 925                        &FGInitialCondition::GetHeadingRadIC);
 926   PropertyManager->Tie("ic/lat-gc-rad", this,
 927                        &FGInitialCondition::GetLatitudeRadIC,
 928                        &FGInitialCondition::SetLatitudeRadIC,
 929                        true);
 930   PropertyManager->Tie("ic/long-gc-rad", this,
 931                        &FGInitialCondition::GetLongitudeRadIC,
 932                        &FGInitialCondition::SetLongitudeRadIC,
 933                        true);
 934 }
 935
 936 //******************************************************************************
 937
 938 void FGInitialCondition::unbind(void){
 939   PropertyManager->Untie("ic/vc-kts");
 940   PropertyManager->Untie("ic/ve-kts");
 941   PropertyManager->Untie("ic/vg-kts");
 942   PropertyManager->Untie("ic/vt-kts");
 943   PropertyManager->Untie("ic/mach-norm");
 944   PropertyManager->Untie("ic/roc-fpm");
 945   PropertyManager->Untie("ic/gamma-deg");
 946   PropertyManager->Untie("ic/alpha-deg");
 947   PropertyManager->Untie("ic/beta-deg");
 948   PropertyManager->Untie("ic/theta-deg");
 949   PropertyManager->Untie("ic/phi-deg");
 950   PropertyManager->Untie("ic/psi-true-deg");
 951   PropertyManager->Untie("ic/lat-gc-deg");
 952   PropertyManager->Untie("ic/long-gc-deg");
 953   PropertyManager->Untie("ic/h-sl-ft");
 954   PropertyManager->Untie("ic/h-agl-ft");
 955   PropertyManager->Untie("ic/sea-level-radius-ft");
 956   PropertyManager->Untie("ic/terrain-altitude-ft");
 957   PropertyManager->Untie("ic/vg-fps");
 958   PropertyManager->Untie("ic/vt-fps");
 959   PropertyManager->Untie("ic/vw-bx-fps");
 960   PropertyManager->Untie("ic/vw-by-fps");
 961   PropertyManager->Untie("ic/vw-bz-fps");
 962   PropertyManager->Untie("ic/vw-north-fps");
 963   PropertyManager->Untie("ic/vw-east-fps");
 964   PropertyManager->Untie("ic/vw-down-fps");
 965   PropertyManager->Untie("ic/vw-mag-fps");
 966   /* PropertyManager->Untie("ic/vw-dir-deg"); */
 967
 968   PropertyManager->Untie("ic/roc-fps");
 969
 970   /*  PropertyManager->Untie("ic/u-fps");
 971   PropertyManager->Untie("ic/v-fps");
 972   PropertyManager->Untie("ic/w-fps"); */
 973
 974   PropertyManager->Untie("ic/gamma-rad");
 975   PropertyManager->Untie("ic/alpha-rad");
 976   PropertyManager->Untie("ic/theta-rad");
 977   PropertyManager->Untie("ic/beta-rad");
 978   PropertyManager->Untie("ic/phi-rad");
 979   PropertyManager->Untie("ic/psi-true-rad");
 980   PropertyManager->Untie("ic/lat-gc-rad");
 981   PropertyManager->Untie("ic/long-gc-rad");
 982 }
 983
 984 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 985 //    The bitmasked value choices are as follows:
 986 //    unset: In this case (the default) JSBSim would only print
 987 //       out the normally expected messages, essentially echoing
 988 //       the config files as they are read. If the environment
 989 //       variable is not set, debug_lvl is set to 1 internally
 990 //    0: This requests JSBSim not to output any messages
 991 //       whatsoever.
 992 //    1: This value explicity requests the normal JSBSim
 993 //       startup messages
 994 //    2: This value asks for a message to be printed out when
 995 //       a class is instantiated
 996 //    4: When this value is set, a message is displayed when a
 997 //       FGModel object executes its Run() method
 998 //    8: When this value is set, various runtime state variables
 999 //       are printed out periodically
1000 //    16: When set various parameters are sanity checked and
1001 //       a message is printed out when they go out of bounds
1002
1003 void FGInitialCondition::Debug(int from)
1004 {
1005   if (debug_lvl <= 0) return;
1006
1007   if (debug_lvl & 1) { // Standard console startup message output
1008   }
1009   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
1010     if (from == 0) cout << "Instantiated: FGInitialCondition" << endl;
1011     if (from == 1) cout << "Destroyed:    FGInitialCondition" << endl;
1012   }
1013   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
1014   }
1015   if (debug_lvl & 8 ) { // Runtime state variables
1016   }
1017   if (debug_lvl & 16) { // Sanity checking
1018   }
1019   if (debug_lvl & 64) {
1020     if (from == 0) { // Constructor
1021       cout << IdSrc << endl;
1022       cout << IdHdr << endl;
1023     }
1024   }
1025 }
1026 }