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