Purpose: Models the Martian atmosphere very simply
Called by: FGFDMExec
- ------------- Copyright (C) 2004 Jon S. Berndt (jsb@hal-pc.org) -------------
+ ------------- Copyright (C) 2004 Jon S. Berndt (jon@jsbsim.org) -------------
This program is free software; you can redistribute it and/or modify it under
- the terms of the GNU General Public License as published by the Free Software
+ the terms of the GNU Lesser General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
details.
- You should have received a copy of the GNU General Public License along with
+ You should have received a copy of the GNU Lesser General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA.
- Further information about the GNU General Public License can also be found on
+ Further information about the GNU Lesser General Public License can also be found on
the world wide web at http://www.gnu.org.
FUNCTIONAL DESCRIPTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
#include "FGMars.h"
-#include "FGState.h"
+#include <iostream>
+
+using namespace std;
namespace JSBSim {
-static const char *IdSrc = "$Id$";
+static const char *IdSrc = "$Id: FGMars.cpp,v 1.10 2010/02/25 05:21:36 jberndt Exp $";
static const char *IdHdr = ID_MARS;
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Name = "FGMars";
Reng = 53.5 * 44.01;
-/*
- lastIndex = 0;
- h = 0.0;
- psiw = 0.0;
-
- MagnitudedAccelDt = MagnitudeAccel = Magnitude = 0.0;
-// turbType = ttNone;
- turbType = ttStandard;
-// turbType = ttBerndt;
- TurbGain = 0.0;
- TurbRate = 1.0;
-*/
-
bind();
Debug(0);
}
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-/*
-FGMars::~FGMars()
-{
- unbind();
- Debug(1);
-}
-*/
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-bool FGMars::InitModel(void)
-{
- FGModel::InitModel();
-
- Calculate(h);
- SLtemperature = intTemperature;
- SLpressure = intPressure;
- SLdensity = intDensity;
- SLsoundspeed = sqrt(SHRatio*Reng*intTemperature);
- rSLtemperature = 1.0/intTemperature;
- rSLpressure = 1.0/intPressure;
- rSLdensity = 1.0/intDensity;
- rSLsoundspeed = 1.0/SLsoundspeed;
- temperature = &intTemperature;
- pressure = &intPressure;
- density = &intDensity;
-
- useExternal=false;
-
- return true;
-}
-
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-bool FGMars::Run(void)
-{
- if (FGModel::Run()) return true;
- if (FDMExec->Holding()) return false;
-
- //do temp, pressure, and density first
- if (!useExternal) {
- h = Propagate->Geth();
- Calculate(h);
- }
-
- if (turbType != ttNone) {
- Turbulence();
- vWindNED += vTurbulence;
- }
-
- if (vWindNED(1) != 0.0) psiw = atan2( vWindNED(2), vWindNED(1) );
-
- if (psiw < 0) psiw += 2*M_PI;
-
- soundspeed = sqrt(SHRatio*Reng*(*temperature));
-
- Debug(2);
-
- return false;
-}
-
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
void FGMars::Calculate(double altitude)
//cout << "Atmosphere: h=" << altitude << " rho= " << intDensity << endl;
}
-//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-// square a value, but preserve the original sign
-
-static inline double
-square_signed (double value)
-{
- if (value < 0)
- return value * value * -1;
- else
- return value * value;
-}
-
-void FGMars::Turbulence(void)
-{
- switch (turbType) {
- case ttStandard: {
- vDirectiondAccelDt(eX) = 1 - 2.0*(double(rand())/double(RAND_MAX));
- vDirectiondAccelDt(eY) = 1 - 2.0*(double(rand())/double(RAND_MAX));
- vDirectiondAccelDt(eZ) = 1 - 2.0*(double(rand())/double(RAND_MAX));
-
- MagnitudedAccelDt = 1 - 2.0*(double(rand())/double(RAND_MAX)) - Magnitude;
- // Scale the magnitude so that it moves
- // away from the peaks
- MagnitudedAccelDt = ((MagnitudedAccelDt - Magnitude) /
- (1 + fabs(Magnitude)));
- MagnitudeAccel += MagnitudedAccelDt*rate*TurbRate*State->Getdt();
- Magnitude += MagnitudeAccel*rate*State->Getdt();
-
- vDirectiondAccelDt.Normalize();
-
- // deemphasise non-vertical forces
- vDirectiondAccelDt(eX) = square_signed(vDirectiondAccelDt(eX));
- vDirectiondAccelDt(eY) = square_signed(vDirectiondAccelDt(eY));
-
- vDirectionAccel += vDirectiondAccelDt*rate*TurbRate*State->Getdt();
- vDirectionAccel.Normalize();
- vDirection += vDirectionAccel*rate*State->Getdt();
-
- vDirection.Normalize();
-
- // Diminish turbulence within three wingspans
- // of the ground
- vTurbulence = TurbGain * Magnitude * vDirection;
- double HOverBMAC = Auxiliary->GetHOverBMAC();
- if (HOverBMAC < 3.0)
- vTurbulence *= (HOverBMAC / 3.0) * (HOverBMAC / 3.0);
-
- vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
-
- vBodyTurbGrad = Propagate->GetTl2b()*vTurbulenceGrad;
- vTurbPQR(eP) = vBodyTurbGrad(eY)/Aircraft->GetWingSpan();
-// if (Aircraft->GetHTailArm() != 0.0)
-// vTurbPQR(eQ) = vBodyTurbGrad(eZ)/Aircraft->GetHTailArm();
-// else
-// vTurbPQR(eQ) = vBodyTurbGrad(eZ)/10.0;
-
- if (Aircraft->GetVTailArm())
- vTurbPQR(eR) = vBodyTurbGrad(eX)/Aircraft->GetVTailArm();
- else
- vTurbPQR(eR) = vBodyTurbGrad(eX)/10.0;
-
- // Clear the horizontal forces
- // actually felt by the plane, now
- // that we've used them to calculate
- // moments.
- vTurbulence(eX) = 0.0;
- vTurbulence(eY) = 0.0;
-
- break;
- }
- case ttBerndt: {
- vDirectiondAccelDt(eX) = 1 - 2.0*(double(rand())/double(RAND_MAX));
- vDirectiondAccelDt(eY) = 1 - 2.0*(double(rand())/double(RAND_MAX));
- vDirectiondAccelDt(eZ) = 1 - 2.0*(double(rand())/double(RAND_MAX));
-
-
- MagnitudedAccelDt = 1 - 2.0*(double(rand())/double(RAND_MAX)) - Magnitude;
- MagnitudeAccel += MagnitudedAccelDt*rate*State->Getdt();
- Magnitude += MagnitudeAccel*rate*State->Getdt();
-
- vDirectiondAccelDt.Normalize();
- vDirectionAccel += vDirectiondAccelDt*rate*State->Getdt();
- vDirectionAccel.Normalize();
- vDirection += vDirectionAccel*rate*State->Getdt();
-
- // Diminish z-vector within two wingspans
- // of the ground
- double HOverBMAC = Auxiliary->GetHOverBMAC();
- if (HOverBMAC < 2.0)
- vDirection(eZ) *= HOverBMAC / 2.0;
-
- vDirection.Normalize();
-
- vTurbulence = TurbGain*Magnitude * vDirection;
- vTurbulenceGrad = TurbGain*MagnitudeAccel * vDirection;
-
- vBodyTurbGrad = Propagate->GetTl2b()*vTurbulenceGrad;
- vTurbPQR(eP) = vBodyTurbGrad(eY)/Aircraft->GetWingSpan();
- if (Aircraft->GetHTailArm() != 0.0)
- vTurbPQR(eQ) = vBodyTurbGrad(eZ)/Aircraft->GetHTailArm();
- else
- vTurbPQR(eQ) = vBodyTurbGrad(eZ)/10.0;
-
- if (Aircraft->GetVTailArm())
- vTurbPQR(eR) = vBodyTurbGrad(eX)/Aircraft->GetVTailArm();
- else
- vTurbPQR(eR) = vBodyTurbGrad(eX)/10.0;
-
- break;
- }
- default:
- break;
- }
-}
-
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// The bitmasked value choices are as follows:
// unset: In this case (the default) JSBSim would only print
if (debug_lvl & 16) { // Sanity checking
}
if (debug_lvl & 32) { // Turbulence
- if (frame == 0 && from == 2) {
- cout << "vTurbulence(X), vTurbulence(Y), vTurbulence(Z), "
+ if (first_pass && from == 2) {
+ cout << "vTurbulenceNED(X), vTurbulenceNED(Y), vTurbulenceNED(Z), "
<< "vTurbulenceGrad(X), vTurbulenceGrad(Y), vTurbulenceGrad(Z), "
<< "vDirection(X), vDirection(Y), vDirection(Z), "
<< "Magnitude, "
<< "vTurbPQR(P), vTurbPQR(Q), vTurbPQR(R), " << endl;
} else if (from == 2) {
- cout << vTurbulence << ", " << vTurbulenceGrad << ", " << vDirection << ", " << Magnitude << ", " << vTurbPQR << endl;
+ cout << vTurbulenceNED << ", " << vTurbulenceGrad << ", " << vDirection << ", " << Magnitude << ", " << vTurbPQR << endl;
}
}
if (debug_lvl & 64) {