diff --git a/src/FDM/YASim/Airplane.cpp b/src/FDM/YASim/Airplane.cpp
index 68ae0d9fb..7cbfe44c1 100644
--- a/src/FDM/YASim/Airplane.cpp
+++ b/src/FDM/YASim/Airplane.cpp
@@ -33,30 +33,7 @@ const float SOLVE_TWEAK = 0.3226;
 
 Airplane::Airplane()
 {
-    _emptyWeight = 0;
-    _pilotPos[0] = _pilotPos[1] = _pilotPos[2] = 0;
-    _wing = 0;
-    _tail = 0;
-    _ballast = 0;
-    _cruiseP = 0;
-    _cruiseT = 0;
-    _cruiseSpeed = 0;
-    _cruiseWeight = 0;
-    _cruiseGlideAngle = 0;
-    _approachP = 0;
-    _approachT = 0;
-    _approachSpeed = 0;
-    _approachAoA = 0;
-    _approachWeight = 0;
-    _approachGlideAngle = 0;
-
-    _dragFactor = 1;
-    _liftRatio = 1;
-    _cruiseAoA = 0;
-    _tailIncidence = 0;
-
-    _failureMsg = 0;
-    _wingsN = 0;
+    _approachConfig.isApproach = true;
 }
 
 Airplane::~Airplane()
@@ -82,10 +59,10 @@ Airplane::~Airplane()
     }
     for(i=0; i<_solveWeights.size(); i++)
         delete (SolveWeight*)_solveWeights.get(i);
-    for(i=0; i<_cruiseControls.size(); i++)
-        delete (Control*)_cruiseControls.get(i);
-    for(i=0; i<_approachControls.size(); i++) {
-        Control* c = (Control*)_approachControls.get(i);
+    for(i=0; i<_cruiseConfig.controls.size(); i++)
+        delete (Control*)_cruiseConfig.controls.get(i);
+    for(i=0; i<_approachConfig.controls.size(); i++) {
+        Control* c = (Control*)_approachConfig.controls.get(i);
         if(c != &_approachElevator)
             delete c;
     }
@@ -148,30 +125,29 @@ void Airplane::updateGearState()
 
 void Airplane::setApproach(float speed, float altitude, float aoa, float fuel, float gla)
 {
-    _approachSpeed = speed;
-    _approachP = Atmosphere::getStdPressure(altitude);
-    _approachT = Atmosphere::getStdTemperature(altitude);
-    _approachAoA = aoa;
-    _approachFuel = fuel;
-    _approachGlideAngle = gla;
+    _approachConfig.speed = speed;
+    _approachConfig.altitude = altitude;
+    _approachConfig.state.setupOrientationFromAoa(aoa); // see runConfig()
+    _approachConfig.aoa = aoa; // not strictly needed, see runConfig()
+    _approachConfig.fuel = fuel;
+    _approachConfig.glideAngle = gla;
 }
  
 void Airplane::setCruise(float speed, float altitude, float fuel, float gla)
 {
-    _cruiseSpeed = speed;
-    _cruiseP = Atmosphere::getStdPressure(altitude);
-    _cruiseT = Atmosphere::getStdTemperature(altitude);
-    _cruiseAoA = 0;
+    _cruiseConfig.speed = speed;
+    _cruiseConfig.altitude = altitude;
+    _cruiseConfig.aoa = 0;
     _tailIncidence = 0;
-    _cruiseFuel = fuel;
-    _cruiseGlideAngle = gla;
+    _cruiseConfig.fuel = fuel;
+    _cruiseConfig.glideAngle = gla;
 }
 
 void Airplane::setElevatorControl(int control)
 {
     _approachElevator.control = control;
     _approachElevator.val = 0;
-    _approachControls.add(&_approachElevator);
+    _approachConfig.controls.add(&_approachElevator);
 }
 
 void Airplane::addApproachControl(int control, float val)
@@ -179,7 +155,7 @@ void Airplane::addApproachControl(int control, float val)
     Control* c = new Control();
     c->control = control;
     c->val = val;
-    _approachControls.add(c);
+    _approachConfig.controls.add(c);
 }
 
 void Airplane::addCruiseControl(int control, float val)
@@ -187,7 +163,7 @@ void Airplane::addCruiseControl(int control, float val)
     Control* c = new Control();
     c->control = control;
     c->val = val;
-    _cruiseControls.add(c);
+    _cruiseConfig.controls.add(c);
 }
 
 void Airplane::addSolutionWeight(bool approach, int idx, float wgt)
@@ -300,31 +276,11 @@ void Airplane::setFuelFraction(float frac)
     }
 }
 
-void Airplane::setupState(const float aoa, const float speed, const float gla, State* s)
-{
-    float cosAoA = Math::cos(aoa);
-    float sinAoA = Math::sin(aoa);
-    s->orient[0] =  cosAoA; s->orient[1] = 0; s->orient[2] = sinAoA;
-    s->orient[3] =       0; s->orient[4] = 1; s->orient[5] =      0;
-    s->orient[6] = -sinAoA; s->orient[7] = 0; s->orient[8] = cosAoA;
-
-    //? what is gla? v[1]=y, v[2]=z? should sin go to v2 instead v1?
-    s->v[0] = speed*Math::cos(gla); s->v[1] = -speed*Math::sin(gla); s->v[2] = 0;
-
-    for(int i=0; i<3; i++)
-	s->pos[i] = s->rot[i] = s->acc[i] = s->racc[i] = 0;
-
-    // Put us 1m above the origin, or else the gravity computation in
-    // Model goes nuts
-    s->pos[2] = 1;
-}
-
 /**
  * @brief add contact point for crash detection 
  * used to add wingtips and fuselage nose and tail
  * 
  * @param pos ...
- * @return void
  */
 
 void Airplane::addContactPoint(float* pos)
@@ -638,9 +594,10 @@ void Airplane::compile()
         t->handle = body->addMass(0, t->pos);
         totalFuel += t->cap;
     }
-    _cruiseWeight = _emptyWeight + totalFuel*_cruiseFuel;
-    _approachWeight = _emptyWeight + totalFuel*_approachFuel;
+    _cruiseConfig.weight = _emptyWeight + totalFuel*_cruiseConfig.fuel;
+    _approachConfig.weight = _emptyWeight + totalFuel*_approachConfig.fuel;
 
+    
     body->recalc();
 
     // Add surfaces for the landing gear.
@@ -723,12 +680,12 @@ void Airplane::solveGear()
     // The force at max compression should be sufficient to stop a
     // plane moving downwards at 2x the approach descent rate.  Assume
     // a 3 degree approach.
-    float descentRate = 2.0f*_approachSpeed/19.1f;
+    float descentRate = 2.0f*_approachConfig.speed/19.1f;
 
     // Spread the kinetic energy according to the gear weights.  This
     // will results in an equal compression fraction (not distance) of
     // each gear.
-    float energy = 0.5f*_approachWeight*descentRate*descentRate;
+    float energy = 0.5f*_approachConfig.weight*descentRate*descentRate;
 
     for(i=0; i<_gears.size(); i++) {
         GearRec* gr = (GearRec*)_gears.get(i);
@@ -743,7 +700,7 @@ void Airplane::solveGear()
         gr->gear->setSpring(k * gr->gear->getSpring());
 
         // Critically damped (too damped, too!)
-        gr->gear->setDamping(2*Math::sqrt(k*_approachWeight*gr->wgt)
+        gr->gear->setDamping(2*Math::sqrt(k*_approachConfig.weight*gr->wgt)
                              * gr->gear->getDamping());
     }
 }
@@ -789,7 +746,7 @@ void Airplane::setupWeights(bool isApproach)
 }
 
 /// load values for controls as defined in cruise configuration
