]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/JSBSim/FGInitialCondition.cpp
Latest JSBSim changes, including some gear tweaking from Jon and some
[flightgear.git] / src / FDM / JSBSim / FGInitialCondition.cpp
index 284cddb23da06323fe8b130c0c9c8a9cc80844dc..29dc1f10f6329f3650d3d34e77b09343c7c7c19b 100644 (file)
@@ -53,25 +53,36 @@ INCLUDES
 #include "FGPosition.h"
 #include "FGAuxiliary.h"
 #include "FGOutput.h"
-#include "FGDefs.h"
+#include "FGConfigFile.h"
 
-static const char *IdSrc = "$Header$";
+static const char *IdSrc = "$Id$";
 static const char *IdHdr = ID_INITIALCONDITION;
 
+//******************************************************************************
 
-FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec) {
-  vt=vc=ve=0;
+FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec)
+{
+  vt=vc=ve=vg=0;
   mach=0;
   alpha=beta=gamma=0;
   theta=phi=psi=0;
   altitude=hdot=0;
   latitude=longitude=0;
-  u=v=w=0;  
+  u=v=w=0;
+  uw=vw=ww=0;
   vnorth=veast=vdown=0;
+  wnorth=weast=wdown=0;
+  whead=wcross=0;
+  wdir=wmag=0;
   lastSpeedSet=setvt;
-  sea_level_radius = EARTHRAD;
-  radius_to_vehicle = EARTHRAD;
+  lastWindSet=setwned;
+  sea_level_radius = FDMExec->GetInertial()->RefRadius();
+  radius_to_vehicle = FDMExec->GetInertial()->RefRadius();
   terrain_altitude = 0;
+
+  salpha=sbeta=stheta=sphi=spsi=sgamma=0;
+  calpha=cbeta=ctheta=cphi=cpsi=cgamma=1;
+
   if(FDMExec != NULL ) {
     fdmex=FDMExec;
     fdmex->GetPosition()->Seth(altitude);
@@ -80,21 +91,27 @@ FGInitialCondition::FGInitialCondition(FGFDMExec *FDMExec) {
     cout << "FGInitialCondition: This class requires a pointer to a valid FGFDMExec object" << endl;
   }
 
+  Debug(0);
 }
 
+//******************************************************************************
 
-FGInitialCondition::~FGInitialCondition(void) {}
+FGInitialCondition::~FGInitialCondition()
+{
+  Debug(1);
+}
 
+//******************************************************************************
 
-void FGInitialCondition::SetVcalibratedKtsIC(float tt) {
+void FGInitialCondition::SetVcalibratedKtsIC(double tt) {
 
-  if(getMachFromVcas(&mach,tt*jsbKTSTOFPS)) {
+  if(getMachFromVcas(&mach,tt*ktstofps)) {
     //cout << "Mach: " << mach << endl;
     lastSpeedSet=setvc;
-    vc=tt*jsbKTSTOFPS;
+    vc=tt*ktstofps;
     vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
     ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
-    //cout << "Vt: " << vt*jsbFPSTOKTS << " Vc: " << vc*jsbFPSTOKTS << endl;
+    //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
   }
   else {
     cout << "Failed to get Mach number for given Vc and altitude, Vc unchanged." << endl;
@@ -102,161 +119,388 @@ void FGInitialCondition::SetVcalibratedKtsIC(float tt) {
   }
 }
 
-void FGInitialCondition::SetVequivalentKtsIC(float tt) {
-  ve=tt*jsbKTSTOFPS;
+//******************************************************************************
+
+void FGInitialCondition::SetVequivalentKtsIC(double tt) {
+  ve=tt*ktstofps;
   lastSpeedSet=setve;
   vt=ve*1/sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
   vc=calcVcas(mach);
 }
 
-void FGInitialCondition::SetVtrueKtsIC(float tt) {
-  vt=tt*jsbKTSTOFPS;
+//******************************************************************************
+
+void FGInitialCondition::SetVgroundFpsIC(double tt) {
+  double ua,va,wa;
+  double vxz;
+
+  vg=tt;
+  lastSpeedSet=setvg;
+  vnorth = vg*cos(psi); veast = vg*sin(psi); vdown = 0;
+  calcUVWfromNED();
+  ua = u + uw; va = v + vw; wa = w + ww;
+  vt = sqrt( ua*ua + va*va + wa*wa );
+  alpha = beta = 0;
+  vxz = sqrt( u*u + w*w );
+  if( w != 0 ) alpha = atan2( w, u );
+  if( vxz != 0 ) beta = atan2( v, vxz );
+  mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
+  vc=calcVcas(mach);
+  ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetVtrueFpsIC(double tt) {
+  vt=tt;
   lastSpeedSet=setvt;
   mach=vt/fdmex->GetAtmosphere()->GetSoundSpeed();
   vc=calcVcas(mach);
   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
 }
 
-void FGInitialCondition::SetMachIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetMachIC(double tt) {
   mach=tt;
   lastSpeedSet=setmach;
   vt=mach*fdmex->GetAtmosphere()->GetSoundSpeed();
   vc=calcVcas(mach);
   ve=vt*sqrt(fdmex->GetAtmosphere()->GetDensityRatio());
-  //cout << "Vt: " << vt*jsbFPSTOKTS << " Vc: " << vc*jsbFPSTOKTS << endl;
+  //cout << "Vt: " << vt*fpstokts << " Vc: " << vc*fpstokts << endl;
 }
 
-void FGInitialCondition::SetClimbRateFpmIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetClimbRateFpmIC(double tt) {
   SetClimbRateFpsIC(tt/60.0);
 }
 
-void FGInitialCondition::SetClimbRateFpsIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetClimbRateFpsIC(double tt) {
 
   if(vt > 0.1) {
     hdot=tt;
     gamma=asin(hdot/vt);
+    sgamma=sin(gamma); cgamma=cos(gamma);
   }
 }
 
-void FGInitialCondition::SetFlightPathAngleRadIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetFlightPathAngleRadIC(double tt) {
   gamma=tt;
+  sgamma=sin(gamma); cgamma=cos(gamma);
   getTheta();
-  hdot=vt*sin(tt);
+  hdot=vt*sgamma;
 }
 
+//******************************************************************************
+
+void FGInitialCondition::SetAlphaRadIC(double tt) {
+  alpha=tt;
+  salpha=sin(alpha); calpha=cos(alpha);
+  getTheta();
+}
 
-void FGInitialCondition::SetUBodyFpsIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetPitchAngleRadIC(double tt) {
+  theta=tt;
+  stheta=sin(theta); ctheta=cos(theta);
+  getAlpha();
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetBetaRadIC(double tt) {
+  beta=tt;
+  sbeta=sin(beta); cbeta=cos(beta);
+  getTheta();
+  
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetRollAngleRadIC(double tt) {
+  phi=tt;
+  sphi=sin(phi); cphi=cos(phi);
+  getTheta();
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetTrueHeadingRadIC(double tt) {
+    psi=tt;
+    spsi=sin(psi); cpsi=cos(psi);
+    calcWindUVW();
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetUBodyFpsIC(double tt) {
   u=tt;
-  vt=sqrt(u*u+v*v+w*w);
+  vt=sqrt(u*u + v*v + w*w);
   lastSpeedSet=setuvw;
 }
 
-  
-void FGInitialCondition::SetVBodyFpsIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetVBodyFpsIC(double tt) {
   v=tt;
-  vt=sqrt(u*u+v*v+w*w);
+  vt=sqrt(u*u + v*v + w*w);
   lastSpeedSet=setuvw;
 }
 
-void FGInitialCondition::SetWBodyFpsIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetWBodyFpsIC(double tt) {
   w=tt;
-  vt=sqrt(u*u+v*v+w*w);
+  vt=sqrt( u*u + v*v + w*w );
   lastSpeedSet=setuvw;
 }
 
+//******************************************************************************
+
+double FGInitialCondition::GetUBodyFpsIC(void) {
+    if(lastSpeedSet == setvg )
+      return u;
+    else
+      return vt*calpha*cbeta - uw;
+}
+
+//******************************************************************************
+
+double FGInitialCondition::GetVBodyFpsIC(void) {
+    if( lastSpeedSet == setvg )
+      return v;
+    else {
+      return vt*sbeta - vw;
+    }  
+}
+
+//******************************************************************************
+
+double FGInitialCondition::GetWBodyFpsIC(void) {
+    if( lastSpeedSet == setvg )
+      return w;
+    else 
+      return vt*salpha*cbeta -ww;
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetWindNEDFpsIC(double wN, double wE, double wD ) {
+  wnorth = wN; weast = wE; wdown = wD;
+  lastWindSet = setwned;
+  calcWindUVW();
+  if(lastSpeedSet == setvg)
+    SetVgroundFpsIC(vg);
+}
+
+//******************************************************************************
+
+// positive from left
+void FGInitialCondition::SetHeadWindKtsIC(double head){ 
+    whead=head*ktstofps;
+    lastWindSet=setwhc; 
+    calcWindUVW();
+    if(lastSpeedSet == setvg)
+      SetVgroundFpsIC(vg);
+
+} 
+
+//******************************************************************************
+
+void FGInitialCondition::SetCrossWindKtsIC(double cross){ 
+    wcross=cross*ktstofps; 
+    lastWindSet=setwhc; 
+    calcWindUVW();
+    if(lastSpeedSet == setvg)
+      SetVgroundFpsIC(vg);
+
+} 
+
+//******************************************************************************
+
+void FGInitialCondition::SetWindDownKtsIC(double wD) { 
+    wdown=wD; 
+    calcWindUVW();
+    if(lastSpeedSet == setvg)
+      SetVgroundFpsIC(vg);
+} 
+
+//******************************************************************************
+
+void FGInitialCondition::SetWindMagKtsIC(double mag) {
+  wmag=mag*ktstofps;
+  lastWindSet=setwmd;
+  calcWindUVW();    
+  if(lastSpeedSet == setvg)
+      SetVgroundFpsIC(vg);
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetWindDirDegIC(double dir) {
+  wdir=dir*degtorad;
+  lastWindSet=setwmd;
+  calcWindUVW();    
+  if(lastSpeedSet == setvg)
+      SetVgroundFpsIC(vg);
+}
+
+
+//******************************************************************************
+
+void FGInitialCondition::calcWindUVW(void) {
+    
+    switch(lastWindSet) {
+      case setwmd:
+        wnorth=wmag*cos(wdir);
+        weast=wmag*sin(wdir);
+      break;
+      case setwhc:
+        wnorth=whead*cos(psi) + wcross*cos(psi+M_PI/2);
+        weast=whead*sin(psi) + wcross*sin(psi+M_PI/2);
+      break;
+      case setwned:
+      break;
+    }    
+    uw=wnorth*ctheta*cpsi +
+       weast*ctheta*spsi -
+       wdown*stheta;
+    vw=wnorth*( sphi*stheta*cpsi - cphi*spsi ) +
+        weast*( sphi*stheta*spsi + cphi*cpsi ) +
+       wdown*sphi*ctheta;
+    ww=wnorth*(cphi*stheta*cpsi + sphi*spsi) +
+       weast*(cphi*stheta*spsi - sphi*cpsi) +
+       wdown*cphi*ctheta;
+            
+   
+    /* cout << "FGInitialCondition::calcWindUVW: wnorth, weast, wdown "
+         << wnorth << ", " << weast << ", " << wdown << endl;
+    cout << "FGInitialCondition::calcWindUVW: theta, phi, psi "
+          << theta << ", " << phi << ", " << psi << endl;
+    cout << "FGInitialCondition::calcWindUVW: uw, vw, ww "
+          << uw << ", " << vw << ", " << ww << endl;   */
+
+}
+
+//******************************************************************************
 
-void FGInitialCondition::SetAltitudeFtIC(float tt) {
+void FGInitialCondition::SetAltitudeFtIC(double tt) {
   altitude=tt;
   fdmex->GetPosition()->Seth(altitude);
   fdmex->GetAtmosphere()->Run();
-
   //lets try to make sure the user gets what they intended
 
   switch(lastSpeedSet) {
   case setned:
   case setuvw:
   case setvt:
-    SetVtrueKtsIC(vt*jsbFPSTOKTS);
+    SetVtrueKtsIC(vt*fpstokts);
     break;
   case setvc:
-    SetVcalibratedKtsIC(vc*jsbFPSTOKTS);
+    SetVcalibratedKtsIC(vc*fpstokts);
     break;
   case setve:
-    SetVequivalentKtsIC(ve*jsbFPSTOKTS);
+    SetVequivalentKtsIC(ve*fpstokts);
     break;
   case setmach:
     SetMachIC(mach);
     break;
+  case setvg:
+    SetVgroundFpsIC(vg);
+    break;
   }
 }
 
-void FGInitialCondition::SetAltitudeAGLFtIC(float tt) {
+//******************************************************************************
+
+void FGInitialCondition::SetAltitudeAGLFtIC(double tt) {
   fdmex->GetPosition()->SetDistanceAGL(tt);
   altitude=fdmex->GetPosition()->Geth();
   SetAltitudeFtIC(altitude);
-} 
+}
+
+//******************************************************************************
+
 void FGInitialCondition::SetSeaLevelRadiusFtIC(double tt) {
   sea_level_radius = tt;
 }
-   
+
+//******************************************************************************
+
 void FGInitialCondition::SetTerrainAltitudeFtIC(double tt) {
   terrain_altitude=tt;
-  fdmex->GetPosition()->SetDistanceAGL(altitude-terrain_altitude);
-  fdmex->GetPosition()->SetRunwayRadius(sea_level_radius + terrain_altitude);
-}  
+}
+
+//******************************************************************************
 
 void FGInitialCondition::calcUVWfromNED(void) {
-  u=vnorth*cos(theta)*cos(psi) + 
-     veast*cos(theta)*sin(psi) - 
-     vdown*sin(theta);
-  v=vnorth*(sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi)) +
-     veast*(sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi)) +
-     vdown*sin(phi)*cos(theta);
-  w=vnorth*(cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi)) +
-     veast*(cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi)) +
-     vdown*cos(phi)*cos(theta);
-}     
-void FGInitialCondition::SetVnorthFpsIC(float tt) {
+  u=vnorth*ctheta*cpsi +
+     veast*ctheta*spsi -
+     vdown*stheta;
+  v=vnorth*( sphi*stheta*cpsi - cphi*spsi ) +
+     veast*( sphi*stheta*spsi + cphi*cpsi ) +
+     vdown*sphi*ctheta;
+  w=vnorth*( cphi*stheta*cpsi + sphi*spsi ) +
+     veast*( cphi*stheta*spsi - sphi*cpsi ) +
+     vdown*cphi*ctheta;
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetVnorthFpsIC(double tt) {
   vnorth=tt;
   calcUVWfromNED();
   vt=sqrt(u*u + v*v + w*w);
   lastSpeedSet=setned;
-}        
-  
-void FGInitialCondition::SetVeastFpsIC(float tt) {
+}
+
+//******************************************************************************
+
+void FGInitialCondition::SetVeastFpsIC(double tt) {
   veast=tt;
   calcUVWfromNED();
   vt=sqrt(u*u + v*v + w*w);
   lastSpeedSet=setned;
-} 
+}
+
+//******************************************************************************
 
-void FGInitialCondition::SetVdownFpsIC(float tt) {
+void FGInitialCondition::SetVdownFpsIC(double tt) {
   vdown=tt;
   calcUVWfromNED();
   vt=sqrt(u*u + v*v + w*w);
   SetClimbRateFpsIC(-1*vdown);
   lastSpeedSet=setned;
-} 
+}
+
+//******************************************************************************
+
+bool FGInitialCondition::getMachFromVcas(double *Mach,double vcas) {
 
-bool FGInitialCondition::getMachFromVcas(float *Mach,float vcas) {
   bool result=false;
-  float guess=1.5;
+  double guess=1.5;
   xlo=xhi=0;
   xmin=0;xmax=50;
   sfunc=&FGInitialCondition::calcVcas;
   if(findInterval(vcas,guess)) {
     if(solve(&mach,vcas))
       result=true;
-  }    
+  }
   return result;
 }
 
+//******************************************************************************
+
 bool FGInitialCondition::getAlpha(void) {
   bool result=false;
-  float guess=theta-gamma;
+  double guess=theta-gamma;
   xlo=xhi=0;
   xmin=fdmex->GetAircraft()->GetAlphaCLMin();
   xmax=fdmex->GetAircraft()->GetAlphaCLMax();
@@ -264,52 +508,70 @@ bool FGInitialCondition::getAlpha(void) {
   if(findInterval(0,guess)){
     if(solve(&alpha,0)){
       result=true;
+      salpha=sin(alpha);
+      calpha=cos(alpha);
     }
   }
+  calcWindUVW();
   return result;
-}      
-    
+}
+
+//******************************************************************************
+
 bool FGInitialCondition::getTheta(void) {
   bool result=false;
-  float guess=alpha+gamma;
+  double guess=alpha+gamma;
   xlo=xhi=0;
   xmin=-89;xmax=89;
   sfunc=&FGInitialCondition::GammaEqOfTheta;
   if(findInterval(0,guess)){
     if(solve(&theta,0)){
       result=true;
+      stheta=sin(theta);
+      ctheta=cos(theta);
     }
   }
+  calcWindUVW();
   return result;
-}      
-  
+}
 
+//******************************************************************************
 
-float FGInitialCondition::GammaEqOfTheta(float Theta) {
-  float a,b,c;
-  
-  a=cos(alpha)*cos(beta)*sin(Theta);
-  b=sin(beta)*sin(phi);
-  c=sin(alpha)*cos(beta)*cos(phi);
-  return sin(gamma)-a+(b+c)*cos(Theta);
+double FGInitialCondition::GammaEqOfTheta(double Theta) {
+  double a,b,c;
+  double sTheta,cTheta;
+
+  //theta=Theta; stheta=sin(theta); ctheta=cos(theta);
+  sTheta=sin(Theta); cTheta=cos(Theta);
+  calcWindUVW();
+  a=wdown + vt*calpha*cbeta + uw;
+  b=vt*sphi*sbeta + vw*sphi;
+  c=vt*cphi*salpha*cbeta + ww*cphi;
+  return vt*sgamma - ( a*sTheta - (b+c)*cTheta);
 }
 
-float FGInitialCondition::GammaEqOfAlpha(float Alpha) {
-  float a,b,c;
-  
-  a=cos(Alpha)*cos(beta)*sin(theta);
-  b=sin(beta)*sin(phi);
-  c=sin(Alpha)*cos(beta)*cos(phi);
-  return sin(gamma)-a+(b+c)*cos(theta);
+//******************************************************************************
+
+double FGInitialCondition::GammaEqOfAlpha(double Alpha) {
+  double a,b,c;
+  double sAlpha,cAlpha;
+
+  sAlpha=sin(Alpha); cAlpha=cos(Alpha);
+  a=wdown + vt*cAlpha*cbeta + uw;
+  b=vt*sphi*sbeta + vw*sphi;
+  c=vt*cphi*sAlpha*cbeta + ww*cphi;
+
+  return vt*sgamma - ( a*stheta - (b+c)*ctheta );
 }
 
-  
-float FGInitialCondition::calcVcas(float Mach) {
+//******************************************************************************
+
+double FGInitialCondition::calcVcas(double Mach) {
 
-  float p=fdmex->GetAtmosphere()->GetPressure();
-  float psl=fdmex->GetAtmosphere()->GetPressureSL();
-  float rhosl=fdmex->GetAtmosphere()->GetDensitySL();
-  float pt,A,B,D,vcas;
+  double p=fdmex->GetAtmosphere()->GetPressure();
+  double psl=fdmex->GetAtmosphere()->GetPressureSL();
+  double rhosl=fdmex->GetAtmosphere()->GetDensitySL();
+  double pt,A,B,D,vcas;
   if(Mach < 0) Mach=0;
   if(Mach < 1)    //calculate total pressure assuming isentropic flow
     pt=p*pow((1 + 0.2*Mach*Mach),3.5);
@@ -339,17 +601,19 @@ float FGInitialCondition::calcVcas(float Mach) {
 
   A = pow(((pt-p)/psl+1),0.28571);
   vcas = sqrt(7*psl/rhosl*(A-1));
-  //cout << "calcVcas: vcas= " << vcas*jsbFPSTOKTS << " mach= " << Mach << " pressure: " << pt << endl;
+  //cout << "calcVcas: vcas= " << vcas*fpstokts << " mach= " << Mach << " pressure: " << pt << endl;
   return vcas;
 }
 
-bool FGInitialCondition::findInterval(float x,float guess) {
-  //void find_interval(inter_params &ip,eqfunc f,float y,float constant, int &flag){
+//******************************************************************************
+
+bool FGInitialCondition::findInterval(double x,double guess) {
+  //void find_interval(inter_params &ip,eqfunc f,double y,double constant, int &flag){
 
   int i=0;
   bool found=false;
-  float flo,fhi,fguess;
-  float lo,hi,step;
+  double flo,fhi,fguess;
+  double lo,hi,step;
   step=0.1;
   fguess=(this->*sfunc)(guess)-x;
   lo=hi=guess;
@@ -380,53 +644,168 @@ bool FGInitialCondition::findInterval(float x,float guess) {
   return found;
 }
 
+//******************************************************************************
 
-
-
-bool FGInitialCondition::solve(float *y,float x) {  
-  float x1,x2,x3,f1,f2,f3,d,d0;
-  float eps=1E-5;
-  float const relax =0.9;
+bool FGInitialCondition::solve(double *y,double x)
+{
+  double x1,x2,x3,f1,f2,f3,d,d0;
+  double eps=1E-5;
+  double const relax =0.9;
   int i;
   bool success=false;
-  
-   //initializations
+
+  //initializations
   d=1;
-  
-    x1=xlo;x3=xhi;
-    f1=(this->*sfunc)(x1)-x;
-    f3=(this->*sfunc)(x3)-x;
-    d0=fabs(x3-x1);
-  
-    //iterations
-    i=0;
-    while ((fabs(d) > eps) && (i < 100)) {
-      d=(x3-x1)/d0;
-      x2=x1-d*d0*f1/(f3-f1);
-      
-      f2=(this->*sfunc)(x2)-x;
-      //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
-      //cout << "                " << f1 << "," << f2 << "," << f3 << endl;
-
-      if(fabs(f2) <= 0.001) {
-        x1=x3=x2;
-      } else if(f1*f2 <= 0.0) {
-        x3=x2;
-        f3=f2;
-        f1=relax*f1;
-      } else if(f2*f3 <= 0) {
-        x1=x2;
-        f1=f2;
-        f3=relax*f3;
-      }
-      //cout << i << endl;
-      i++;
-    }//end while
-    if(i < 100) {
-      success=true;
-      *y=x2;
+  x2 = 0;
+  x1=xlo;x3=xhi;
+  f1=(this->*sfunc)(x1)-x;
+  f3=(this->*sfunc)(x3)-x;
+  d0=fabs(x3-x1);
+
+  //iterations
+  i=0;
+  while ((fabs(d) > eps) && (i < 100)) {
+    d=(x3-x1)/d0;
+    x2 = x1-d*d0*f1/(f3-f1);
+    
+    f2=(this->*sfunc)(x2)-x;
+    //cout << "solve x1,x2,x3: " << x1 << "," << x2 << "," << x3 << endl;
+    //cout << "                " << f1 << "," << f2 << "," << f3 << endl;
+
+    if(fabs(f2) <= 0.001) {
+      x1=x3=x2;
+    } else if(f1*f2 <= 0.0) {
+      x3=x2;
+      f3=f2;
+      f1=relax*f1;
+    } else if(f2*f3 <= 0) {
+      x1=x2;
+      f1=f2;
+      f3=relax*f3;
     }
+    //cout << i << endl;
+    i++;
+  }//end while
+  if(i < 100) {
+    success=true;
+    *y=x2;
+  }
 
-  //cout << "Success= " << success << " Vcas: " << vcas*jsbFPSTOKTS << " Mach: " << x2 << endl;
+  //cout << "Success= " << success << " Vcas: " << vcas*fpstokts << " Mach: " << x2 << endl;
   return success;
 }
+
+//******************************************************************************
+
+double FGInitialCondition::GetWindDirDegIC(void) {
+  if(weast != 0.0) 
+    return atan2(weast,wnorth)*radtodeg;
+  else if(wnorth > 0) 
+    return 0.0;
+  else
+    return 180.0;
+}        
+
+//******************************************************************************
+
+bool FGInitialCondition::Load(string acpath, string acname, string rstfile)
+{
+  string resetDef;
+  string token="";
+
+  double temp;
+
+# ifndef macintosh
+  resetDef = acpath + "/" + acname + "/" + rstfile + ".xml";
+# else
+  resetDef = acpath + ";" + acname + ";" + rstfile + ".xml";
+# endif
+
+  FGConfigFile resetfile(resetDef);
+  if (!resetfile.IsOpen()) {
+    cerr << "Failed to open reset file: " << resetDef << endl;
+    return false;
+  }  
+
+  resetfile.GetNextConfigLine();
+  token = resetfile.GetValue();
+  if (token != string("initialize")) {
+    cerr << "The reset file " << resetDef
+         << " does not appear to be a reset file" << endl;
+    return false;
+  }
+  
+  resetfile.GetNextConfigLine();
+  resetfile >> token;
+  while (token != string("/initialize") && token != string("EOF")) {
+    if (token == "UBODY" ) { resetfile >> temp; SetUBodyFpsIC(temp); } 
+    if (token == "VBODY" ) { resetfile >> temp; SetVBodyFpsIC(temp); } 
+    if (token == "WBODY" ) { resetfile >> temp; SetWBodyFpsIC(temp); }  
+    if (token == "LATITUDE" ) { resetfile >> temp; SetLatitudeDegIC(temp); }
+    if (token == "LONGITUDE" ) { resetfile >> temp; SetLongitudeDegIC(temp); }
+    if (token == "PHI" ) { resetfile >> temp; SetRollAngleDegIC(temp); }
+    if (token == "THETA" ) { resetfile >> temp; SetPitchAngleDegIC(temp); }
+    if (token == "PSI" ) { resetfile >> temp; SetTrueHeadingDegIC(temp); }
+    if (token == "ALPHA" ) { resetfile >> temp; SetAlphaDegIC(temp); }
+    if (token == "BETA" ) { resetfile >> temp; SetBetaDegIC(temp); }
+    if (token == "GAMMA" ) { resetfile >> temp; SetFlightPathAngleDegIC(temp); }
+    if (token == "ROC" ) { resetfile >> temp; SetClimbRateFpmIC(temp); }
+    if (token == "ALTITUDE" ) { resetfile >> temp; SetAltitudeFtIC(temp); }
+    if (token == "WINDDIR" ) { resetfile >> temp; SetWindDirDegIC(temp); }
+    if (token == "VWIND" ) { resetfile >> temp; SetWindMagKtsIC(temp); }
+    if (token == "HWIND" ) { resetfile >> temp; SetHeadWindKtsIC(temp); }
+    if (token == "XWIND" ) { resetfile >> temp; SetCrossWindKtsIC(temp); }
+    if (token == "VC" ) { resetfile >> temp; SetVcalibratedKtsIC(temp); }
+    if (token == "MACH" ) { resetfile >> temp; SetMachIC(temp); }
+    if (token == "VGROUND" ) { resetfile >> temp; SetVgroundKtsIC(temp); }
+    resetfile >> token;
+  }
+
+  fdmex->RunIC(this);
+  
+  return true;
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//    The bitmasked value choices are as follows:
+//    unset: In this case (the default) JSBSim would only print
+//       out the normally expected messages, essentially echoing
+//       the config files as they are read. If the environment
+//       variable is not set, debug_lvl is set to 1 internally
+//    0: This requests JSBSim not to output any messages
+//       whatsoever.
+//    1: This value explicity requests the normal JSBSim
+//       startup messages
+//    2: This value asks for a message to be printed out when
+//       a class is instantiated
+//    4: When this value is set, a message is displayed when a
+//       FGModel object executes its Run() method
+//    8: When this value is set, various runtime state variables
+//       are printed out periodically
+//    16: When set various parameters are sanity checked and
+//       a message is printed out when they go out of bounds
+
+void FGInitialCondition::Debug(int from)
+{
+  if (debug_lvl <= 0) return;
+
+  if (debug_lvl & 1) { // Standard console startup message output
+  }
+  if (debug_lvl & 2 ) { // Instantiation/Destruction notification
+    if (from == 0) cout << "Instantiated: FGInitialCondition" << endl;
+    if (from == 1) cout << "Destroyed:    FGInitialCondition" << endl;
+  }
+  if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
+  }
+  if (debug_lvl & 8 ) { // Runtime state variables
+  }
+  if (debug_lvl & 16) { // Sanity checking
+  }
+  if (debug_lvl & 64) {
+    if (from == 0) { // Constructor
+      cout << IdSrc << endl;
+      cout << IdHdr << endl;
+    }
+  }
+}
+