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