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