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