]> git.mxchange.org Git - flightgear.git/blobdiff - src/FDM/JSBSim/models/FGAuxiliary.cpp
Merge branch 'next' into durk-atc
[flightgear.git] / src / FDM / JSBSim / models / FGAuxiliary.cpp
index d8a8d77c0f2ce22d40e3356e607865d12b465bff..36c1e084abdf2fb4f99819e0b0e62fb6c6a6f978 100755 (executable)
@@ -4,7 +4,7 @@
  Author:       Tony Peden, Jon Berndt
  Date started: 01/26/99
  Purpose:      Calculates additional parameters needed by the visual system, etc.
- Called by:    FGSimExec
+ Called by:    FGFDMExec
 
  ------------- Copyright (C) 1999  Jon S. Berndt (jon@jsbsim.org) -------------
 
@@ -59,7 +59,7 @@ using namespace std;
 
 namespace JSBSim {
 
-static const char *IdSrc = "$Id$";
+static const char *IdSrc = "$Id: FGAuxiliary.cpp,v 1.47 2011/03/29 11:49:27 jberndt Exp $";
 static const char *IdHdr = ID_AUXILIARY;
 
 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -142,22 +142,31 @@ bool FGAuxiliary::Run()
   if (FGModel::Run()) return true; // return true if error returned from base class
   if (FDMExec->Holding()) return false;
 
-  const FGColumnVector3& vPQR = Propagate->GetPQR();
-  const FGColumnVector3& vUVW = Propagate->GetUVW();
-  const FGColumnVector3& vUVWdot = Propagate->GetUVWdot();
-  const FGColumnVector3& vVel = Propagate->GetVel();
+  RunPreFunctions();
 
-  p = Atmosphere->GetPressure();
-  rhosl = Atmosphere->GetDensitySL();
-  psl = Atmosphere->GetPressureSL();
-  sat = Atmosphere->GetTemperature();
+  const double density = FDMExec->GetAtmosphere()->GetDensity();
+  const double soundspeed = FDMExec->GetAtmosphere()->GetSoundSpeed();
+  const double DistanceAGL = FDMExec->GetPropagate()->GetDistanceAGL();
+  const double wingspan = FDMExec->GetAircraft()->GetWingSpan();
+  const FGMatrix33& Tl2b = FDMExec->GetPropagate()->GetTl2b();
+  const FGMatrix33& Tb2l = FDMExec->GetPropagate()->GetTb2l();
+
+  const FGColumnVector3& vPQR = FDMExec->GetPropagate()->GetPQR();
+  const FGColumnVector3& vUVW = FDMExec->GetPropagate()->GetUVW();
+  const FGColumnVector3& vUVWdot = FDMExec->GetPropagate()->GetUVWdot();
+  const FGColumnVector3& vVel = FDMExec->GetPropagate()->GetVel();
+
+  p = FDMExec->GetAtmosphere()->GetPressure();
+  rhosl = FDMExec->GetAtmosphere()->GetDensitySL();
+  psl = FDMExec->GetAtmosphere()->GetPressureSL();
+  sat = FDMExec->GetAtmosphere()->GetTemperature();
 
 // Rotation
 
-  double cTht = Propagate->GetCosEuler(eTht);
-  double sTht = Propagate->GetSinEuler(eTht);
-  double cPhi = Propagate->GetCosEuler(ePhi);
-  double sPhi = Propagate->GetSinEuler(ePhi);
+  double cTht = FDMExec->GetPropagate()->GetCosEuler(eTht);
+  double sTht = FDMExec->GetPropagate()->GetSinEuler(eTht);
+  double cPhi = FDMExec->GetPropagate()->GetCosEuler(ePhi);
+  double sPhi = FDMExec->GetPropagate()->GetSinEuler(ePhi);
 
   vEulerRates(eTht) = vPQR(eQ)*cPhi - vPQR(eR)*sPhi;
   if (cTht != 0.0) {
@@ -165,57 +174,41 @@ bool FGAuxiliary::Run()
     vEulerRates(ePhi) = vPQR(eP) + vEulerRates(ePsi)*sTht;
   }
 
-// 12/16/2005, JSB: For ground handling purposes, at this time, let's ramp
-// in the effects of wind from 10 fps to 30 fps when there is weight on the
-// landing gear wheels.
-
-  if (GroundReactions->GetWOW() && vUVW(eU) < 10) {
-    vAeroPQR = vPQR;
-    vAeroUVW = vUVW;
-  } else if (GroundReactions->GetWOW() && vUVW(eU) < 30) {
-    double factor = (vUVW(eU) - 10.0)/20.0;
-    vAeroPQR = vPQR - factor*Atmosphere->GetTurbPQR();
-    vAeroUVW = vUVW - factor*Propagate->GetTl2b()*Atmosphere->GetTotalWindNED();
-  } else {
-    FGColumnVector3 wind = Propagate->GetTl2b()*Atmosphere->GetTotalWindNED();
-    vAeroPQR = vPQR - Atmosphere->GetTurbPQR();
-    vAeroUVW = vUVW - wind;
-  }
+// Combine the wind speed with aircraft speed to obtain wind relative speed
+  FGColumnVector3 wind = Tl2b*FDMExec->GetAtmosphere()->GetTotalWindNED();
+  vAeroPQR = vPQR - FDMExec->GetAtmosphere()->GetTurbPQR();
+  vAeroUVW = vUVW - wind;
 
   Vt = vAeroUVW.Magnitude();
-  if ( Vt > 0.05) {
+  double Vt2 = Vt*Vt;
+  alpha = beta = adot = bdot = 0;
+  double mUW = (vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eW)*vAeroUVW(eW));
+
+  if ( Vt > 1.0 ) {
     if (vAeroUVW(eW) != 0.0)
       alpha = vAeroUVW(eU)*vAeroUVW(eU) > 0.0 ? atan2(vAeroUVW(eW), vAeroUVW(eU)) : 0.0;
     if (vAeroUVW(eV) != 0.0)
-      beta = vAeroUVW(eU)*vAeroUVW(eU)+vAeroUVW(eW)*vAeroUVW(eW) > 0.0 ? atan2(vAeroUVW(eV),
-             sqrt(vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eW)*vAeroUVW(eW))) : 0.0;
+      beta = mUW > 0.0 ? atan2(vAeroUVW(eV), sqrt(mUW)) : 0.0;
 
-    double mUW = (vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eW)*vAeroUVW(eW));
     double signU=1;
