Use the density values from the environment subsystem, to properly handle

density variations due to humidity.

src/FDM/YASim/Airplane.cpp

@@ -667,7 +667,8 @@ void Airplane::runCruise()
 {
     setupState(_cruiseAoA, _cruiseSpeed, &_cruiseState);
     _model.setState(&_cruiseState);
-    _model.setAir(_cruiseP, _cruiseT);
+    _model.setAir(_cruiseP, _cruiseT,
+                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
 
     // The control configuration
     _controls.reset();
@@ -691,7 +692,8 @@ void Airplane::runCruise()
     for(i=0; i<_thrusters.size(); i++) {
 	Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
 	t->setWind(wind);
-	t->setAir(_cruiseP, _cruiseT);
+	t->setAir(_cruiseP, _cruiseT,
+                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
     }
     stabilizeThrust();
 
@@ -708,7 +710,8 @@ void Airplane::runApproach()
 {
     setupState(_approachAoA, _approachSpeed, &_approachState);
     _model.setState(&_approachState);
-    _model.setAir(_approachP, _approachT);
+    _model.setAir(_approachP, _approachT,
+                  Atmosphere::calcStdDensity(_approachP, _approachT));
 
     // The control configuration
     _controls.reset();
@@ -732,7 +735,8 @@ void Airplane::runApproach()
     for(i=0; i<_thrusters.size(); i++) {
 	Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
 	t->setWind(wind);
-	t->setAir(_approachP, _approachT);
+	t->setAir(_approachP, _approachT,
+                  Atmosphere::calcStdDensity(_approachP, _approachT));
     }
     stabilizeThrust();
 

src/FDM/YASim/Atmosphere.cpp

@@ -53,10 +53,11 @@ float Atmosphere::getStdDensity(float alt)
     return getRecord(alt, 3);
 }
 
-float Atmosphere::calcVEAS(float spd, float pressure, float temp)
+float Atmosphere::calcVEAS(float spd,
+                           float pressure, float temp, float density)
 {
     static float rho0 = getStdDensity(0);
-    float densityRatio = calcDensity(pressure, temp) / rho0;
+    float densityRatio = density / rho0;
     return spd * Math::sqrt(densityRatio);
 }
 
@@ -91,7 +92,7 @@ float Atmosphere::calcVCAS(float spd, float pressure, float temp)
     return Math::sqrt((7*p0/rho0)*(tmp-1));
 }
 
-float Atmosphere::calcDensity(float pressure, float temp)
+float Atmosphere::calcStdDensity(float pressure, float temp)
 {
     return pressure / (R * temp);
 }

src/FDM/YASim/Atmosphere.hpp

@@ -10,9 +10,9 @@ public:
     static float getStdDensity(float alt);
 
     static float calcVCAS(float spd, float pressure, float temp);
-    static float calcVEAS(float spd, float pressure, float temp);
+    static float calcVEAS(float spd, float pressure, float temp, float density);
     static float calcMach(float spd, float temp);
-    static float calcDensity(float pressure, float temp);
+    static float calcStdDensity(float pressure, float temp);
     
     // Given ambient ("0") pressure/density/temperature values,
     // calculate the properties of static air (air accelerated to the

src/FDM/YASim/Model.cpp

@@ -84,7 +84,7 @@ void Model::initIteration()
         localWind(pos, _s, v);
 
 	t->setWind(v);
-	t->setAir(_pressure, _temp);
+	t->setAir(_pressure, _temp, _rho);
 	t->integrate(_integrator.getInterval());
 
 	t->getTorque(v);
@@ -191,11 +191,11 @@ void Model::setGroundPlane(double* planeNormal, double fromOrigin)
     _ground[3] = fromOrigin;
 }
 
-void Model::setAir(float pressure, float temp)
+void Model::setAir(float pressure, float temp, float density)
 {
     _pressure = pressure;
     _temp = temp;
-    _rho = Atmosphere::calcDensity(pressure, temp);
+    _rho = density;
 }
 
 void Model::setWind(float* wind)

src/FDM/YASim/Model.hpp

@@ -50,7 +50,7 @@ public:
     void setGroundPlane(double* planeNormal, double fromOrigin);
     void setGroundEffect(float* pos, float span, float mul);
     void setWind(float* wind);
-    void setAir(float pressure, float temp);
+    void setAir(float pressure, float temp, float density);
 
     // BodyEnvironment callbacks
     virtual void calcForces(State* s);

src/FDM/YASim/Thruster.cpp

@@ -61,11 +61,11 @@ void Thruster::setWind(float* wind)
     for(i=0; i<3; i++) _wind[i] = wind[i];
 }
 
