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->GetAerodynamics()->GetAlphaCLMin();
 513   xmax=fdmex->GetAerodynamics()->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
 780 void FGInitialCondition::bind(void){
 781   PropertyManager->Tie("ic/vc-kts", this,
 782                        &FGInitialCondition::GetVcalibratedKtsIC,
 783                        &FGInitialCondition::SetVcalibratedKtsIC,
 784                        true);
 785   PropertyManager->Tie("ic/ve-kts", this,
 786                        &FGInitialCondition::GetVequivalentKtsIC,
 787                        &FGInitialCondition::SetVequivalentKtsIC,
 788                        true);
 789   PropertyManager->Tie("ic/vg-kts", this,
 790                        &FGInitialCondition::GetVgroundKtsIC,
 791                        &FGInitialCondition::SetVgroundKtsIC,
 792                        true);
 793   PropertyManager->Tie("ic/vt-kts", this,
 794                        &FGInitialCondition::GetVtrueKtsIC,
 795                        &FGInitialCondition::SetVtrueKtsIC,
 796                        true);
 797   PropertyManager->Tie("ic/mach-norm", this,
 798                        &FGInitialCondition::GetMachIC,
 799                        &FGInitialCondition::SetMachIC,
 800                        true);
 801   PropertyManager->Tie("ic/roc-fpm", this,
 802                        &FGInitialCondition::GetClimbRateFpmIC,
 803                        &FGInitialCondition::SetClimbRateFpmIC,
 804                        true);
 805   PropertyManager->Tie("ic/gamma-deg", this,
 806                        &FGInitialCondition::GetFlightPathAngleDegIC,
 807                        &FGInitialCondition::SetFlightPathAngleDegIC,
 808                        true);
 809   PropertyManager->Tie("ic/alpha-deg", this,
 810                        &FGInitialCondition::GetAlphaDegIC,
 811                        &FGInitialCondition::SetAlphaDegIC,
 812                        true);
 813   PropertyManager->Tie("ic/beta-deg", this,
 814                        &FGInitialCondition::GetBetaDegIC,
 815                        &FGInitialCondition::SetBetaDegIC,
 816                        true);
 817   PropertyManager->Tie("ic/theta-deg", this,
 818                        &FGInitialCondition::GetPitchAngleDegIC,
 819                        &FGInitialCondition::SetPitchAngleDegIC,
 820                        true);
 821   PropertyManager->Tie("ic/phi-deg", this,
 822                        &FGInitialCondition::GetRollAngleDegIC,
 823                        &FGInitialCondition::SetRollAngleDegIC,
 824                        true);
 825   PropertyManager->Tie("ic/psi-true-deg", this,
 826                        &FGInitialCondition::GetHeadingDegIC );
 827   PropertyManager->Tie("ic/lat-gc-deg", this,
 828                        &FGInitialCondition::GetLatitudeDegIC,
 829                        &FGInitialCondition::SetLatitudeDegIC,
 830                        true);
 831   PropertyManager->Tie("ic/long-gc-deg", this,
 832                        &FGInitialCondition::GetLongitudeDegIC,
 833                        &FGInitialCondition::SetLongitudeDegIC,
 834                        true);
 835   PropertyManager->Tie("ic/h-sl-ft", this,
 836                        &FGInitialCondition::GetAltitudeFtIC,
 837                        &FGInitialCondition::SetAltitudeFtIC,
 838                        true);
 839   PropertyManager->Tie("ic/h-agl-ft", this,
 840                        &FGInitialCondition::GetAltitudeAGLFtIC,
 841                        &FGInitialCondition::SetAltitudeAGLFtIC,
 842                        true);
 843   PropertyManager->Tie("ic/sea-level-radius-ft", this,
 844                        &FGInitialCondition::GetSeaLevelRadiusFtIC,
 845                        &FGInitialCondition::SetSeaLevelRadiusFtIC,
 846                        true);
 847   PropertyManager->Tie("ic/terrain-altitude-ft", this,
 848                        &FGInitialCondition::GetTerrainAltitudeFtIC,
 849                        &FGInitialCondition::SetTerrainAltitudeFtIC,
 850                        true);
 851   PropertyManager->Tie("ic/vg-fps", this,
 852                        &FGInitialCondition::GetVgroundFpsIC,
 853                        &FGInitialCondition::SetVgroundFpsIC,
 854                        true);
 855   PropertyManager->Tie("ic/vt-fps", this,
 856                        &FGInitialCondition::GetVtrueFpsIC,
 857                        &FGInitialCondition::SetVtrueFpsIC,
 858                        true);
 859   PropertyManager->Tie("ic/vw-bx-fps", this,
 860                        &FGInitialCondition::GetWindUFpsIC);
 861   PropertyManager->Tie("ic/vw-by-fps", this,
 862                        &FGInitialCondition::GetWindVFpsIC);
 863   PropertyManager->Tie("ic/vw-bz-fps", this,
 864                        &FGInitialCondition::GetWindWFpsIC);
 865   PropertyManager->Tie("ic/vw-north-fps", this,
 866                        &FGInitialCondition::GetWindNFpsIC);
 867   PropertyManager->Tie("ic/vw-east-fps", this,
 868                        &FGInitialCondition::GetWindEFpsIC);
 869   PropertyManager->Tie("ic/vw-down-fps", this,
 870                        &FGInitialCondition::GetWindDFpsIC);
 871   PropertyManager->Tie("ic/vw-mag-fps", this,
 872                        &FGInitialCondition::GetWindFpsIC);
 873  /*  PropertyManager->Tie("ic/vw-dir-deg", this,
 874                        &FGInitialCondition::GetWindDirDegIC,
 875                        &FGInitialCondition::SetWindDirDegIC,
 876                        true); */
 877
 878   PropertyManager->Tie("ic/roc-fps", this,
 879                        &FGInitialCondition::GetClimbRateFpsIC,
 880                        &FGInitialCondition::SetClimbRateFpsIC,
 881                        true);
 882   /* PropertyManager->Tie("ic/u-fps", this,
 883                        &FGInitialCondition::GetUBodyFpsIC,
 884                        &FGInitialCondition::SetUBodyFpsIC,
 885                        true);
 886   PropertyManager->Tie("ic/v-fps", this,
 887                        &FGInitialCondition::GetVBodyFpsIC,
 888                        &FGInitialCondition::SetVBodyFpsIC,
 889                        true);
 890   PropertyManager->Tie("ic/w-fps", this,
 891                        &FGInitialCondition::GetWBodyFpsIC,
 892                        &FGInitialCondition::SetWBodyFpsIC,
 893                        true); */
 894
 895   PropertyManager->Tie("ic/gamma-rad", this,
 896                        &FGInitialCondition::GetFlightPathAngleRadIC,
 897                        &FGInitialCondition::SetFlightPathAngleRadIC,
 898                        true);
 899   PropertyManager->Tie("ic/alpha-rad", this,
 900                        &FGInitialCondition::GetAlphaRadIC,
 901                        &FGInitialCondition::SetAlphaRadIC,
 902                        true);
 903   PropertyManager->Tie("ic/theta-rad", this,
 904                        &FGInitialCondition::GetPitchAngleRadIC,
 905                        &FGInitialCondition::SetPitchAngleRadIC,
 906                        true);
 907   PropertyManager->Tie("ic/beta-rad", this,
 908                        &FGInitialCondition::GetBetaRadIC,
 909                        &FGInitialCondition::SetBetaRadIC,
 910                        true);
 911   PropertyManager->Tie("ic/phi-rad", this,
 912                        &FGInitialCondition::GetRollAngleRadIC,
 913                        &FGInitialCondition::SetRollAngleRadIC,
 914                        true);
 915   PropertyManager->Tie("ic/psi-true-rad", this,
 916                        &FGInitialCondition::GetHeadingRadIC);
 917   PropertyManager->Tie("ic/lat-gc-rad", this,
 918                        &FGInitialCondition::GetLatitudeRadIC,
 919                        &FGInitialCondition::SetLatitudeRadIC,
 920                        true);
 921   PropertyManager->Tie("ic/long-gc-rad", this,
 922                        &FGInitialCondition::GetLongitudeRadIC,
 923                        &FGInitialCondition::SetLongitudeRadIC,
 924                        true);
 925 }
 926
 927 //******************************************************************************
 928
 929 void FGInitialCondition::unbind(void){
 930   PropertyManager->Untie("ic/vc-kts");
 931   PropertyManager->Untie("ic/ve-kts");
 932   PropertyManager->Untie("ic/vg-kts");
 933   PropertyManager->Untie("ic/vt-kts");
 934   PropertyManager->Untie("ic/mach-norm");
 935   PropertyManager->Untie("ic/roc-fpm");
 936   PropertyManager->Untie("ic/gamma-deg");
 937   PropertyManager->Untie("ic/alpha-deg");
 938   PropertyManager->Untie("ic/beta-deg");
 939   PropertyManager->Untie("ic/theta-deg");
 940   PropertyManager->Untie("ic/phi-deg");
 941   PropertyManager->Untie("ic/psi-true-deg");
 942   PropertyManager->Untie("ic/lat-gc-deg");
 943   PropertyManager->Untie("ic/long-gc-deg");
 944   PropertyManager->Untie("ic/h-sl-ft");
 945   PropertyManager->Untie("ic/h-agl-ft");
 946   PropertyManager->Untie("ic/sea-level-radius-ft");
 947   PropertyManager->Untie("ic/terrain-altitude-ft");
 948   PropertyManager->Untie("ic/vg-fps");
 949   PropertyManager->Untie("ic/vt-fps");
 950   PropertyManager->Untie("ic/vw-bx-fps");
 951   PropertyManager->Untie("ic/vw-by-fps");
 952   PropertyManager->Untie("ic/vw-bz-fps");
 953   PropertyManager->Untie("ic/vw-north-fps");
 954   PropertyManager->Untie("ic/vw-east-fps");
 955   PropertyManager->Untie("ic/vw-down-fps");
 956   PropertyManager->Untie("ic/vw-mag-fps");
 957   /* PropertyManager->Untie("ic/vw-dir-deg"); */
 958
 959   PropertyManager->Untie("ic/roc-fps");
 960
 961   /*  PropertyManager->Untie("ic/u-fps");
 962   PropertyManager->Untie("ic/v-fps");
 963   PropertyManager->Untie("ic/w-fps"); */
 964
 965   PropertyManager->Untie("ic/gamma-rad");
 966   PropertyManager->Untie("ic/alpha-rad");
 967   PropertyManager->Untie("ic/theta-rad");
 968   PropertyManager->Untie("ic/beta-rad");
 969   PropertyManager->Untie("ic/phi-rad");
 970   PropertyManager->Untie("ic/psi-true-rad");
 971   PropertyManager->Untie("ic/lat-gc-rad");
 972   PropertyManager->Untie("ic/long-gc-rad");
 973 }
 974
 975 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 976 //    The bitmasked value choices are as follows:
 977 //    unset: In this case (the default) JSBSim would only print
 978 //       out the normally expected messages, essentially echoing
 979 //       the config files as they are read. If the environment
 980 //       variable is not set, debug_lvl is set to 1 internally
 981 //    0: This requests JSBSim not to output any messages
 982 //       whatsoever.
 983 //    1: This value explicity requests the normal JSBSim
 984 //       startup messages
 985 //    2: This value asks for a message to be printed out when
 986 //       a class is instantiated
 987 //    4: When this value is set, a message is displayed when a
 988 //       FGModel object executes its Run() method
 989 //    8: When this value is set, various runtime state variables
 990 //       are printed out periodically
 991 //    16: When set various parameters are sanity checked and
 992 //       a message is printed out when they go out of bounds
 993
 994 void FGInitialCondition::Debug(int from)
 995 {
 996   if (debug_lvl <= 0) return;
 997
 998   if (debug_lvl & 1) { // Standard console startup message output
 999   }
1000   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
1001     if (from == 0) cout << "Instantiated: FGInitialCondition" << endl;
1002     if (from == 1) cout << "Destroyed:    FGInitialCondition" << endl;
1003   }
1004   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
1005   }
1006   if (debug_lvl & 8 ) { // Runtime state variables
1007   }
1008   if (debug_lvl & 16) { // Sanity checking
1009   }
1010   if (debug_lvl & 64) {
1011     if (from == 0) { // Constructor
1012       cout << IdSrc << endl;
1013       cout << IdHdr << endl;
1014     }
1015   }
1016 }
1017