-void Airplane::loadControls(Vector& controls)
+void Airplane::loadControls(const Vector& controls)
 {
   _controls.reset();
   for(int i=0; i < controls.size(); i++) {
@@ -800,81 +757,43 @@ void Airplane::loadControls(Vector& controls)
 }
 
 /// Helper for solve()
-void Airplane::runCruise()
+void Airplane::runConfig(Config &cfg)
 {
-    setupState(_cruiseAoA, _cruiseSpeed,_cruiseGlideAngle, &_cruiseState);
-    _model.setState(&_cruiseState);
-    _model.setAir(_cruiseP, _cruiseT,
-                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
-
-    // The control configuration
-    loadControls(_cruiseControls);
-    
-    // The local wind
-    float wind[3];
-    Math::mul3(-1, _cruiseState.v, wind);
-    Math::vmul33(_cruiseState.orient, wind, wind);
- 
-    setFuelFraction(_cruiseFuel);
-    setupWeights(false);
-   
-    // Set up the thruster parameters and iterate until the thrust
-    // stabilizes.
-    for(int i=0; i<_thrusters.size(); i++) {
-	Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
-	t->setWind(wind);
-	t->setAir(_cruiseP, _cruiseT,
-                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
-    }
-    stabilizeThrust();
-
-    updateGearState();
-
-    // Precompute thrust in the model, and calculate aerodynamic forces
-    _model.getBody()->recalc();
-    _model.getBody()->reset();
-    _model.initIteration();
-    _model.calcForces(&_cruiseState);
+  // aoa is consider to be given for approach so we calculate orientation 
+  // only once in setApproach()
+  if (!cfg.isApproach) {
+    cfg.state.setupOrientationFromAoa(cfg.aoa);
+  }
+  cfg.state.setupSpeedAndPosition(cfg.speed, cfg.glideAngle);
+  _model.setState(&cfg.state);
+  _model.setStandardAtmosphere(cfg.altitude);
+  loadControls(cfg.controls);
+  
+  // The local wind
+  float wind[3];
+  Math::mul3(-1, cfg.state.v, wind);
+  cfg.state.globalToLocal(wind, wind);
+  
+  setFuelFraction(cfg.fuel);
+  setupWeights(cfg.isApproach);
+  
+  // Set up the thruster parameters and iterate until the thrust
+  // stabilizes.
+  for(int i=0; i<_thrusters.size(); i++) {
+    Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
+    t->setWind(wind);
+    t->setStandardAtmosphere(cfg.altitude);
+  }
+  
+  stabilizeThrust();
+  updateGearState();
+  
+  // Precompute thrust in the model, and calculate aerodynamic forces
+  _model.getBody()->recalc();
+  _model.getBody()->reset();
+  _model.initIteration();
+  _model.calcForces(&cfg.state);
 }
-
-/// Helper for solve()
-void Airplane::runApproach()
-{
-    setupState(_approachAoA, _approachSpeed,_approachGlideAngle, &_approachState);
-    _model.setState(&_approachState);
-    _model.setAir(_approachP, _approachT,
-                  Atmosphere::calcStdDensity(_approachP, _approachT));
-
-    // The control configuration
-    loadControls(_approachControls);
-    
-    // The local wind
-    float wind[3];
-    Math::mul3(-1, _approachState.v, wind);
-    Math::vmul33(_approachState.orient, wind, wind);
-    
-    setFuelFraction(_approachFuel);
-    setupWeights(true);
-
-    // Run the thrusters until they get to a stable setting.  FIXME:
-    // this is lots of wasted work.
-    for(int i=0; i<_thrusters.size(); i++) {
-	Thruster* t = ((ThrustRec*)_thrusters.get(i))->thruster;
-	t->setWind(wind);
-	t->setAir(_approachP, _approachT,
-                  Atmosphere::calcStdDensity(_approachP, _approachT));
-    }
-    stabilizeThrust();
-
-    updateGearState();
-
-    // Precompute thrust in the model, and calculate aerodynamic forces
-    _model.getBody()->recalc();
-    _model.getBody()->reset();
-    _model.initIteration();
-    _model.calcForces(&_approachState);
-}
-
 /// Used only in Airplane::solve() and solveHelicopter(), not at runtime
 void Airplane::applyDragFactor(float factor)
 {
@@ -960,51 +879,51 @@ void Airplane::solve()
         }
 
 	// Run an iteration at cruise, and extract the needed numbers:
-	runCruise();
+	runConfig(_cruiseConfig);
 
 	_model.getThrust(tmp);
-        float thrust = tmp[0] + _cruiseWeight * Math::sin(_cruiseGlideAngle) * 9.81;
+        float thrust = tmp[0] + _cruiseConfig.weight * Math::sin(_cruiseConfig.glideAngle) * 9.81;
 
 	_model.getBody()->getAccel(tmp);
-        Math::tmul33(_cruiseState.orient, tmp, tmp);
-	float xforce = _cruiseWeight * tmp[0];
-	float clift0 = _cruiseWeight * tmp[2];
+        _cruiseConfig.state.localToGlobal(tmp, tmp);
+	float xforce = _cruiseConfig.weight * tmp[0];
+	float clift0 = _cruiseConfig.weight * tmp[2];
 
 	_model.getBody()->getAngularAccel(tmp);
-        Math::tmul33(_cruiseState.orient, tmp, tmp);
+        _cruiseConfig.state.localToGlobal(tmp, tmp);
 	float pitch0 = tmp[1];
 
 	// Run an approach iteration, and do likewise
-	runApproach();
+        runConfig(_approachConfig);
 
 	_model.getBody()->getAngularAccel(tmp);
-        Math::tmul33(_approachState.orient, tmp, tmp);
+        _approachConfig.state.localToGlobal(tmp, tmp);
 	double apitch0 = tmp[1];
 
 	_model.getBody()->getAccel(tmp);
-        Math::tmul33(_approachState.orient, tmp, tmp);
-	float alift = _approachWeight * tmp[2];
+        _approachConfig.state.localToGlobal(tmp, tmp);
+	float alift = _approachConfig.weight * tmp[2];
 
 	// Modify the cruise AoA a bit to get a derivative
-	_cruiseAoA += ARCMIN;
-	runCruise();
-	_cruiseAoA -= ARCMIN;
-
+	_cruiseConfig.aoa += ARCMIN;
+        runConfig(_cruiseConfig);
+        _cruiseConfig.aoa -= ARCMIN;
+        
 	_model.getBody()->getAccel(tmp);
-        Math::tmul33(_cruiseState.orient, tmp, tmp);
-	float clift1 = _cruiseWeight * tmp[2];
+        _cruiseConfig.state.localToGlobal(tmp, tmp);
+	float clift1 = _cruiseConfig.weight * tmp[2];
 
 	// Do the same with the tail incidence
 	_tail->setIncidence(_tailIncidence + ARCMIN);
-	runCruise();
-	_tail->setIncidence(_tailIncidence);
+        runConfig(_cruiseConfig);
+        _tail->setIncidence(_tailIncidence);
 
 	_model.getBody()->getAngularAccel(tmp);
-        Math::tmul33(_cruiseState.orient, tmp, tmp);
+        _cruiseConfig.state.localToGlobal(tmp, tmp);
 	float pitch1 = tmp[1];
 
 	// Now calculate:
-	float awgt = 9.8f * _approachWeight;
+	float awgt = 9.8f * _approachConfig.weight;
 
 	float dragFactor = thrust / (thrust-xforce);
 	float liftFactor = awgt / (awgt+alift);
@@ -1021,11 +940,11 @@ void Airplane::solve()
         // variable).
         const float ELEVDIDDLE = 0.001f;
         _approachElevator.val += ELEVDIDDLE;
-        runApproach();
+        runConfig(_approachConfig);
         _approachElevator.val -= ELEVDIDDLE;
 
 	_model.getBody()->getAngularAccel(tmp);
-        Math::tmul33(_approachState.orient, tmp, tmp);
+        _approachConfig.state.localToGlobal(tmp, tmp);
 	double apitch1 = tmp[1];
         float elevDelta = -apitch0 * (ELEVDIDDLE/(apitch1-apitch0));
 
@@ -1044,14 +963,14 @@ void Airplane::solve()
 	}
 
 	// OK, now we can adjust the minor variables:
-	_cruiseAoA += SOLVE_TWEAK*aoaDelta;
+	_cruiseConfig.aoa += SOLVE_TWEAK*aoaDelta;
 	_tailIncidence += SOLVE_TWEAK*tailDelta;
 	
-	_cruiseAoA = Math::clamp(_cruiseAoA, -0.175f, 0.175f);
+	_cruiseConfig.aoa = Math::clamp(_cruiseConfig.aoa, -0.175f, 0.175f);
 	_tailIncidence = Math::clamp(_tailIncidence, -0.175f, 0.175f);
 
-        if(abs(xforce/_cruiseWeight) < STHRESH*0.0001 &&
-           abs(alift/_approachWeight) < STHRESH*0.0001 &&
+        if(abs(xforce/_cruiseConfig.weight) < STHRESH*0.0001 &&
+           abs(alift/_approachConfig.weight) < STHRESH*0.0001 &&
            abs(aoaDelta) < STHRESH*.000017 &&
            abs(tailDelta) < STHRESH*.000017)
         {
@@ -1074,7 +993,7 @@ void Airplane::solve()
     } else if(_liftRatio < 1e-04 || _liftRatio > 1e4) {
 	_failureMsg = "Lift ratio beyond reasonable bounds.";
 	return;
-    } else if(Math::abs(_cruiseAoA) >= .17453293) {
+    } else if(Math::abs(_cruiseConfig.aoa) >= .17453293) {
 	_failureMsg = "Cruise AoA > 10 degrees";
 	return;
     } else if(Math::abs(_tailIncidence) >= .17453293) {
@@ -1102,14 +1021,12 @@ void Airplane::solveHelicopter()
         applyDragFactor(Math::pow(15.7/1000, 1/SOLVE_TWEAK));
         applyLiftRatio(Math::pow(104, 1/SOLVE_TWEAK));
     }
-    setupState(0,0,0, &_cruiseState);
-    _model.setState(&_cruiseState);
+    _cruiseConfig.state.setupState(0,0,0);
+    _model.setState(&_cruiseConfig.state);
     setupWeights(true);
     _controls.reset();
     _model.getBody()->reset();
-    _model.setAir(_cruiseP, _cruiseT,
-                  Atmosphere::calcStdDensity(_cruiseP, _cruiseT));
-    
+    _model.setStandardAtmosphere(_cruiseConfig.altitude);    
 }
 
 float Airplane::getCGMAC()
diff --git a/src/FDM/YASim/Airplane.hpp b/src/FDM/YASim/Airplane.hpp
index 8ee625fd5..3aee7c224 100644
--- a/src/FDM/YASim/Airplane.hpp
+++ b/src/FDM/YASim/Airplane.hpp
@@ -95,14 +95,13 @@ public:
     int getSolutionIterations() const { return _solutionIterations; }
     float getDragCoefficient() const { return _dragFactor; }
     float getLiftRatio() const { return _liftRatio; }
-    float getCruiseAoA() const { return _cruiseAoA; }
+    float getCruiseAoA() const { return _cruiseConfig.aoa; }
     float getTailIncidence() const { return _tailIncidence; }
     float getApproachElevator() const { return _approachElevator.val; }
     const char* getFailureMsg() const { return _failureMsg; }
 
-    static void setupState(const float aoa, const float speed, const float gla, yasim::State* s); // utility
-    void loadApproachControls() { loadControls(_approachControls); }
-    void loadCruiseControls() { loadControls(_cruiseControls); }
+    void loadApproachControls() { loadControls(_approachConfig.controls); }
+    void loadCruiseControls() { loadControls(_cruiseConfig.controls); }
     
     float getCGHardLimitXMin() const { return _cgMin; } // get min x-coordinate for c.g (from main gear)
     float getCGHardLimitXMax() const { return _cgMax; } // get max x-coordinate for c.g (from nose gear)
@@ -114,21 +113,60 @@ public:
     void  setAutoBallast(bool allowed) { _autoBallast = allowed; }
     
 private:
-    struct Tank { float pos[3]; float cap; float fill;
-        float density; int handle; };
-    struct Fuselage { float front[3], back[3], width, taper, mid, _cx, _cy, _cz, _idrag;
-        Vector surfs; };
-    struct GearRec { Gear* gear; Surface* surf; float wgt; };
-    struct ThrustRec { Thruster* thruster;
-        int handle; float cg[3]; float mass; };
-    struct Control { int control; float val; };
-    struct WeightRec { int handle; Surface* surf; };
-    struct SolveWeight { bool approach; int idx; float wgt; };
-    struct ContactRec { Gear* gear; float p[3]; };
+    struct Tank { 
+      float pos[3]; 
+      float cap, fill, density;
+      int handle;       
+    };
+    struct Fuselage { 
+      float  front[3], back[3];
+      float  width, taper, mid, _cx, _cy, _cz, _idrag;
+      Vector surfs;      
+    };
+    struct GearRec { 
+      Gear* gear;
+      Surface* surf;
+      float wgt;
+    };
+    struct ThrustRec { 
+      Thruster* thruster;
+      int handle;
+      float cg[3];
+      float mass;
+    };
+    struct Control { 
+      int control;
+      float val;
+    };
+    struct WeightRec { 
+      int handle;
+      Surface* surf;
+    };
+    struct SolveWeight { 
+      bool approach;
+      int idx;
+      float wgt;
+    };
+    struct ContactRec {
+      Gear* gear;
+      float p[3];
+    };
+    struct Config {
+      bool isApproach {false};
+      float speed {0};
+      float fuel {0};
+      float glideAngle {0};      
+      float aoa {0};
+      float altitude;
+      float weight;
+      State state;
+      Vector controls;
+    };
+    Config _cruiseConfig;
+    Config _approachConfig;
 
-    void loadControls(Vector &controls);
-    void runCruise();
-    void runApproach();
+    void loadControls(const Vector& controls);
+    void runConfig(Config &cfg);
     void solveGear();
     void solve();
     void solveHelicopter();
@@ -148,17 +186,17 @@ private:
     Model _model;
     ControlMap _controls;
 
-    float _emptyWeight;
-    float _pilotPos[3];
+    float _emptyWeight {0};
+    float _pilotPos[3] {0, 0, 0};
 
-    Wing* _wing;
-    Wing* _tail;
+    Wing* _wing {nullptr};
+    Wing* _tail {nullptr};
 
     Vector _fuselages;
     Vector _vstabs;
     Vector _tanks;
     Vector _thrusters;
-    float _ballast;
+    float _ballast {0};
 
     Vector _gears;
     Vector _contacts; // non-gear ground contact points
@@ -168,39 +206,19 @@ private:
 
     Vector _solveWeights;
 
-    Vector _cruiseControls;
-    State _cruiseState;
-    float _cruiseP;
-    float _cruiseT;
-    float _cruiseSpeed;
-    float _cruiseWeight;
-    float _cruiseFuel;
-    float _cruiseGlideAngle;
-
-    Vector _approachControls;
-    State _approachState;
-    float _approachP;
-    float _approachT;
-    float _approachSpeed;
-    float _approachAoA;
-    float _approachWeight;
-    float _approachFuel;
-    float _approachGlideAngle;
-
-    int _solutionIterations;
-    float _dragFactor;
-    float _liftRatio;
-    float _cruiseAoA;
-    float _tailIncidence;
+    int _solutionIterations {0};
+    float _dragFactor {1};
+    float _liftRatio {1};
+    float _tailIncidence {0};
     Control _approachElevator;
-    const char* _failureMsg;
+    const char* _failureMsg {0};
     
-    float _cgMax = -1e6;         // hard limits for cg from gear position
-    float _cgMin = 1e6;          // hard limits for cg from gear position
-    float _cgDesiredMax = 0.3f;  // desired cg max in %MAC from config
-    float _cgDesiredMin = 0.25f; // desired cg min in %MAC from config
-    float _cgDesiredFront; // calculated desired cg x max
-    float _cgDesiredAft;   // calculated desired cg x min 
+    float _cgMax {-1e6};         // hard limits for cg from gear position
+    float _cgMin {1e6};          // hard limits for cg from gear position
+    float _cgDesiredMax {0.3f};  // desired cg max in %MAC from config
+    float _cgDesiredMin {0.25f}; // desired cg min in %MAC from config
+    float _cgDesiredFront {0};   // calculated desired cg x max
+    float _cgDesiredAft {0};     // calculated desired cg x min 
     bool _autoBallast = false;
 };
 
diff --git a/src/FDM/YASim/Atmosphere.cpp b/src/FDM/YASim/Atmosphere.cpp
index 5fb7347f9..5e751dbce 100644
--- a/src/FDM/YASim/Atmosphere.cpp
+++ b/src/FDM/YASim/Atmosphere.cpp
@@ -1,7 +1,8 @@
+#include <string>
 #include "Math.hpp"
 #include "Atmosphere.hpp"
-namespace yasim {
 
+namespace yasim {
 // Copied from McCormick, who got it from "The ARDC Model Atmosphere"
 // Note that there's an error in the text in the first entry,
 // McCormick lists 299.16/101325/1.22500, but those don't agree with
@@ -9,7 +10,7 @@ namespace yasim {
 // pretty hot for a "standard" atmosphere.
 // Numbers above 19000 meters calculated from src/Environment/environment.cxx
 //                             meters   kelvin      Pa   kg/m^3
-float Atmosphere::data[][4] = {{ -900.0f, 293.91f, 111679.0f, 1.32353f },
+float Atmosphere::data[][Atmosphere::numColumns] = {{ -900.0f, 293.91f, 111679.0f, 1.32353f },
                                {    0.0f, 288.11f, 101325.0f, 1.22500f },
 			       {   900.0f, 282.31f,  90971.0f, 1.12260f },
 			       {  1800.0f, 276.46f,  81494.0f, 1.02690f },
@@ -48,28 +49,34 @@ float Atmosphere::data[][4] = {{ -900.0f, 293.91f, 111679.0f, 1.32353f },
 
 // Universal gas constant for air, in SI units.  P = R * rho * T.
 // P in pascals (N/m^2), rho is kg/m^3, T in kelvin.
-const float R = 287.1f;
+const float R = 287.058f;
 
 // Specific heat ratio for air, at "low" temperatures.  
 const float GAMMA = 1.4f;
 
+void Atmosphere::setStandard(float altitude)
+{
+    _density = getStdDensity(altitude);
+    _pressure = getStdPressure(altitude);
+    _temperature = getStdTemperature(altitude);
+}
+
 float Atmosphere::getStdTemperature(float alt)
 {
-    return getRecord(alt, 1);
+    return getRecord(alt, TEMPERATURE);
 }
 
 float Atmosphere::getStdPressure(float alt)
 {
-    return getRecord(alt, 2);
+    return getRecord(alt, PRESSURE);
 }
 
 float Atmosphere::getStdDensity(float alt)
 {
-    return getRecord(alt, 3);
+    return getRecord(alt, DENSITY);
 }
 
-float Atmosphere::calcVEAS(float spd,
-                           float pressure, float temp, float density)
+float Atmosphere::calcVEAS(float spd, float pressure, float temp, float density)
 {
     static float rho0 = getStdDensity(0);
     float densityRatio = density / rho0;
@@ -122,9 +129,14 @@ float Atmosphere::spdFromMach(float mach, float temp)
     return mach * Math::sqrt(GAMMA * R * temp);
 }
 
+float Atmosphere::spdFromMach(float mach)
+{
+    return spdFromMach(mach, _temperature);
+}
+
 float Atmosphere::spdFromVCAS(float vcas, float pressure, float temp)
 {
-                                // FIXME: does not account for supersonic
+    // FIXME: does not account for supersonic
     float p0 = getStdPressure(0);
     float rho0 = getStdDensity(0);
 
@@ -136,6 +148,11 @@ float Atmosphere::spdFromVCAS(float vcas, float pressure, float temp)
     return vtas;
 }
 
+float Atmosphere::spdFromVCAS(float vcas)
+{
+    return spdFromVCAS(vcas, _pressure, _temperature);
+}
+
 void Atmosphere::calcStaticAir(float p0, float t0, float d0, float v,
                                float* pOut, float* tOut, float* dOut)
 {
@@ -147,14 +164,24 @@ void Atmosphere::calcStaticAir(float p0, float t0, float d0, float v,
     *pOut = (*dOut) * R * (*tOut);
 }
 
-float Atmosphere::getRecord(float alt, int recNum)
+void Atmosphere::calcStaticAir(float v, float* pOut, float* tOut, float* dOut)
 {
-    int hi = (sizeof(data) / (4*sizeof(float))) - 1;
+    return calcStaticAir(_pressure, _temperature, _density, v, pOut, tOut, dOut);
+}
+
+
+float Atmosphere::getRecord(float alt, Column recNum)
+{
+    int hi = maxTableIndex();
     int lo = 0;
 
     // safety valve, clamp to the edges of the table
-    if(alt < data[0][0])       hi=1;
-    else if(alt > data[hi][0]) lo = hi-1;
+    if(alt < data[0][ALTITUDE]) {
+      hi = 1;
+    }
+    else if(alt > data[hi][ALTITUDE]) {
+      lo = hi-1;
+    }
 
     // binary search
     while(1) {
@@ -165,10 +192,38 @@ float Atmosphere::getRecord(float alt, int recNum)
     }
 
     // interpolate
-    float frac = (alt - data[lo][0])/(data[hi][0] - data[lo][0]);
+    float frac = (alt - data[lo][ALTITUDE])/(data[hi][ALTITUDE] - data[lo][ALTITUDE]);
     float a = data[lo][recNum];
     float b = data[hi][recNum];
     return a + frac * (b-a);
 }
 
+int Atmosphere::maxTableIndex() {
+    return (sizeof(data) / (numColumns * sizeof(float))) - 1;
+}
+
+bool Atmosphere::test() {
+    bool passed = true;
+    int rows = maxTableIndex() + 1;
+    const float maxDeviation = 0.0002f;
+    
+    fprintf(stderr, "Atmosphere::test()\n");
+    fprintf(stderr, "Columns = %d\n", numColumns);
+    fprintf(stderr, "Rows = %d\n", rows);
+    
+    for (int alt = 0; alt < maxTableIndex(); alt++) {
+      float density = calcStdDensity(data[alt][PRESSURE], data[alt][TEMPERATURE]);
+      float delta = data[alt][DENSITY] - density;
+      fprintf(stderr, "%d : %f \n", alt, delta);
+      if (Math::abs(delta) > maxDeviation) { 
+        passed = false;       
+        fprintf(stderr,"FAIL: Deviation above limit of %1.6f\n", maxDeviation);
+      }
+    }
+    if (passed) {
+      fprintf(stderr,"Deviation below %1.6f for all rows.\n", maxDeviation);
+    }
+    return passed;
+}
+
 }; // namespace yasim
diff --git a/src/FDM/YASim/Atmosphere.hpp b/src/FDM/YASim/Atmosphere.hpp
index b0cfd3d86..bcac90728 100644
--- a/src/FDM/YASim/Atmosphere.hpp
+++ b/src/FDM/YASim/Atmosphere.hpp
@@ -3,8 +3,28 @@
 
 namespace yasim {
 
+//constexpr int Atmosphere::numColumns {4};
+  
 class Atmosphere {
+    enum Column {
+      ALTITUDE,
+      TEMPERATURE,
+      PRESSURE,
+      DENSITY
+    };
+    static const int numColumns {4};
 public:
+    void setTemperature(float t) { _temperature = t; }
+    void setPressure(float p) { _pressure = p; }
+    void setDensity(float d) { _density = d; }
+
+    //set temperature, pressure and density to standard values for given altitude
+    void setStandard(float altitude);
+    
+    float getTemperature() const { return _temperature; }
+    float getPressure() const { return _pressure; }
+    float getDensity() const { return _density; }
+    
     static float getStdTemperature(float alt);
     static float getStdPressure(float alt);
     static float getStdDensity(float alt);
@@ -16,6 +36,8 @@ public:
 
     static float spdFromMach(float mach, float temp);
     static float spdFromVCAS(float vcas, float pressure, float temp);
+    float spdFromMach(float mach);
+    float spdFromVCAS(float vcas);
     
     // Given ambient ("0") pressure/density/temperature values,
     // calculate the properties of static air (air accelerated to the
@@ -23,10 +45,17 @@ public:
     // compressibility, but not shock effects.
     static void calcStaticAir(float p0, float t0, float d0, float v,
                               float* pOut, float* tOut, float* dOut);
-
+    void calcStaticAir(float v, float* pOut, float* tOut, float* dOut);
+    static bool test();
+    
 private:
-    static float getRecord(float alt, int idx);
-    static float data[][4];
+    static float getRecord(float alt, Atmosphere::Column recNum);
+    static float data[][numColumns];
+    static int maxTableIndex();
+    
+    float _temperature = 288.11f;
+    float _pressure = 101325.0f;
+    float _density = 1.22500f;
 };
 
 }; // namespace yasim
diff --git a/src/FDM/YASim/BodyEnvironment.hpp b/src/FDM/YASim/BodyEnvironment.hpp
index 0ff096a42..892894d23 100644
--- a/src/FDM/YASim/BodyEnvironment.hpp
+++ b/src/FDM/YASim/BodyEnvironment.hpp
@@ -13,45 +13,41 @@ namespace yasim {
 // You can get local->global transformations by calling Math::tmul33()
 // and using the same matrix.
 struct State {
-    double pos[3];    // position
-    float  orient[9]; // global->local xform matrix
-    float  v[3];      // velocity
-    float  rot[3];    // rotational velocity
-    float  acc[3];    // acceleration
-    float  racc[3];   // rotational acceleration
+    double pos[3] {0, 0, 0};    // position
+    float  orient[9]  {1,0,0,  0,1,0,  0,0,1}; // global->local xform matrix
+    float  v[3] {0, 0, 0};      // velocity
+    float  rot[3] {0, 0, 0};    // rotational velocity
+    float  acc[3] {0, 0, 0};    // acceleration
+    float  racc[3] {0, 0, 0};   // rotational acceleration
 
     // Simple initialization
     State() {
-        int i;
-        for(i=0; i<3; i++) {
-            pos[i] = v[i] = rot[i] = acc[i] = racc[i] = 0;
-            int j;
-            for(j=0; j<3; j++)
-                orient[3*i+j] = i==j ? 1.0f : 0.0f;
-        }
     }
 
-    void posLocalToGlobal(float* lpos, double *gpos) {
+    void posLocalToGlobal(const float* lpos, double *gpos) const {
         float tmp[3];
         Math::tmul33(orient, lpos, tmp);
         gpos[0] = tmp[0] + pos[0];
         gpos[1] = tmp[1] + pos[1];
         gpos[2] = tmp[2] + pos[2];
     }
-    void posGlobalToLocal(double* gpos, float *lpos) {
+    
+    void posGlobalToLocal(const double* gpos, float *lpos) const {
         lpos[0] = (float)(gpos[0] - pos[0]);
         lpos[1] = (float)(gpos[1] - pos[1]);
         lpos[2] = (float)(gpos[2] - pos[2]);
         Math::vmul33(orient, lpos, lpos);
     }
-    void velLocalToGlobal(float* lvel, float *gvel) {
-        Math::tmul33(orient, lvel, gvel);
+    
+    void localToGlobal(const float* local, float *global) const {
+        Math::tmul33(orient, local, global);
     }
-    void velGlobalToLocal(float* gvel, float *lvel) {
-        Math::vmul33(orient, gvel, lvel);
+    
+    void globalToLocal(const float* global, float *local) const {
+        Math::vmul33(orient, global, local);
     }
 
-    void planeGlobalToLocal(double* gplane, float *lplane) {
+    void planeGlobalToLocal(const double* gplane, float *lplane) const {
       // First the normal vector transformed to local coordinates.
       lplane[0] = (float)-gplane[0];
       lplane[1] = (float)-gplane[1];
@@ -62,7 +58,35 @@ struct State {
       lplane[3] = (float)(pos[0]*gplane[0] + pos[1]*gplane[1]
                           + pos[2]*gplane[2] - gplane[3]);
     }
+    
+    // used by Airplane::runCruise, runApproach, solveHelicopter and in yasim-test
+    void setupOrientationFromAoa(float aoa) 
+    {
+      float cosAoA = Math::cos(aoa);
+      float sinAoA = Math::sin(aoa);
+      orient[0] =  cosAoA; orient[1] = 0; orient[2] = sinAoA;
+      orient[3] =       0; orient[4] = 1; orient[5] =      0;
+      orient[6] = -sinAoA; orient[7] = 0; orient[8] = cosAoA;
+    }
+    
+    void setupSpeedAndPosition(float speed, float gla) 
+    {
+      
+      // FIXME check axis, guess sin should go to 2 instead of 1?
+      v[0] = speed*Math::cos(gla); 
+      v[1] = -speed*Math::sin(gla); 
+      v[2] = 0;
+      for(int i=0; i<3; i++) {
+        pos[i] = rot[i] = acc[i] = racc[i] = 0;
+      }
+      
+      pos[2] = 1;
+    }
 
+    void setupState(float aoa, float speed, float gla) {
+        setupOrientationFromAoa(aoa);
+        setupSpeedAndPosition(speed, gla);
+    }
 };
 
 //
diff --git a/src/FDM/YASim/CMakeLists.txt b/src/FDM/YASim/CMakeLists.txt
index 8492a3d44..0d7efa881 100644
--- a/src/FDM/YASim/CMakeLists.txt
+++ b/src/FDM/YASim/CMakeLists.txt
@@ -22,7 +22,6 @@ set(COMMON
 	Rotorpart.cpp
 	SimpleJet.cpp
 	Surface.cpp
-	Thruster.cpp
 	TurbineEngine.cpp
 	Turbulence.cpp
 	Wing.cpp
@@ -40,6 +39,7 @@ flightgear_component(YASim  "${SOURCES}")
 if(ENABLE_TESTS)
 add_executable(yasim yasim-test.cpp ${COMMON})
 add_executable(yasim-proptest proptest.cpp ${COMMON})
+add_executable(yasim-atmotest yasim-atmotest.cpp Atmosphere.cpp )
 
 target_link_libraries(yasim SimGearCore)
 target_link_libraries(yasim-proptest SimGearCore)
diff --git a/src/FDM/YASim/FGFDM.cpp b/src/FDM/YASim/FGFDM.cpp
index 8bcd6437f..e5e3f51f5 100644
--- a/src/FDM/YASim/FGFDM.cpp
+++ b/src/FDM/YASim/FGFDM.cpp
@@ -27,25 +27,6 @@
 
 namespace yasim {
 
-// Some conversion factors
-static const float KTS2MPS = 0.514444444444;
-static const float KMH2MPS = 1/3.6;
-static const float FT2M = 0.3048;
-static const float DEG2RAD = 0.0174532925199;
-static const float RPM2RAD = 0.10471975512;
-static const float LBS2N = 4.44822;
-static const float LBS2KG = 0.45359237;
-static const float KG2LBS = 2.2046225;
-static const float CM2GALS = 264.172037284;
-static const float HP2W = 745.700;
-static const float INHG2PA = 3386.389;
-static const float K2DEGF = 1.8;
-static const float K2DEGFOFFSET = -459.4;
-static const float CIN2CM = 1.6387064e-5;
-static const float YASIM_PI = 3.14159265358979323846;
-
-static const float NM2FTLB = (1/(LBS2N*FT2M));
-
 // Stubs, so that this can be compiled without the FlightGear
 // binary.  What's the best way to handle this?
 
diff --git a/src/FDM/YASim/FGFDM.hpp b/src/FDM/YASim/FGFDM.hpp
index c888202d2..28eabba21 100644
--- a/src/FDM/YASim/FGFDM.hpp
+++ b/src/FDM/YASim/FGFDM.hpp
@@ -4,6 +4,7 @@
 #include <simgear/xml/easyxml.hxx>
 #include <simgear/props/props.hxx>
 
+#include "yasim-common.hpp"
 #include "Airplane.hpp"
 #include "Vector.hpp"
 
diff --git a/src/FDM/YASim/Gear.cpp b/src/FDM/YASim/Gear.cpp
index 74c2a9a04..7cb066660 100644
--- a/src/FDM/YASim/Gear.cpp
+++ b/src/FDM/YASim/Gear.cpp
@@ -165,7 +165,7 @@ void Gear::calcForce(RigidBody* body, State *s, float* v, float* rot)
 
     // The velocity of the contact patch transformed to local coordinates.
     float glvel[3];
-    s->velGlobalToLocal(_global_vel, glvel);
+    s->globalToLocal(_global_vel, glvel);
 
     // First off, make sure that the gear "tip" is below the ground.
     // If it's not, there's no force.
diff --git a/src/FDM/YASim/Hitch.cpp b/src/FDM/YASim/Hitch.cpp
index b3d3bfb1c..cb4abd56e 100644
--- a/src/FDM/YASim/Hitch.cpp
+++ b/src/FDM/YASim/Hitch.cpp
@@ -413,7 +413,7 @@ void Hitch::calcForce(Ground *g_cb, RigidBody* body, State* s)
                                   //With this trick, both player in aerotow get the same length
     Math::unit3(delta,deltaN);
     float lvel[3];
-    s->velGlobalToLocal(s->v,lvel);
+    s->globalToLocal(s->v,lvel);
     _speed_in_tow_direction=Math::dot3(lvel,deltaN);
     if (_towEndIsConnectedToProperty && _nodeIsMultiplayer)
     {
@@ -421,7 +421,7 @@ void Hitch::calcForce(Ground *g_cb, RigidBody* body, State* s)
         _timeLagCorrectedDist=_dist+mp_delta_dist_due_to_time_lag;
         if(_forceIsCalculatedByMaster && !_open)
         {
-            s->velGlobalToLocal(_mp_force,_force);
+            s->globalToLocal(_mp_force,_force);
             return;
         }
     }
@@ -528,7 +528,7 @@ void Hitch::calcForce(Ground *g_cb, RigidBody* body, State* s)
     //the same for the tow end:
     Math::mul3(grav_frac_tow_end*grav_force,ground,grav_force_v);
     Math::add3(grav_force_v,_towEndForce,_towEndForce);
-    s->velLocalToGlobal(_towEndForce,_towEndForce);
+    s->localToGlobal(_towEndForce,_towEndForce);
 
     if(_forceMagnitude>=_towBrakeForce)
     {
diff --git a/src/FDM/YASim/Hook.cpp b/src/FDM/YASim/Hook.cpp
index 092404883..0536e3432 100644
--- a/src/FDM/YASim/Hook.cpp
+++ b/src/FDM/YASim/Hook.cpp
@@ -216,8 +216,8 @@ void Hook::calcForce(Ground* g_cb, RigidBody* body, State* s, float* lv, float*
     float wire_lpos[2][3];
     s->posGlobalToLocal(dpos[0], wire_lpos[0]);
     s->posGlobalToLocal(dpos[1], wire_lpos[1]);
-    s->velGlobalToLocal(wire_vel[0], wire_vel[0]);
-    s->velGlobalToLocal(wire_vel[1], wire_vel[1]);
+    s->globalToLocal(wire_vel[0], wire_vel[0]);
+    s->globalToLocal(wire_vel[1], wire_vel[1]);
 
     // Compute the velocity of the hooks mount point in the local frame.
     float mount_vel[3];
diff --git a/src/FDM/YASim/Integrator.cpp b/src/FDM/YASim/Integrator.cpp
index efb080c04..57bf1c93c 100644
--- a/src/FDM/YASim/Integrator.cpp
+++ b/src/FDM/YASim/Integrator.cpp
@@ -2,43 +2,6 @@
 #include "Integrator.hpp"
 namespace yasim {
 
-void Integrator::setBody(RigidBody* body)
-{
-    _body = body;
-}
-
-void Integrator::setEnvironment(BodyEnvironment* env)
-{
-    _env = env;
-}
-
-void Integrator::setInterval(float dt)
-{
-    _dt = dt;
-}
-
-float Integrator::getInterval()
-{
-    return _dt;
-}
-
-void Integrator::setState(State* s)
-{
-    _s = *s;
-}
-
-State* Integrator::getState()
-{
-    return &_s;
-}
-
-// Transforms a "local" vector to a "global" vector (not coordinate!)
-// using the specified orientation.
-void Integrator::l2gVector(float* orient, float* v, float* out)
-{
-    Math::tmul33(orient, v, out);
-}
-
 // Updates a position vector for a body c.g. motion with velocity v
 // over time dt, from orientation o0 to o1.  Because the position
 // references the local coordinate origin, but the velocity is that of
diff --git a/src/FDM/YASim/Integrator.hpp b/src/FDM/YASim/Integrator.hpp
index 323177ced..6a84aedf7 100644
--- a/src/FDM/YASim/Integrator.hpp
+++ b/src/FDM/YASim/Integrator.hpp
@@ -17,24 +17,24 @@ class Integrator
 {
 public:
     // Sets the RigidBody that will be integrated.
-    void setBody(RigidBody* body);
+    void setBody(RigidBody* body) { _body = body; }  
 
     // Sets the BodyEnvironment object used to calculate the second
     // derivatives.
-    void setEnvironment(BodyEnvironment* env);
+    void setEnvironment(BodyEnvironment* env) { _env = env; }
 
     // Sets the time interval between steps.  Units can be anything,
     // but they must match the other values (if you specify speed in
     // m/s, then angular acceleration had better be in 1/s^2 and the
     // time interval should be in seconds, etc...)
-    void setInterval(float dt);
-    float getInterval();
+    void setInterval(float dt) { _dt = dt; }
+    float getInterval() const { return _dt; }
 
     // The current state, i.e. initial conditions for the next
     // integration iteration.  Note that the acceleration parameters
     // in the State object are ignored.
-    State* getState();
-    void setState(State* s);
+    State* getState() {  return &_s; }
+    void setState(State* s) {  _s = *s; }
 
     // Do a 4th order Runge-Kutta integration over one time interval.
     // This is the top level of the simulation.
@@ -43,9 +43,13 @@ public:
 private:
     void orthonormalize(float* m);
     void rotMatrix(float* r, float dt, float* out);
-    void l2gVector(float* orient, float* v, float* out);
-    void extrapolatePosition(double* pos, float* v, float dt,
-                             float* o1, float* o2);
+    
+    // Transforms a "local" vector to a "global" vector (not coordinate!)
+    // using the specified orientation.    
+    void l2gVector(const float* orient, const float* v, float* out) { 
+      Math::tmul33(orient, v, out); }
+    
+    void extrapolatePosition(double* pos, float* v, float dt, float* o1, float* o2);
 
     BodyEnvironment* _env;
     RigidBody* _body;
diff --git a/src/FDM/YASim/Jet.cpp b/src/FDM/YASim/Jet.cpp
index ca69b88a5..4d66b5db6 100644
--- a/src/FDM/YASim/Jet.cpp
+++ b/src/FDM/YASim/Jet.cpp
@@ -59,21 +59,6 @@ void Jet::setMaxThrust(float thrust, float afterburner)
     }
 }
 
-void Jet::setVMax(float spd)
-{
-    _vMax = spd;
-}
-
-void Jet::setTSFC(float tsfc)
-{
-    _tsfc = tsfc;
-}
-
-void Jet::setATSFC(float atsfc)
-{
-    _atsfc = atsfc;
-}
-
 void Jet::setRPMs(float idleN1, float maxN1, float idleN2, float maxN2)
 {
     _n1Min = idleN1;
@@ -82,16 +67,6 @@ void Jet::setRPMs(float idleN1, float maxN1, float idleN2, float maxN2)
     _n2Max = maxN2;
 }
 
-void Jet::setEGT(float takeoffEGT)
-{
-    _egt0 = takeoffEGT;
-}
-
-void Jet::setEPR(float takeoffEPR)
-{
-    _epr0 = takeoffEPR;
-}
-
 void Jet::setSpooling(float time)
 {
     // 2.3 = -ln(0.1), i.e. x=2.3 is the 90% point we're defining
@@ -100,11 +75,6 @@ void Jet::setSpooling(float time)
     _decay = 1.5f * 2.3f / time;
 }
 
-void Jet::setVectorAngle(float angle)
-{
-    _maxRot = angle;
-}
-
 void Jet::setReheat(float reheat)
 {
     _reheat = Math::clamp(reheat, 0, 1);
@@ -112,9 +82,7 @@ void Jet::setReheat(float reheat)
 
 void Jet::setRotation(float rot)
 {
-    if(rot < 0) rot = 0;
-    if(rot > 1) rot = 1;
-    _rotControl = rot;
+    _rotControl = Math::clamp(rot, 0, 1);
 }
 
 float Jet::getN1()
@@ -127,16 +95,12 @@ float Jet::getN2()
     return _n2 * _tempCorrect;
 }
 
-float Jet::getEPR()
-{
-    return _epr;
-}
-
 float Jet::getEGT()
 {
     // Exactly zero means "off" -- return the ambient temperature
-    if(_egt == 0) return _temp;
-
+    if(_egt == 0) {
+      return _atmo.getTemperature();
+    }
     return _egt * _tempCorrect * _tempCorrect;
 }
 
@@ -155,8 +119,7 @@ void Jet::integrate(float dt)
     float spd = -Math::dot3(_wind, _dir);
 
     float statT, statP, statD;
-    Atmosphere::calcStaticAir(_pressure, _temp, _rho, spd,
-                              &statP, &statT, &statD);
+    _atmo.calcStaticAir(spd, &statP, &statT, &statD);
     _pressureCorrect = statP/P0;
     _tempCorrect = Math::sqrt(statT/T0);
 
@@ -175,8 +138,8 @@ void Jet::integrate(float dt)
     // Now get a "beta" (i.e. EPR - 1) value.  The output values are
     // expressed as functions of beta.
     float ibeta0 = 1/(_epr0 - 1);
-    float betaTarget = (_epr0 - 1) * (setThrust/_maxThrust) * (P0/_pressure)
-	* (_temp/statT);
+    float betaTarget = (_epr0 - 1) * (setThrust/_maxThrust) * (P0/_atmo.getPressure())
+	* (_atmo.getTemperature()/statT);
     float n2Target = _n2Min + (betaTarget*ibeta0) * (_n2Max - _n2Min);
 
     // Note that this "first" beta value is used to compute a target
@@ -191,7 +154,7 @@ void Jet::integrate(float dt)
     // The actual thrust produced is keyed to the N1 speed.  Add the
     // afterburners in at the end.
     float betaN1 =  (_epr0-1) * (_n1 - _n1Min) / (_n1Max - _n1Min);
-    _thrust = _maxThrust * betaN1/((_epr0-1)*(P0/_pressure)*(_temp/statT));
+    _thrust = _maxThrust * betaN1/((_epr0-1)*(P0/_atmo.getPressure())*(_atmo.getTemperature()/statT));
     _thrust *= 1 + _reheat*(_abFactor-1);
 
     // Finally, calculate the output variables.   Use a 80/20 mix of
@@ -214,16 +177,6 @@ void Jet::integrate(float dt)
     if(_reverseThrust) _thrust *= -_reverseEff;
 }
 
-bool Jet::isRunning()
-{
-    return _running;
-}
-
-bool Jet::isCranking()
-{
-    return _cranking;
-}
-
 void Jet::getThrust(float* out)
 {
     Math::mul3(_thrust, _dir, out);
@@ -240,13 +193,11 @@ void Jet::getThrust(float* out)
 void Jet::getTorque(float* out)
 {
     out[0] = out[1] = out[2] = 0;
-    return;
 }
 
 void Jet::getGyro(float* out)
 {
     out[0] = out[1] = out[2] = 0;
-    return;
 }
 
 }; // namespace yasim
diff --git a/src/FDM/YASim/Jet.hpp b/src/FDM/YASim/Jet.hpp
index 227407223..938367371 100644
--- a/src/FDM/YASim/Jet.hpp
+++ b/src/FDM/YASim/Jet.hpp
@@ -12,13 +12,13 @@ public:
     virtual Jet* getJet() { return this; }
     
     void setMaxThrust(float thrust, float afterburner=0);
-    void setVMax(float spd);
-    void setTSFC(float tsfc);
-    void setATSFC(float atsfc);
+    void setVMax(float spd) { _vMax = spd; };
+    void setTSFC(float tsfc) { _tsfc = tsfc; };
+    void setATSFC(float atsfc) { _atsfc = atsfc; };
     void setRPMs(float idleN1, float maxN1, float idleN2, float maxN2);
-    void setEGT(float takeoffEGT);
-    void setEPR(float takeoffEPR);
-    void setVectorAngle(float angle);
+    void setEGT(float takeoffEGT) { _egt0 = takeoffEGT; };
+    void setEPR(float takeoffEPR) { _epr0 = takeoffEPR; };
+    void setVectorAngle(float angle) { _maxRot = angle; };
 
     // The time it takes the engine to reach 90% thrust from idle
     void setSpooling(float time);
@@ -37,15 +37,15 @@ public:
 
     float getN1();
     float getN2();
-    float getEPR();
+    float getEPR() const { return _epr; };
     float getEGT();
 
     // Normalized "performance" number.  Used for fuzzy numbers in FGFDM
     float getPerfNorm() { return (_n1 - _n1Min) / (_n1Max - _n1Min); }
 
     // From Thruster:
-    virtual bool isRunning();
-    virtual bool isCranking();
+    virtual bool isRunning() { return _running; };
+    virtual bool isCranking() { return _cranking; };
     virtual void getThrust(float* out);
     virtual void getTorque(float* out);
     virtual void getGyro(float* out);
diff --git a/src/FDM/YASim/Launchbar.cpp b/src/FDM/YASim/Launchbar.cpp
index 26d2ab58e..8a1a276fe 100644
--- a/src/FDM/YASim/Launchbar.cpp
+++ b/src/FDM/YASim/Launchbar.cpp
@@ -374,8 +374,8 @@ void Launchbar::calcForce(Ground *g_cb, RigidBody* body, State* s, float* lv, fl
     // Transform the velocities of the endpoints to the
     // local coordinate sytem.
     float lvel[2][3];
-    s->velGlobalToLocal(vel[0], lvel[0]);
-    s->velGlobalToLocal(vel[1], lvel[1]);
+    s->globalToLocal(vel[0], lvel[0]);
+    s->globalToLocal(vel[1], lvel[1]);
 
     // Compute the position of the launchbar tip relative to the cat.
     float tip_pos_on_cat = getPercentPosOnCat(llb_mount, 0.0, lend);
diff --git a/src/FDM/YASim/Model.cpp b/src/FDM/YASim/Model.cpp
index 7a9da1519..4ca6d95da 100644
--- a/src/FDM/YASim/Model.cpp
+++ b/src/FDM/YASim/Model.cpp
@@ -50,28 +50,13 @@ void printState(State* s)
 
 Model::Model()
 {
-    int i;
-    for(i=0; i<3; i++) _wind[i] = 0;
-
     _integrator.setBody(&_body);
     _integrator.setEnvironment(this);
-
     // Default value of 30 Hz
     _integrator.setInterval(1.0f/30.0f);
 
-    _agl = 0;
-    _crashed = false;
-    _turb = 0;
     _ground_cb = new Ground();
-    _hook = 0;
-    _launchbar = 0;
 
-    _wingSpan = 0;
-    _groundEffect = 0;
-    for(i=0; i<3; i++) _geRefPoint[i] = 0;
-
-    _global_ground[0] = 0; _global_ground[1] = 0; _global_ground[2] = 1;
-    _global_ground[3] = -100000;
     _modelN = fgGetNode("/fdm/yasim/forces", true);
     _fAeroXN = _modelN->getNode("f-x-drag", true);
     _fAeroYN = _modelN->getNode("f-y-side", true);
@@ -98,11 +83,10 @@ void Model::getThrust(float* out) const
 {
     float tmp[3];
     out[0] = out[1] = out[2] = 0;
-    int i;
-    for(i=0; i<_thrusters.size(); i++) {
-	Thruster* t = (Thruster*)_thrusters.get(i);
-	t->getThrust(tmp);
-	Math::add3(tmp, out, out);
+    for(int i=0; i<_thrusters.size(); i++) {
+        Thruster* t = (Thruster*)_thrusters.get(i);
+        t->getThrust(tmp);
+        Math::add3(tmp, out, out);
     }
 }
 
@@ -127,7 +111,7 @@ void Model::initIteration()
         localWind(pos, _s, v, alt);
 
 	t->setWind(v);
-	t->setAir(_pressure, _temp, _rho);
+	t->setAtmosphere(_atmo);
 	t->integrate(_integrator.getInterval());
 
 	t->getTorque(v);
@@ -153,8 +137,6 @@ void Model::initIteration()
         Hitch* h = (Hitch*)_hitches.get(i);
         h->integrate(_integrator.getInterval());
     }
-
-    
 }
 
 // This function initializes some variables for the rotor calculation
@@ -186,7 +168,7 @@ void Model::iterate()
 void Model::setState(State* s)
 {
     _integrator.setState(s);
-    _s = _integrator.getState();
+    _s = s;
 }
 
 
@@ -203,20 +185,6 @@ void Model::setGroundEffect(const float* pos, float span, float mul)
     _groundEffect = mul;
 }
 
-void Model::setAir(float pressure, float temp, float density)
-{
-    _pressure = pressure;
-    _temp = temp;
-    _rho = density;
-}
-
-void Model::setAirFromStandardAtmosphere(float altitude)
-{
-    _pressure = Atmosphere::getStdPressure(altitude);
-    _temp = Atmosphere::getStdTemperature(altitude);
-    _rho = Atmosphere::getStdDensity(altitude);
-}
-
 void Model::updateGround(State* s)
 {
     float dummy[3];
@@ -336,7 +304,7 @@ void Model::calcForces(State* s)
       localWind(pos, s, vs, alt);
 
       float force[3], torque[3];
-      sf->calcForce(vs, _rho, force, torque);
+      sf->calcForce(vs, _atmo.getDensity(), force, torque);
       Math::add3(faero, force, faero);
 
       _body.addForce(pos, force);
@@ -349,7 +317,7 @@ void Model::calcForces(State* s)
         float vs[3], pos[3];
         r->getPosition(pos);
         localWind(pos, s, vs, alt);
-        r->calcLiftFactor(vs, _rho,s);
+        r->calcLiftFactor(vs, _atmo.getDensity(), s);
         float tq=0; 
         // total torque of rotor (scalar) for calculating new rotor rpm
 
@@ -363,7 +331,7 @@ void Model::calcForces(State* s)
             localWind(pos, s, vs, alt,true);
 
             float force[3], torque[3];
-            rp->calcForce(vs, _rho, force, torque, &torque_scalar);
+            rp->calcForce(vs, _atmo.getDensity(), force, torque, &torque_scalar);
             tq+=torque_scalar;
             rp->getPositionForceAttac(pos);
 
@@ -486,7 +454,7 @@ void Model::newState(State* s)
 
 // Calculates the airflow direction at the given point and for the
 // specified aircraft velocity.
-void Model::localWind(const float* pos, yasim::State* s, float* out, float alt, bool is_rotor)
+void Model::localWind(const float* pos, const yasim::State* s, float* out, float alt, bool is_rotor)
 {
     float tmp[3], lwind[3], lrot[3], lv[3];
 
@@ -521,8 +489,6 @@ void Model::localWind(const float* pos, yasim::State* s, float* out, float alt,
         _rotorgear.getDownWash(pos,lv,tmp);
         Math::add3(out,tmp, out);    //  + downwash
     }
-
-
 }
 
 }; // namespace yasim
diff --git a/src/FDM/YASim/Model.hpp b/src/FDM/YASim/Model.hpp
index 5ea3b358a..75e2de16e 100644
--- a/src/FDM/YASim/Model.hpp
+++ b/src/FDM/YASim/Model.hpp
@@ -7,6 +7,7 @@
 #include "Vector.hpp"
 #include "Turbulence.hpp"
 #include "Rotor.hpp"
+#include "Atmosphere.hpp"
 #include <simgear/props/props.hxx>
 
 namespace yasim {
@@ -68,8 +69,8 @@ public:
     //
     void setGroundEffect(const float* pos, float span, float mul);
     void setWind(float* wind) { Math::set3(wind, _wind); }
-    void setAir(float pressure, float temp, float density);
-    void setAirFromStandardAtmosphere(float altitude);
+    void setAtmosphere(Atmosphere a) { _atmo = a; };
+    void setStandardAtmosphere(float altitude) { _atmo.setStandard(altitude); };
 
     void updateGround(State* s);
 
@@ -81,39 +82,38 @@ private:
     void initRotorIteration();
     void calcGearForce(Gear* g, float* v, float* rot, float* ground);
     float gearFriction(float wgt, float v, Gear* g);
-    void localWind(const float* pos, State* s, float* out, float alt, bool is_rotor = false);
+    void localWind(const float* pos, const yasim::State* s, float* out, float alt, bool is_rotor = false);
 
     Integrator _integrator;
     RigidBody _body;
 
-    Turbulence* _turb;
+    Turbulence* _turb {nullptr};
 
     Vector _thrusters;
     Vector _surfaces;
     Rotorgear _rotorgear;
     Vector _gears;
-    Hook* _hook;
-    Launchbar* _launchbar;
+    Hook* _hook {nullptr};
+    Launchbar* _launchbar {nullptr};
     Vector _hitches;
 
-    float _wingSpan;
-    float _groundEffect;
-    float _geRefPoint[3];
+    float _wingSpan {0};
+    float _groundEffect {0};
+    float _geRefPoint[3] {0,0,0};
 
     Ground* _ground_cb;
-    double _global_ground[4];
-    float _pressure;
-    float _temp;
-    float _rho;
-    float _wind[3];
+    double _global_ground[4] {0,0,1, -1e5};
+    Atmosphere _atmo;
+    float _wind[3] {0,0,0};
+    
 
     // Accumulators for the total internal gyro and engine torque
-    float _gyro[3];
-    float _torque[3];
+    float _gyro[3] {0,0,0};
+    float _torque[3] {0,0,0};
 
     State* _s;
-    bool _crashed;
-    float _agl;
+    bool _crashed {false};
+    float _agl {0};
     SGPropertyNode_ptr _modelN;  
     SGPropertyNode_ptr _fAeroXN;
     SGPropertyNode_ptr _fAeroYN;
diff --git a/src/FDM/YASim/PropEngine.cpp b/src/FDM/YASim/PropEngine.cpp
index c6581f445..cb357e9fc 100644
--- a/src/FDM/YASim/PropEngine.cpp
+++ b/src/FDM/YASim/PropEngine.cpp
@@ -123,8 +123,8 @@ void PropEngine::stabilize()
     // If we cannot manage this in 100 iterations, give up.
     for (int n = 0; n < 100; n++) {
 	float ptau, thrust;
-	_prop->calc(_rho, speed, _omega * _gearRatio, &thrust, &ptau);
-	_eng->calc(_pressure, _temp, _omega);
+	_prop->calc(_atmo.getDensity(), speed, _omega * _gearRatio, &thrust, &ptau);
+	_eng->calc(_atmo.getPressure(), _atmo.getTemperature(), _omega);
         _eng->stabilize();
 
         // Do it again -- the turbo sets the target MP in the first
@@ -133,7 +133,7 @@ void PropEngine::stabilize()
         // it works without side effects (other than solver
         // performance).  In the future, the Engine objects should
         // store state to allow them to do the work themselves.
-	_eng->calc(_pressure, _temp, _omega);
+	_eng->calc(_atmo.getPressure(), _atmo.getTemperature(), _omega);
 
         // Compute torque as seen by the engine's end of the gearbox.
         // The propeller will be moving more slowly (for gear ratios
@@ -185,11 +185,11 @@ void PropEngine::integrate(float dt)
     _eng->setMixture(_mixture);
     _eng->setFuelState(_fuel);
     
-    _prop->calc(_rho, speed, _omega * _gearRatio, &thrust, &propTorque);
+    _prop->calc(_atmo.getDensity(), speed, _omega * _gearRatio, &thrust, &propTorque);
     if(_omega == 0.0)
         _omega = 0.001; // hack to get around reports of NaNs somewhere...
     propTorque *= _gearRatio;
-    _eng->calc(_pressure, _temp, _omega);
+    _eng->calc(_atmo.getPressure(), _atmo.getTemperature(), _omega);
     _eng->integrate(dt);
     engTorque = _eng->getTorque();
     _fuelFlow = _eng->getFuelFlow();
diff --git a/src/FDM/YASim/Rotor.cpp b/src/FDM/YASim/Rotor.cpp
index 14a0e3edc..f4c338562 100644
--- a/src/FDM/YASim/Rotor.cpp
+++ b/src/FDM/YASim/Rotor.cpp
@@ -829,7 +829,7 @@ void Rotor::calcLiftFactor(float* v, float rho, State *s)
 
     //store the gravity direction
     Glue::geodUp(s->pos, _grav_direction);
-    s->velGlobalToLocal(_grav_direction, _grav_direction);
+    s->globalToLocal(_grav_direction, _grav_direction);
 }
 
 void Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s)
diff --git a/src/FDM/YASim/Surface.cpp b/src/FDM/YASim/Surface.cpp
index 475c9ec90..6f01dc836 100644
--- a/src/FDM/YASim/Surface.cpp
+++ b/src/FDM/YASim/Surface.cpp
@@ -7,35 +7,12 @@ int Surface::s_idGenerator = 0;
 Surface::Surface( Version * version ) :
     _version(version)
 {
-    // create id for surface
     _id = s_idGenerator++;
-    // Start in a "sane" mode, so unset stuff doesn't freak us out
-    _c0 = 1;
-    _cx = _cy = _cz = 1;
-    _cz0 = 0;
-    _peaks[0] = _peaks[1] = 1;
-    int i;
-    for(i=0; i<4; i++) {
-	_stalls[i] = 0;
-        _widths[i] = 0.01; // half a degree
-    }
+
     _orient[0] = 1; _orient[1] = 0; _orient[2] = 0;
     _orient[3] = 0; _orient[4] = 1; _orient[5] = 0;
     _orient[6] = 0; _orient[7] = 0; _orient[8] = 1;
     
-    _chord = 0;
-    _incidence = 0;
-    _twist = 0;
-    _slatPos = _spoilerPos = _flapPos = 0;
-    _slatDrag = _spoilerDrag = _flapDrag = 1;
-
-    _flapLift = 0;
-    _flapEffectiveness = 1;
-    _slatAlpha = 0;
-    _spoilerLift = 1;
-    _inducedDrag = 1;
-    _stallAlpha = 0;
-    _alpha = 0;
     _surfN = fgGetNode("/fdm/yasim/debug/surfaces", true);
     if (_surfN != 0) {
       _surfN = _surfN->getChild("surface", _id, true);
diff --git a/src/FDM/YASim/Surface.hpp b/src/FDM/YASim/Surface.hpp
index f65890cfe..6114dc26f 100644
--- a/src/FDM/YASim/Surface.hpp
+++ b/src/FDM/YASim/Surface.hpp
@@ -92,36 +92,36 @@ private:
     float flapLift(float alpha);
     float controlDrag(float lift, float drag);
 
-    float _chord;     // X-axis size
-    float _c0;        // total force coefficient
-    float _cx;        // X-axis force coefficient
-    float _cy;        // Y-axis force coefficient
-    float _cz;        // Z-axis force coefficient
-    float _cz0;       // Z-axis force offset
-    float _peaks[2];  // Stall peak coefficients (fwd, back)
-    float _stalls[4]; // Stall angles (fwd/back, pos/neg)
-    float _widths[4]; // Stall widths  " "
+    float _chord {0};     // X-axis size
+    float _c0 {1};        // total force coefficient
+    float _cx {1};        // X-axis force coefficient
+    float _cy {1};        // Y-axis force coefficient
+    float _cz {1};        // Z-axis force coefficient
+    float _cz0 {0};       // Z-axis force offset
+    float _peaks[2] {1, 1};  // Stall peak coefficients (fwd, back)
+    float _stalls[4] {0, 0, 0, 0}; // Stall angles (fwd/back, pos/neg)
+    float _widths[4] {0.01, 0.01, 0.01, 0.01}; // Stall widths  " "
     float _pos[3];    // position in local coords
     float _orient[9]; // local->surface orthonormal matrix
 
-    float _slatAlpha;
-    float _slatDrag;
-    float _flapLift;
-    float _flapDrag;
-    float _flapEffectiveness;
-    float _spoilerLift;
-    float _spoilerDrag;
+    float _slatAlpha {0};
+    float _slatDrag {1};
+    float _flapLift {0};
+    float _flapDrag {1};
+    float _flapEffectiveness {1};
+    float _spoilerLift {1};
+    float _spoilerDrag {1};
 
-    float _slatPos;
-    float _flapPos;
-    float _spoilerPos;
-    float _incidence;
-    float _twist;
-    float _inducedDrag;
+    float _slatPos {0};
+    float _flapPos {0};
+    float _spoilerPos {0};
+    float _incidence {0};
+    float _twist {0};
+    float _inducedDrag {1};
     
     // used during calculations
-    float _stallAlpha;
-    float _alpha;
+    float _stallAlpha {0};
+    float _alpha {0};
 
     Version * _version;
     SGPropertyNode* _fxN;
diff --git a/src/FDM/YASim/Thruster.cpp b/src/FDM/YASim/Thruster.cpp
deleted file mode 100644
index cfc31a7a9..000000000
--- a/src/FDM/YASim/Thruster.cpp
+++ /dev/null
@@ -1,72 +0,0 @@
-#include "Math.hpp"
-#include "Thruster.hpp"
-namespace yasim {
-
-Thruster::Thruster()
-{
-    _dir[0] = 1; _dir[1] = 0; _dir[2] = 0;
-    int i;
-    for(i=0; i<3; i++) _pos[i] = _wind[i] = 0;
-    _throttle = 0;
-    _mixture = 0;
-    _starter = false;
-    _pressure = _temp = _rho = 0;
-}
-
-Thruster::~Thruster()
-{
-}
-
-void Thruster::getPosition(float* out)
-{
-    int i;
-    for(i=0; i<3; i++) out[i] = _pos[i];
-}
-
-void Thruster::setPosition(float* pos)
-{
-    int i;
-    for(i=0; i<3; i++) _pos[i] = pos[i];
-}
-
-void Thruster::getDirection(float* out)
-{
-    int i;
-    for(i=0; i<3; i++) out[i] = _dir[i];
-}
-
-void Thruster::setDirection(float* dir)
-{
-    Math::unit3(dir, _dir);
-}
-
-void Thruster::setThrottle(float throttle)
-{
-    _throttle = Math::clamp(throttle, -1, 1);
-}
-
-void Thruster::setMixture(float mixture)
-{
-    _mixture = Math::clamp(mixture, 0, 1);
-}
-
-
-void Thruster::setStarter(bool starter)
-{
-    _starter = starter;
-}
-
-void Thruster::setWind(float* wind)
-{
-    int i;
-    for(i=0; i<3; i++) _wind[i] = wind[i];
-}
-
-void Thruster::setAir(float pressure, float temp, float density)
-{
-    _pressure = pressure;
-    _temp = temp;
-    _rho = density;
-}
-
-}; // namespace yasim
diff --git a/src/FDM/YASim/Thruster.hpp b/src/FDM/YASim/Thruster.hpp
index c022ebaa4..210e88342 100644
--- a/src/FDM/YASim/Thruster.hpp
+++ b/src/FDM/YASim/Thruster.hpp
@@ -1,6 +1,9 @@
 #ifndef _THRUSTER_HPP
 #define _THRUSTER_HPP
 
+#include "Atmosphere.hpp"
+#include "Math.hpp"
+
 namespace yasim {
 
 class Jet;
@@ -10,8 +13,8 @@ class Engine;
 
 class Thruster {
 public:
-    Thruster();
-    virtual ~Thruster();
+    Thruster() {};
+    virtual ~Thruster() {};
 
     // Typing info, these are the possible sub-type (or sub-parts)
     // that a thruster might have.  Any might return null.  A little
@@ -22,15 +25,15 @@ public:
     virtual Engine* getEngine() { return 0; }
     
     // Static data
-    void getPosition(float* out);
-    void setPosition(float* pos);
-    void getDirection(float* out);
-    void setDirection(float* dir);
+    void getPosition(float* out) const { Math::set3(_pos, out); }
+    void setPosition(const float* pos) { Math::set3(pos, _pos); }
+    void getDirection(float* out) const { Math::set3(_dir, out); }
+    void setDirection(const float* dir) { Math::unit3(dir, _dir); }
 
     // Controls
-    void setThrottle(float throttle);
-    void setMixture(float mixture);
-    void setStarter(bool starter);
+    void setThrottle(float throttle) { _throttle = Math::clamp(throttle, -1, 1); }
+    void setMixture(float mixture) { _mixture = Math::clamp(mixture, 0, 1); }
+    void setStarter(bool starter) { _starter = starter; }
     void setFuelState(bool hasFuel) { _fuel = hasFuel; }
 
     // Dynamic output
@@ -42,26 +45,24 @@ public:
     virtual float getFuelFlow()=0; // in kg/s
 
     // Runtime instructions
-    void setWind(float* wind);
-    void setAir(float pressure, float temp, float density);
+    void setWind(const float* wind) { Math::set3(wind, _wind); };
+    void setAtmosphere(Atmosphere a) { _atmo = a; };
+    void setStandardAtmosphere(float altitude) { _atmo.setStandard(altitude); };
     virtual void init() {}
     virtual void integrate(float dt)=0;
     virtual void stabilize()=0;
 
 protected:
-    float _pos[3];
-    float _dir[3];
-    float _throttle;
-    float _mixture;
-    bool _starter; // true=engaged, false=disengaged
+    float _pos[3] {0, 0, 0};
+    float _dir[3] {1, 0, 0};
+    float _throttle = 0;
+    float _mixture = 0;
+    bool _starter = false; // true=engaged, false=disengaged
     bool _fuel; // true=available, false=out
 
-    float _wind[3];
-    float _pressure;
-    float _temp;
-    float _rho;
+    float _wind[3] {0, 0, 0};
+    Atmosphere _atmo;
 };
 
 }; // namespace yasim
 #endif // _THRUSTER_HPP
-
diff --git a/src/FDM/YASim/Version.hpp b/src/FDM/YASim/Version.hpp
index 0c56f06a6..10e7e3e01 100644
--- a/src/FDM/YASim/Version.hpp
+++ b/src/FDM/YASim/Version.hpp
@@ -16,20 +16,20 @@ public:
   } YASIM_VERSION;
 
   void setVersion( const char * version );
-  int getVersion() { return _version; }
-  bool isVersion( YASIM_VERSION version );
-  bool isVersionOrNewer( YASIM_VERSION version );
+  int getVersion() const { return _version; }
+  bool isVersion( YASIM_VERSION version ) const;
+  bool isVersionOrNewer( YASIM_VERSION version ) const;
 
 private:
   YASIM_VERSION _version;
 };
 
-inline bool Version::isVersion( YASIM_VERSION version )
+inline bool Version::isVersion( YASIM_VERSION version ) const
 {
   return _version == version;
 }
 
-inline bool Version::isVersionOrNewer( YASIM_VERSION version )
+inline bool Version::isVersionOrNewer( YASIM_VERSION version ) const
 {
   return _version >= version;
 }
diff --git a/src/FDM/YASim/Wing.cpp b/src/FDM/YASim/Wing.cpp
index 68489231a..d060561df 100644
--- a/src/FDM/YASim/Wing.cpp
+++ b/src/FDM/YASim/Wing.cpp
@@ -7,47 +7,12 @@ static const float RAD2DEG = 57.2957795131;
 Wing::Wing( Version * version ) :
   _version(version)
 {
-    _mirror = false;
-    _base[0] = _base[1] = _base[2] = 0;
-    _length = 0;
-    _chord = 0;
-    _taper = 0;
-    _sweep = 0;
-    _dihedral = 0;
-    _stall = 0;
-    _stallWidth = 0;
-    _stallPeak = 0;
-    _twist = 0;
-    _camber = 0;
-    _incidence = 0;
-    _inducedDrag = 1;
-    _dragScale = 1;
-    _liftRatio = 1;
-    _flap0Start = 0;
-    _flap0End = 0;
-    _flap0Lift = 0;
-    _flap0Drag = 0;
-    _flap1Start = 0;
-    _flap1End = 0;
-    _flap1Lift = 0;
-    _flap1Drag = 0;
-    _spoilerStart = 0;
-    _spoilerEnd = 0;
-    _spoilerLift = 0;
-    _spoilerDrag = 0;
-    _slatStart = 0;
-    _slatEnd = 0;
-    _slatAoA = 0;
-    _slatDrag = 0;
-    _meanChord = 0;
-    _wingspan = 0;
-    _aspectRatio = 1;
+
 }
 
 Wing::~Wing()
 {
-    int i;
-    for(i=0; i<_surfs.size(); i++) {
+    for(int i=0; i<_surfs.size(); i++) {
         SurfRec* s = (SurfRec*)_surfs.get(i);
         delete s->surface;
         delete s;
@@ -57,8 +22,7 @@ Wing::~Wing()
 void Wing::setIncidence(float incidence)
 {
     _incidence = incidence;
-    int i;
-    for(i=0; i<_surfs.size(); i++)
+    for(int i=0; i<_surfs.size(); i++)
         ((SurfRec*)_surfs.get(i))->surface->setIncidence(incidence);
 }
 
diff --git a/src/FDM/YASim/Wing.hpp b/src/FDM/YASim/Wing.hpp
index df289f2f6..bdaf128a1 100644
--- a/src/FDM/YASim/Wing.hpp
+++ b/src/FDM/YASim/Wing.hpp
@@ -106,55 +106,55 @@ private:
     Vector _slatSurfs;
     Vector _spoilerSurfs;
 
-    bool _mirror;
+    bool _mirror {false};
 
-    float _base[3];
-    float _length;
-    float _chord;
-    float _taper;
-    float _sweep;
-    float _dihedral;
+    float _base[3] {0,0,0};
+    float _length {0};
+    float _chord {0};
+    float _taper {1};
+    float _sweep {0};
+    float _dihedral {0};
     
     // calculated from above
-    float _tip[3];
-    float _meanChord; // std. mean chord
-    float _mac; // mean aerodynamic chord length
-    float _macRootDistance; // y-distance of mac from root
-    float _macX; // x-coordinate of mac (leading edge)
-    float _netSpan;
-    float _wingspan;
-    float _aspectRatio;
+    float _tip[3] {0,0,0};
+    float _meanChord {0}; // std. mean chord
+    float _mac {0}; // mean aerodynamic chord length
+    float _macRootDistance {0}; // y-distance of mac from root
+    float _macX {0}; // x-coordinate of mac (leading edge)
+    float _netSpan {0};
+    float _wingspan {0};
+    float _aspectRatio {1};
 
-    float _stall;
-    float _stallWidth;
-    float _stallPeak;
-    float _twist;
-    float _camber;
-    float _incidence;
-    float _inducedDrag;
+    float _stall {0};
+    float _stallWidth {0};
+    float _stallPeak {0};
+    float _twist {0};
+    float _camber {0};
+    float _incidence {0};
+    float _inducedDrag {1};
 
-    float _dragScale;
-    float _liftRatio;
+    float _dragScale {1};
+    float _liftRatio {1};
 
-    float _flap0Start;
-    float _flap0End;
-    float _flap0Lift;
-    float _flap0Drag;
+    float _flap0Start {0};
+    float _flap0End {0};
+    float _flap0Lift {0};
+    float _flap0Drag {0};
 
-    float _flap1Start;
-    float _flap1End;
-    float _flap1Lift;
-    float _flap1Drag;
+    float _flap1Start {0};
+    float _flap1End {0};
+    float _flap1Lift {0};
+    float _flap1Drag {0};
 
-    float _spoilerStart;
-    float _spoilerEnd;
-    float _spoilerLift;
-    float _spoilerDrag;
+    float _spoilerStart {0};
+    float _spoilerEnd {0};
+    float _spoilerLift {0};
+    float _spoilerDrag {0};
 
-    float _slatStart;
-    float _slatEnd;
-    float _slatAoA;
-    float _slatDrag;
+    float _slatStart {0};
+    float _slatEnd {0};
+    float _slatAoA {0};
+    float _slatDrag {0};
 
     Version * _version;
 };
diff --git a/src/FDM/YASim/YASim.cxx b/src/FDM/YASim/YASim.cxx
index 868f65f88..cfdf70d3f 100644
--- a/src/FDM/YASim/YASim.cxx
+++ b/src/FDM/YASim/YASim.cxx
@@ -15,6 +15,7 @@
 #include <Main/globals.hxx>
 #include <Main/fg_props.hxx>
 
+#include "yasim-common.hpp"
 #include "FGFDM.hpp"
 #include "Atmosphere.hpp"
 #include "Math.hpp"
@@ -34,22 +35,6 @@
 using namespace yasim;
 using std::string;
 
-static const float YASIM_PI = 3.14159265358979323846;
-static const float RAD2DEG = 180/YASIM_PI;
-static const float PI2 = YASIM_PI*2;
-static const float M2FT = 3.2808399;
-static const float FT2M = 0.3048;
-static const float MPS2KTS = 3600.0/1852.0;
-static const float INHG2PA = 3386.389;
-static const float SLUG2KG = 14.59390;
-
-#if 0 // unused constant values
-static const float RAD2RPM = 9.54929658551;
-static const float CM2GALS = 264.172037284; // gallons/cubic meter
-static const float KG2LBS = 2.20462262185;
-static const float W2HP = 1.3416e-3;
-#endif
-
 YASim::YASim(double dt) :
     _simTime(0)
 {
@@ -312,7 +297,11 @@ void YASim::copyToYASim(bool copyState)
     float temp     = _temp_degc->getFloatValue() + 273.15;
     float dens     = _density_slugft3->getFloatValue() *
                             SLUG2KG * M2FT*M2FT*M2FT;
-
+    Atmosphere atmo;
+    atmo.setDensity(dens);
+    atmo.setTemperature(temp);
+    atmo.setPressure(pressure);
+                            
     // Convert and set:
     Model* model = _fdm->getAirplane()->getModel();
     State s;
@@ -332,7 +321,7 @@ void YASim::copyToYASim(bool copyState)
     Math::mmul33(s.orient, xyz2ned, s.orient);
 
     // Velocity
-    float v[3];
+    float v[3] {0, 0, 0};
     bool needCopy = false;
     switch (_speed_set) {
     case NED:
@@ -347,15 +336,14 @@ void YASim::copyToYASim(bool copyState)
         Math::tmul33(s.orient, v, v);
         break;
     case KNOTS:
-        v[0] = Atmosphere::spdFromVCAS(get_V_calibrated_kts()/MPS2KTS,
-                                       pressure, temp);
+        v[0] = atmo.spdFromVCAS(get_V_calibrated_kts()/MPS2KTS);
         v[1] = 0;
         v[2] = 0;
         Math::tmul33(s.orient, v, v);
         needCopy = true;
         break;
     case MACH:
-        v[0] = Atmosphere::spdFromMach(get_Mach_number(), temp);
+        v[0] = atmo.spdFromMach(get_Mach_number());
         v[1] = 0;
         v[2] = 0;
         Math::tmul33(s.orient, v, v);
@@ -374,20 +362,18 @@ void YASim::copyToYASim(bool copyState)
     if(copyState || needCopy)
         model->setState(&s);
 
-    // wind
     Math::tmul33(xyz2ned, wind, wind);
     model->setWind(wind);
-
-    // air
-    model->setAir(pressure, temp, dens);
+    model->setAtmosphere(atmo);
 
     // Query a ground plane for each gear/hook/launchbar and
     // write that value into the corresponding class.
     _fdm->getAirplane()->getModel()->updateGround(&s);
 
     Launchbar* l = model->getLaunchbar();
-    if (l)
+    if (l) {
         l->setLaunchCmd(0.0 < _catapult_launch_cmd->getFloatValue());
+    }
 }
 
 // All the settables:
diff --git a/src/FDM/YASim/proptest.cpp b/src/FDM/YASim/proptest.cpp
index c0482b8fb..c32c849e3 100644
--- a/src/FDM/YASim/proptest.cpp
+++ b/src/FDM/YASim/proptest.cpp
@@ -29,23 +29,6 @@ float fgGetFloat (const char * name, float defaultValue) { return 0; }
 double fgGetDouble (const char * name, double defaultValue) { return 0; }
 bool fgSetDouble (const char * name, double defaultValue) { return false; }
 
-static const float KTS2MPS = 0.514444444444;
-static const float RPM2RAD = 0.10471975512;
-static const float HP2W = 745.700;
-static const float FT2M = 0.3048;
-static const float N2LB = 0.224809;
-
-// static const float DEG2RAD = 0.0174532925199;
-// static const float LBS2N = 4.44822;
-// static const float LBS2KG = 0.45359237;
-// static const float KG2LBS = 2.2046225;
-// static const float CM2GALS = 264.172037284;
-// static const float INHG2PA = 3386.389;
-// static const float K2DEGF = 1.8;
-// static const float K2DEGFOFFSET = -459.4;
-// static const float CIN2CM = 1.6387064e-5;
-// static const float YASIM_PI = 3.14159265358979323846;
-
 const int COUNT = 100;
 
 int main(int argc, char** argv)
@@ -71,9 +54,7 @@ int main(int argc, char** argv)
     eng->setBoost(1);
 
     float alt = (argc > 2 ? atof(argv[2]) : 0) * FT2M;
-    pe->setAir(Atmosphere::getStdPressure(alt),
-               Atmosphere::getStdTemperature(alt),
-               Atmosphere::getStdDensity(alt));
+    pe->setStandardAtmosphere(alt);
  
     float speed = (argc > 3 ? atof(argv[3]) : 0) * KTS2MPS;
     float wind[3];
diff --git a/src/FDM/YASim/yasim-atmotest.cpp b/src/FDM/YASim/yasim-atmotest.cpp
new file mode 100644
index 000000000..384cbd522
--- /dev/null
+++ b/src/FDM/YASim/yasim-atmotest.cpp
@@ -0,0 +1,12 @@
+#include "Atmosphere.hpp"
+
+using namespace yasim;
+
+int main(int argc, char** argv)
+{
+  Atmosphere a;
+  if (a.test()) {
+    return 0;
+  }
+  return 1;
+}
diff --git a/src/FDM/YASim/yasim-common.hpp b/src/FDM/YASim/yasim-common.hpp
new file mode 100644
index 000000000..6f2678087
--- /dev/null
+++ b/src/FDM/YASim/yasim-common.hpp
@@ -0,0 +1,34 @@
+#ifndef _YASIM_COMMON_HPP
+#define _YASIM_COMMON_HPP
+
+namespace yasim {
+  static const float YASIM_PI = 3.14159265358979323846;
+  static const float PI2 = YASIM_PI*2;
+  static const float RAD2DEG = 180/YASIM_PI;
+  static const float DEG2RAD = YASIM_PI/180;
+  static const float RPM2RAD = YASIM_PI/30;
+
+  static const float KTS2MPS = 1852.0/3600.0;
+  static const float MPS2KTS = 3600.0/1852.0;
+  static const float KMH2MPS = 1/3.6;
+
+  static const float FT2M = 0.3048;
+  static const float M2FT = 1/FT2M;
+
+  static const float LBS2N = 4.44822;
+  static const float N2LB = 1/LBS2N;
+  static const float LBS2KG = 0.45359237;
+  static const float KG2LBS = 1/LBS2KG;
+  static const float CM2GALS = 264.172037284;
+  static const float HP2W = 745.700;
+  static const float INHG2PA = 3386.389;
+  static const float K2DEGF = 1.8;
+  static const float K2DEGFOFFSET = -459.4;
+  static const float CIN2CM = 1.6387064e-5;
+  
+  static const float NM2FTLB = (1/(LBS2N*FT2M));
+  static const float SLUG2KG = 14.59390;
+    
+};
+
+#endif // ifndef _YASIM_COMMON_HPP
diff --git a/src/FDM/YASim/yasim-test.cpp b/src/FDM/YASim/yasim-test.cpp
index 9e40dc233..58efbfafe 100644
--- a/src/FDM/YASim/yasim-test.cpp
+++ b/src/FDM/YASim/yasim-test.cpp
@@ -7,6 +7,7 @@
 #include <simgear/xml/easyxml.hxx>
 #include <simgear/misc/sg_path.hxx>
 
+#include "yasim-common.hpp"
 #include "FGFDM.hpp"
 #include "Atmosphere.hpp"
 #include "RigidBody.hpp"
@@ -26,11 +27,6 @@ float fgGetFloat (const char * name, float defaultValue) { return 0; }
 double fgGetDouble (const char * name, double defaultValue = 0.0) { return 0; }
 bool fgSetDouble (const char * name, double defaultValue = 0.0) { return 0; }
 
-static const float RAD2DEG = 57.2957795131;
-static const float DEG2RAD = 0.0174532925199;
-/// knots 2 meters per second
-static const float KTS2MPS = 0.514444444444;
-
 
 enum Config 
 {
@@ -54,7 +50,7 @@ void yasim_graph(Airplane* a, const float alt, const float kts, int cfg = CONFIG
   Model* m = a->getModel();
   State s;
 
-  m->setAirFromStandardAtmosphere(alt);
+  m->setStandardAtmosphere(alt);
 
   switch (cfg) {
     case CONFIG_APPROACH:
@@ -74,7 +70,7 @@ void yasim_graph(Airplane* a, const float alt, const float kts, int cfg = CONFIG
 
   for(int deg=-15; deg<=90; deg++) {
     float aoa = deg * DEG2RAD;
-    Airplane::setupState(aoa, kts * KTS2MPS, 0 ,&s);
+    s.setupState(aoa, kts * KTS2MPS, 0);
     m->getBody()->reset();
     m->initIteration();
     m->calcForces(&s);
@@ -129,7 +125,7 @@ void yasim_drag(Airplane* a, const float aoa, const float alt, int cfg = CONFIG_
   Model* m = a->getModel();
   State s;
   
-  m->setAirFromStandardAtmosphere(alt);
+  m->setStandardAtmosphere(alt);
   
   switch (cfg) {
     case CONFIG_APPROACH:
@@ -148,7 +144,7 @@ void yasim_drag(Airplane* a, const float aoa, const float alt, int cfg = CONFIG_
   printf("#kts, drag\n");
   
   for(int kts=15; kts<=150; kts++) {
-    Airplane::setupState(aoa, kts * KTS2MPS, 0 ,&s);
+    s.setupState(aoa, kts * KTS2MPS, 0);
     m->getBody()->reset();
     m->initIteration();
     m->calcForces(&s);
@@ -254,11 +250,11 @@ int main(int argc, char** argv)
     float MACy = a->getWing()->getMACy();
     
     printf("       Iterations: %d\n", a->getSolutionIterations());
-    printf(" Drag Coefficient: %f\n", drag);
-    printf("       Lift Ratio: %f\n", a->getLiftRatio());
-    printf("       Cruise AoA: %f deg\n", aoa);
-    printf("   Tail Incidence: %f deg\n", tail);
-    printf("Approach Elevator: %f\n\n", a->getApproachElevator());
+    printf(" Drag Coefficient: %.3f\n", drag);
+    printf("       Lift Ratio: %.3f\n", a->getLiftRatio());
+    printf("       Cruise AoA: %.2f deg\n", aoa);
+    printf("   Tail Incidence: %.2f deg\n", tail);
+    printf("Approach Elevator: %.3f\n\n", a->getApproachElevator());
     printf("               CG: x:%.3f, y:%.3f, z:%.3f\n", cg[0], cg[1], cg[2]);
     printf("    Wing MAC (*1): x:%.2f, y:%.2f, length:%.1f \n", MACx, MACy, MAC);
     printf("    CG-x rel. MAC: %.3f\n", a->getCGMAC());