-void Thruster::setAir(float pressure, float temp)
+void Thruster::setAir(float pressure, float temp, float density)
 {
     _pressure = pressure;
     _temp = temp;
-    _rho = _pressure / (287.1f * _temp);
+    _rho = density;
 }
 
 }; // namespace yasim

src/FDM/YASim/Thruster.hpp

@@ -42,7 +42,7 @@ public:
 
     // Runtime instructions
     void setWind(float* wind);
-    void setAir(float pressure, float temp);
+    void setAir(float pressure, float temp, float density);
     virtual void init() {}
     virtual void integrate(float dt)=0;
     virtual void stabilize()=0;

src/FDM/YASim/YASim.cxx

@@ -28,6 +28,7 @@ static const float CM2GALS = 264.172037284; // gallons/cubic meter
 static const float KG2LBS = 2.20462262185;
 static const float W2HP = 1.3416e-3;
 static const float INHG2PA = 3386.389;
+static const float SLUG2KG = 14.59390;
 
 void YASim::printDEBUG()
 {
@@ -236,8 +237,10 @@ void YASim::copyToYASim(bool copyState)
     double ground = getACModel()->get3DModel()->getFGLocation()->get_cur_elev_m();
     // cout << "YASIM: ground = " << ground << endl;
 
-    float pressure = fgGetDouble("/environment/pressure-inhg") * INHG2PA;
+    float pressure = fgGetFloat("/environment/pressure-inhg") * INHG2PA;
-    float temp = fgGetDouble("/environment/temperature-degc") + 273.15;
+    float temp = fgGetFloat("/environment/temperature-degc") + 273.15;
+    float dens = fgGetFloat("/environment/density-slugft3") 
+        * SLUG2KG * M2FT*M2FT*M2FT;
 
     // Convert and set:
     Model* model = _fdm->getAirplane()->getModel();
@@ -274,7 +277,7 @@ void YASim::copyToYASim(bool copyState)
     model->setGroundPlane(gplane, rad);
 
     // air
-    model->setAir(pressure, temp);
+    model->setAir(pressure, temp, dens);
 }
 
 // All the settables:
@@ -366,10 +369,11 @@ void YASim::copyFromYASim()
     _set_Velocities_Wind_Body(v[0]*M2FT, -v[1]*M2FT, -v[2]*M2FT);
     _set_V_rel_wind(Math::mag3(v)*M2FT); // units?
 
-
     float P = fgGetDouble("/environment/pressure-inhg") * INHG2PA;
     float T = fgGetDouble("/environment/temperature-degC") + 273.15;
-    _set_V_equiv_kts(Atmosphere::calcVEAS(v[0], P, T)*MPS2KTS);
+    float D = fgGetFloat("/environment/density-slugft3")
+        *SLUG2KG * M2FT*M2FT*M2FT;
+    _set_V_equiv_kts(Atmosphere::calcVEAS(v[0], P, T, D)*MPS2KTS);
     _set_V_calibrated_kts(Atmosphere::calcVCAS(v[0], P, T)*MPS2KTS);
     _set_Mach_number(Atmosphere::calcMach(v[0], T));