-    if (vAeroUVW(eU) != 0.0)
-      signU = vAeroUVW(eU)/fabs(vAeroUVW(eU));
+    if (vAeroUVW(eU) < 0.0) signU=-1;
 
-    if ( (mUW == 0.0) || (Vt == 0.0) ) {
-      adot = 0.0;
-      bdot = 0.0;
-    } else {
+    if ( mUW >= 1.0 ) {
       adot = (vAeroUVW(eU)*vUVWdot(eW) - vAeroUVW(eW)*vUVWdot(eU))/mUW;
-      bdot = (signU*mUW*vUVWdot(eV) - vAeroUVW(eV)*(vAeroUVW(eU)*vUVWdot(eU)
-              + vAeroUVW(eW)*vUVWdot(eW)))/(Vt*Vt*sqrt(mUW));
+      bdot = (signU*mUW*vUVWdot(eV)
+             - vAeroUVW(eV)*(vAeroUVW(eU)*vUVWdot(eU) + vAeroUVW(eW)*vUVWdot(eW)))/(Vt2*sqrt(mUW));
     }
-  } else {
-    alpha = beta = adot = bdot = 0;
   }
 
-  Re = Vt * Aircraft->Getcbar() / Atmosphere->GetKinematicViscosity();
+  Re = Vt * FDMExec->GetAircraft()->Getcbar() / FDMExec->GetAtmosphere()->GetKinematicViscosity();
 
