src/FDM/JSBSim/FGInitialCondition.cpp

   1 /*******************************************************************************
   2
   3  Header:       FGInitialCondition.cpp
   4  Author:       Tony Peden
   5  Date started: 7/1/99
   6
   7  ------------- Copyright (C) 1999  Anthony K. Peden (apeden@earthlink.net) -------------
   8
   9  This program is free software; you can redistribute it and/or modify it under
  10  the terms of the GNU General Public License as published by the Free Software
  11  Foundation; either version 2 of the License, or (at your option) any later
  12  version.
  13
  14  This program is distributed in the hope that it will be useful, but WITHOUT
  15  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  16  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  17  details.
  18
  19  You should have received a copy of the GNU General Public License along with
  20  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  21  Place - Suite 330, Boston, MA  02111-1307, USA.
  22
  23  Further information about the GNU General Public License can also be found on
  24  the world wide web at http://www.gnu.org.
  25
  26
  27  HISTORY
  28 --------------------------------------------------------------------------------
  29 7/1/99   TP   Created
  30
  31
  32 FUNCTIONAL DESCRIPTION
  33 --------------------------------------------------------------------------------
  34
  35 The purpose of this class is to take a set of initial conditions and provide
  36 a kinematically consistent set of body axis velocity components, euler
  37 angles, and altitude.  This class does not attempt to trim the model i.e.
  38 the sim will most likely start in a very dynamic state (unless, of course,
  39 you have chosen your IC's wisely) even after setting it up with this class.
  40
  41 ********************************************************************************
  42 INCLUDES
  43 *******************************************************************************/
  44
  45 #include "FGInitialCondition.h"
  46 #include "FGFDMExec.h"
  47 #include "FGState.h"
  48 #include "FGAtmosphere.h"
  49 #include "FGFCS.h"
  50 #include "FGAircraft.h"
  51 #include "FGTranslation.h"
  52 #include "FGRotation.h"
  53 #include "FGPosition.h"
  54 #include "FGAuxiliary.h"
  55 #include "FGOutput.h"
  56 #include "FGConfigFile.h"
  57 #include "FGPropertyManager.h"
  58
  59 static const char *IdSrc = "$Id$";
  60 static const char *IdHdr = ID_INITIALCONDITION;
  61
  62 //******************************************************************************
  63
  64 FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec)
  65 {
  66   vt=vc=ve=vg=0;
  67   mach=0;
  68   alpha=beta=gamma=0;
  69   theta=phi=psi=0;
  70   altitude=hdot=0;
  71   latitude=longitude=0;
  72   u=v=w=0;
  73   uw=vw=ww=0;
  74   vnorth=veast=vdown=0;
  75   wnorth=weast=wdown=0;
  76   whead=wcross=0;
  77   wdir=wmag=0;
  78   lastSpeedSet=setvt;
  79   lastWindSet=setwned;
  80   sea_level_radius = FDMExec->GetInertial()->RefRadius();
  81   radius_to_vehicle = FDMExec->GetInertial()->RefRadius();
  82   terrain_altitude = 0;
  83
  84   salpha=sbeta=stheta=sphi=spsi=sgamma=0;
  85   calpha=cbeta=ctheta=cphi=cpsi=cgamma=1;
  86
  87   if(FDMExec != NULL ) {
  88     fdmex=FDMExec;
  89     fdmex->GetPosition()->Seth(altitude);
  90     fdmex->GetAtmosphere()->Run();
  91     PropertyManager=fdmex->GetPropertyManager();
  92     bind();
  93   } else {
  94     cout << "FGInitialCondition: This class requires a pointer to a valid FGFDMExec object" << endl;
  95   }
  96
  97   Debug(0);
  98 }
  99
 100 //******************************************************************************
 101
 102 FGInitialCondition::~FGInitialCondition()
 103 {
 104   unbind();
 105   Debug(1);
 106 }
 107
 108 //******************************************************************************
 109
 110 void FGInitialCondition::SetVcalibratedKtsIC(double tt) {
 111
 112   if(getMachFromVcas(&mach,tt*ktstofps)) {
 113
 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 acpath, string acname, string rstfile)
 722 {
 723   string resetDef;
 724   string token="";
 725
 726   double temp;
 727
 728 # ifndef macintosh
 729   resetDef = acpath + "/" + acname + "/" + rstfile + ".xml";
 730 # else
 731   resetDef = acpath + ";" + acname + ";" + rstfile + ".xml";
 732 # endif
 733
 734   FGConfigFile resetfile(resetDef);
 735   if (!resetfile.IsOpen()) {
 736     cerr << "Failed to open reset file: " << resetDef << endl;
 737     return false;
 738   }
 739
 740   resetfile.GetNextConfigLine();
 741   token = resetfile.GetValue();
 742   if (token != string("initialize")) {
 743     cerr << "The reset file " << resetDef
 744          << " does not appear to be a reset file" << endl;
 745     return false;
 746   }
 747
 748   resetfile.GetNextConfigLine();
 749   resetfile >> token;
 750   while (token != string("/initialize") && token != string("EOF")) {
 751     if (token == "UBODY" ) { resetfile >> temp; SetUBodyFpsIC(temp); }
 752     if (token == "VBODY" ) { resetfile >> temp; SetVBodyFpsIC(temp); }
 753     if (token == "WBODY" ) { resetfile >> temp; SetWBodyFpsIC(temp); }
 754     if (token == "LATITUDE" ) { resetfile >> temp; SetLatitudeDegIC(temp); }
 755     if (token == "LONGITUDE" ) { resetfile >> temp; SetLongitudeDegIC(temp); }
 756     if (token == "PHI" ) { resetfile >> temp; SetRollAngleDegIC(temp); }
 757     if (token == "THETA" ) { resetfile >> temp; SetPitchAngleDegIC(temp); }
 758     if (token == "PSI" ) { resetfile >> temp; SetTrueHeadingDegIC(temp); }
 759     if (token == "ALPHA" ) { resetfile >> temp; SetAlphaDegIC(temp); }
 760     if (token == "BETA" ) { resetfile >> temp; SetBetaDegIC(temp); }
 761     if (token == "GAMMA" ) { resetfile >> temp; SetFlightPathAngleDegIC(temp); }
 762     if (token == "ROC" ) { resetfile >> temp; SetClimbRateFpmIC(temp); }
 763     if (token == "ALTITUDE" ) { resetfile >> temp; SetAltitudeFtIC(temp); }
 764     if (token == "WINDDIR" ) { resetfile >> temp; SetWindDirDegIC(temp); }
 765     if (token == "VWIND" ) { resetfile >> temp; SetWindMagKtsIC(temp); }
 766     if (token == "HWIND" ) { resetfile >> temp; SetHeadWindKtsIC(temp); }
 767     if (token == "XWIND" ) { resetfile >> temp; SetCrossWindKtsIC(temp); }
 768     if (token == "VC" ) { resetfile >> temp; SetVcalibratedKtsIC(temp); }
 769     if (token == "MACH" ) { resetfile >> temp; SetMachIC(temp); }
 770     if (token == "VGROUND" ) { resetfile >> temp; SetVgroundKtsIC(temp); }
 771     resetfile >> token;
 772   }
 773
 774   fdmex->RunIC(this);
 775
 776   return true;
 777 }
 778
 779 //******************************************************************************
 780
 781 void FGInitialCondition::bind(void){
 782   PropertyManager->Tie("ic/vc-kts", this,
 783                        &FGInitialCondition::GetVcalibratedKtsIC,
 784                        &FGInitialCondition::SetVcalibratedKtsIC,
 785                        true);
 786   PropertyManager->Tie("ic/ve-kts", this,
 787                        &FGInitialCondition::GetVequivalentKtsIC,
 788                        &FGInitialCondition::SetVequivalentKtsIC,
 789                        true);
 790   PropertyManager->Tie("ic/vg-kts", this,
 791                        &FGInitialCondition::GetVgroundKtsIC,
 792                        &FGInitialCondition::SetVgroundKtsIC,
 793                        true);
 794   PropertyManager->Tie("ic/vt-kts", this,
 795                        &FGInitialCondition::GetVtrueKtsIC,
 796                        &FGInitialCondition::SetVtrueKtsIC,
 797                        true);
 798   PropertyManager->Tie("ic/mach-norm", this,
 799                        &FGInitialCondition::GetMachIC,
 800                        &FGInitialCondition::SetMachIC,
 801                        true);
 802   PropertyManager->Tie("ic/roc-fpm", this,
 803                        &FGInitialCondition::GetClimbRateFpmIC,
 804                        &FGInitialCondition::SetClimbRateFpmIC,
 805                        true);
 806   PropertyManager->Tie("ic/gamma-deg", this,
 807                        &FGInitialCondition::GetFlightPathAngleDegIC,
 808                        &FGInitialCondition::SetFlightPathAngleDegIC,
 809                        true);
 810   PropertyManager->Tie("ic/alpha-deg", this,
 811                        &FGInitialCondition::GetAlphaDegIC,
 812                        &FGInitialCondition::SetAlphaDegIC,
 813                        true);
 814   PropertyManager->Tie("ic/beta-deg", this,
 815                        &FGInitialCondition::GetBetaDegIC,
 816                        &FGInitialCondition::SetBetaDegIC,
 817                        true);
 818   PropertyManager->Tie("ic/theta-deg", this,
 819                        &FGInitialCondition::GetPitchAngleDegIC,
 820                        &FGInitialCondition::SetPitchAngleDegIC,
 821                        true);
 822   PropertyManager->Tie("ic/phi-deg", this,
 823                        &FGInitialCondition::GetRollAngleDegIC,
 824                        &FGInitialCondition::SetRollAngleDegIC,
 825                        true);
 826   PropertyManager->Tie("ic/psi-true-deg", this,
 827                        &FGInitialCondition::GetHeadingDegIC );
 828   PropertyManager->Tie("ic/lat-gc-deg", this,
 829                        &FGInitialCondition::GetLatitudeDegIC,
 830                        &FGInitialCondition::SetLatitudeDegIC,
 831                        true);
 832   PropertyManager->Tie("ic/long-gc-deg", this,
 833                        &FGInitialCondition::GetLongitudeDegIC,
 834                        &FGInitialCondition::SetLongitudeDegIC,
 835                        true);
 836   PropertyManager->Tie("ic/h-sl-ft", this,
 837                        &FGInitialCondition::GetAltitudeFtIC,
 838                        &FGInitialCondition::SetAltitudeFtIC,
 839                        true);
 840   PropertyManager->Tie("ic/h-agl-ft", this,
 841                        &FGInitialCondition::GetAltitudeAGLFtIC,
 842                        &FGInitialCondition::SetAltitudeAGLFtIC,
 843                        true);
 844   PropertyManager->Tie("ic/sea-level-radius-ft", this,
 845                        &FGInitialCondition::GetSeaLevelRadiusFtIC,
 846                        &FGInitialCondition::SetSeaLevelRadiusFtIC,
 847                        true);
 848   PropertyManager->Tie("ic/terrain-altitude-ft", this,
 849                        &FGInitialCondition::GetTerrainAltitudeFtIC,
 850                        &FGInitialCondition::SetTerrainAltitudeFtIC,
 851                        true);
 852   PropertyManager->Tie("ic/vg-fps", this,
 853                        &FGInitialCondition::GetVgroundFpsIC,
 854                        &FGInitialCondition::SetVgroundFpsIC,
 855                        true);
 856   PropertyManager->Tie("ic/vt-fps", this,
 857                        &FGInitialCondition::GetVtrueFpsIC,
 858                        &FGInitialCondition::SetVtrueFpsIC,
 859                        true);
 860   PropertyManager->Tie("ic/vw-bx-fps", this,
 861                        &FGInitialCondition::GetWindUFpsIC);
 862   PropertyManager->Tie("ic/vw-by-fps", this,
 863                        &FGInitialCondition::GetWindVFpsIC);
 864   PropertyManager->Tie("ic/vw-bz-fps", this,
 865                        &FGInitialCondition::GetWindWFpsIC);
 866   PropertyManager->Tie("ic/vw-north-fps", this,
 867                        &FGInitialCondition::GetWindNFpsIC);
 868   PropertyManager->Tie("ic/vw-east-fps", this,
 869                        &FGInitialCondition::GetWindEFpsIC);
 870   PropertyManager->Tie("ic/vw-down-fps", this,
 871                        &FGInitialCondition::GetWindDFpsIC);
 872   PropertyManager->Tie("ic/vw-mag-fps", this,
 873                        &FGInitialCondition::GetWindFpsIC);
 874  /*  PropertyManager->Tie("ic/vw-dir-deg", this,
 875                        &FGInitialCondition::GetWindDirDegIC,
 876                        &FGInitialCondition::SetWindDirDegIC,
 877                        true); */
 878
 879   PropertyManager->Tie("ic/roc-fps", this,
 880                        &FGInitialCondition::GetClimbRateFpsIC,
 881                        &FGInitialCondition::SetClimbRateFpsIC,
 882                        true);
 883   /* PropertyManager->Tie("ic/u-fps", this,
 884                        &FGInitialCondition::GetUBodyFpsIC,
 885                        &FGInitialCondition::SetUBodyFpsIC,
 886                        true);
 887   PropertyManager->Tie("ic/v-fps", this,
 888                        &FGInitialCondition::GetVBodyFpsIC,
 889                        &FGInitialCondition::SetVBodyFpsIC,
 890                        true);
 891   PropertyManager->Tie("ic/w-fps", this,
 892                        &FGInitialCondition::GetWBodyFpsIC,
 893                        &FGInitialCondition::SetWBodyFpsIC,
 894                        true); */
 895
 896   PropertyManager->Tie("ic/gamma-rad", this,
 897                        &FGInitialCondition::GetFlightPathAngleRadIC,
 898                        &FGInitialCondition::SetFlightPathAngleRadIC,
 899                        true);
 900   PropertyManager->Tie("ic/alpha-rad", this,
 901                        &FGInitialCondition::GetAlphaRadIC,
 902                        &FGInitialCondition::SetAlphaRadIC,
 903                        true);
 904   PropertyManager->Tie("ic/theta-rad", this,
 905                        &FGInitialCondition::GetPitchAngleRadIC,
 906                        &FGInitialCondition::SetPitchAngleRadIC,
 907                        true);
 908   PropertyManager->Tie("ic/beta-rad", this,
 909                        &FGInitialCondition::GetBetaRadIC,
 910                        &FGInitialCondition::SetBetaRadIC,
 911                        true);
 912   PropertyManager->Tie("ic/phi-rad", this,
 913                        &FGInitialCondition::GetRollAngleRadIC,
 914                        &FGInitialCondition::SetRollAngleRadIC,
 915                        true);
 916   PropertyManager->Tie("ic/psi-true-rad", this,
 917                        &FGInitialCondition::GetHeadingRadIC);
 918   PropertyManager->Tie("ic/lat-gc-rad", this,
 919                        &FGInitialCondition::GetLatitudeRadIC,
 920                        &FGInitialCondition::SetLatitudeRadIC,
 921                        true);
 922   PropertyManager->Tie("ic/long-gc-rad", this,
 923                        &FGInitialCondition::GetLongitudeRadIC,
 924                        &FGInitialCondition::SetLongitudeRadIC,
 925                        true);
 926 }
 927
 928 //******************************************************************************
 929
 930 void FGInitialCondition::unbind(void){
 931   PropertyManager->Untie("ic/vc-kts");
 932   PropertyManager->Untie("ic/ve-kts");
 933   PropertyManager->Untie("ic/vg-kts");
 934   PropertyManager->Untie("ic/vt-kts");
 935   PropertyManager->Untie("ic/mach-norm");
 936   PropertyManager->Untie("ic/roc-fpm");
 937   PropertyManager->Untie("ic/gamma-deg");
 938   PropertyManager->Untie("ic/alpha-deg");
 939   PropertyManager->Untie("ic/beta-deg");
 940   PropertyManager->Untie("ic/theta-deg");
 941   PropertyManager->Untie("ic/phi-deg");
 942   PropertyManager->Untie("ic/psi-true-deg");
 943   PropertyManager->Untie("ic/lat-gc-deg");
 944   PropertyManager->Untie("ic/long-gc-deg");
 945   PropertyManager->Untie("ic/h-sl-ft");
 946   PropertyManager->Untie("ic/h-agl-ft");
 947   PropertyManager->Untie("ic/sea-level-radius-ft");
 948   PropertyManager->Untie("ic/terrain-altitude-ft");
 949   PropertyManager->Untie("ic/vg-fps");
 950   PropertyManager->Untie("ic/vt-fps");
 951   PropertyManager->Untie("ic/vw-bx-fps");
 952   PropertyManager->Untie("ic/vw-by-fps");
 953   PropertyManager->Untie("ic/vw-bz-fps");
 954   PropertyManager->Untie("ic/vw-north-fps");
 955   PropertyManager->Untie("ic/vw-east-fps");
 956   PropertyManager->Untie("ic/vw-down-fps");
 957   PropertyManager->Untie("ic/vw-mag-fps");
 958   /* PropertyManager->Untie("ic/vw-dir-deg"); */
 959
 960   PropertyManager->Untie("ic/roc-fps");
 961
 962   /*  PropertyManager->Untie("ic/u-fps");
 963   PropertyManager->Untie("ic/v-fps");
 964   PropertyManager->Untie("ic/w-fps"); */
 965
 966   PropertyManager->Untie("ic/gamma-rad");
 967   PropertyManager->Untie("ic/alpha-rad");
 968   PropertyManager->Untie("ic/theta-rad");
 969   PropertyManager->Untie("ic/beta-rad");
 970   PropertyManager->Untie("ic/phi-rad");
 971   PropertyManager->Untie("ic/psi-true-rad");
 972   PropertyManager->Untie("ic/lat-gc-rad");
 973   PropertyManager->Untie("ic/long-gc-rad");
 974 }
 975
 976 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 977 //    The bitmasked value choices are as follows:
 978 //    unset: In this case (the default) JSBSim would only print
 979 //       out the normally expected messages, essentially echoing
 980 //       the config files as they are read. If the environment
 981 //       variable is not set, debug_lvl is set to 1 internally
 982 //    0: This requests JSBSim not to output any messages
 983 //       whatsoever.
 984 //    1: This value explicity requests the normal JSBSim
 985 //       startup messages
 986 //    2: This value asks for a message to be printed out when
 987 //       a class is instantiated
 988 //    4: When this value is set, a message is displayed when a
 989 //       FGModel object executes its Run() method
 990 //    8: When this value is set, various runtime state variables
 991 //       are printed out periodically
 992 //    16: When set various parameters are sanity checked and
 993 //       a message is printed out when they go out of bounds
 994
 995 void FGInitialCondition::Debug(int from)
 996 {
 997   if (debug_lvl <= 0) return;
 998
 999   if (debug_lvl & 1) { // Standard console startup message output
1000   }
1001   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
1002     if (from == 0) cout << "Instantiated: FGInitialCondition" << endl;
1003     if (from == 1) cout << "Destroyed:    FGInitialCondition" << endl;
1004   }
1005   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
1006   }
1007   if (debug_lvl & 8 ) { // Runtime state variables
1008   }
1009   if (debug_lvl & 16) { // Sanity checking
1010   }
1011   if (debug_lvl & 64) {
1012     if (from == 0) { // Constructor
1013       cout << IdSrc << endl;
1014       cout << IdHdr << endl;
1015     }
1016   }
1017 }
1018