-  qbar = 0.5*Atmosphere->GetDensity()*Vt*Vt;
-  qbarUW = 0.5*Atmosphere->GetDensity()*(vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eW)*vAeroUVW(eW));
-  qbarUV = 0.5*Atmosphere->GetDensity()*(vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eV)*vAeroUVW(eV));
-  Mach = Vt / Atmosphere->GetSoundSpeed();
-  MachU = vMachUVW(eU) = vAeroUVW(eU) / Atmosphere->GetSoundSpeed();
-  vMachUVW(eV) = vAeroUVW(eV) / Atmosphere->GetSoundSpeed();
-  vMachUVW(eW) = vAeroUVW(eW) / Atmosphere->GetSoundSpeed();
+  qbar = 0.5*density*Vt2;
+  qbarUW = 0.5*density*(mUW);
+  qbarUV = 0.5*density*(vAeroUVW(eU)*vAeroUVW(eU) + vAeroUVW(eV)*vAeroUVW(eV));
+  Mach = Vt / soundspeed;
+  MachU = vMachUVW(eU) = vAeroUVW(eU) / soundspeed;
+  vMachUVW(eV) = vAeroUVW(eV) / soundspeed;
+  vMachUVW(eW) = vAeroUVW(eW) / soundspeed;
 
 // Position
 
@@ -245,19 +238,15 @@ bool FGAuxiliary::Run()
     vcas = veas = 0.0;
   }
 
+  const double SLgravity = FDMExec->GetInertial()->SLgravity();
+
   vPilotAccel.InitMatrix();
   if ( Vt > 1.0 ) {
-     vAircraftAccel = Aerodynamics->GetForces()
-                    + Propulsion->GetForces()
-                    + GroundReactions->GetForces()
-                    + ExternalReactions->GetForces()
-                    + BuoyantForces->GetForces();
-
-     vAircraftAccel /= MassBalance->GetMass();
+     vAircraftAccel = FDMExec->GetAircraft()->GetBodyAccel();
      // Nz is Acceleration in "g's", along normal axis (-Z body axis)
-     Nz = -vAircraftAccel(eZ)/Inertial->gravity();
-     vToEyePt = MassBalance->StructuralToBody(Aircraft->GetXYZep());
-     vPilotAccel = vAircraftAccel + Propagate->GetPQRdot() * vToEyePt;
+     Nz = -vAircraftAccel(eZ)/SLgravity;
+     vToEyePt = FDMExec->GetMassBalance()->StructuralToBody(FDMExec->GetAircraft()->GetXYZep());
+     vPilotAccel = vAircraftAccel + FDMExec->GetPropagate()->GetPQRdot() * vToEyePt;
      vPilotAccel += vPQR * (vPQR * vToEyePt);
   } else {
      // The line below handles low velocity (and on-ground) cases, basically
@@ -266,30 +255,32 @@ bool FGAuxiliary::Run()
      // any jitter that could be introduced by the landing gear. Theoretically,
      // this branch could be eliminated, with a penalty of having a short
      // transient at startup (lasting only a fraction of a second).
-     vPilotAccel = Propagate->GetTl2b() * FGColumnVector3( 0.0, 0.0, -Inertial->gravity() );
-     Nz = -vPilotAccel(eZ)/Inertial->gravity();
+     vPilotAccel = Tl2b * FGColumnVector3( 0.0, 0.0, -SLgravity );
+     Nz = -vPilotAccel(eZ)/SLgravity;
   }
 
-  vPilotAccelN = vPilotAccel/Inertial->gravity();
+  vPilotAccelN = vPilotAccel/SLgravity;
 
   // VRP computation
-  const FGLocation& vLocation = Propagate->GetLocation();
-  FGColumnVector3 vrpStructural = Aircraft->GetXYZvrp();
-  FGColumnVector3 vrpBody = MassBalance->StructuralToBody( vrpStructural );
-  FGColumnVector3 vrpLocal = Propagate->GetTb2l() * vrpBody;
+  const FGLocation& vLocation = FDMExec->GetPropagate()->GetLocation();
+  const FGColumnVector3& vrpStructural = FDMExec->GetAircraft()->GetXYZvrp();
+  const FGColumnVector3 vrpBody = FDMExec->GetMassBalance()->StructuralToBody( vrpStructural );
+  const FGColumnVector3 vrpLocal = Tb2l * vrpBody;
   vLocationVRP = vLocation.LocalToLocation( vrpLocal );
 
   // Recompute some derived values now that we know the dependent parameters values ...
-  hoverbcg = Propagate->GetDistanceAGL() / Aircraft->GetWingSpan();
+  hoverbcg = DistanceAGL / wingspan;
 
-  FGColumnVector3 vMac = Propagate->GetTb2l()*MassBalance->StructuralToBody(Aircraft->GetXYZrp());
-  hoverbmac = (Propagate->GetDistanceAGL() + vMac(3)) / Aircraft->GetWingSpan();
+  FGColumnVector3 vMac = Tb2l*FDMExec->GetMassBalance()->StructuralToBody(FDMExec->GetAircraft()->GetXYZrp());
+  hoverbmac = (DistanceAGL + vMac(3)) / wingspan;
 
   // when all model are executed, 
   // please calculate the distance from the initial point
 
   CalculateRelativePosition();
 
+  RunPostFunctions();
+
   return false;
 }
 
@@ -297,15 +288,16 @@ bool FGAuxiliary::Run()
 //
 // A positive headwind is blowing with you, a negative headwind is blowing against you.
 // psi is the direction the wind is blowing *towards*.
+// ToDo: should this simply be in the atmosphere class? Same with Get Crosswind.
 
 double FGAuxiliary::GetHeadWind(void) const
 {
   double psiw,vw;
 
-  psiw = Atmosphere->GetWindPsi();
-  vw = Atmosphere->GetTotalWindNED().Magnitude();
+  psiw = FDMExec->GetAtmosphere()->GetWindPsi();
+  vw = FDMExec->GetAtmosphere()->GetTotalWindNED().Magnitude();
 
-  return vw*cos(psiw - Propagate->GetEuler(ePsi));
+  return vw*cos(psiw - FDMExec->GetPropagate()->GetEuler(ePsi));
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -318,10 +310,17 @@ double FGAuxiliary::GetCrossWind(void) const
 {
   double psiw,vw;
 
-  psiw = Atmosphere->GetWindPsi();
-  vw = Atmosphere->GetTotalWindNED().Magnitude();
+  psiw = FDMExec->GetAtmosphere()->GetWindPsi();
+  vw = FDMExec->GetAtmosphere()->GetTotalWindNED().Magnitude();
+
+  return  vw*sin(psiw - FDMExec->GetPropagate()->GetEuler(ePsi));
+}
 
-  return  vw*sin(psiw - Propagate->GetEuler(ePsi));
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+double FGAuxiliary::GethVRP(void) const
+{
+  return vLocationVRP.GetRadius() - FDMExec->GetPropagate()->GetSeaLevelRadius();
 }
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -390,7 +389,7 @@ void FGAuxiliary::bind(void)
 
 void FGAuxiliary::CalculateRelativePosition(void)
 { 
-  const double earth_radius_mt = Inertial->GetRefRadius()*fttom;
+  const double earth_radius_mt = FDMExec->GetInertial()->GetRefRadius()*fttom;
   lat_relative_position=(FDMExec->GetPropagate()->GetLatitude()  - FDMExec->GetIC()->GetLatitudeDegIC() *degtorad)*earth_radius_mt;
   lon_relative_position=(FDMExec->GetPropagate()->GetLongitude() - FDMExec->GetIC()->GetLongitudeDegIC()*degtorad)*earth_radius_mt;
   relative_position = sqrt(lat_relative_position*lat_relative_position + lon_relative_position*lon_relative_position);