src/Cockpit/steam.cxx

   1 // steam.cxx - Steam Gauge Calculations
   2 //
   3 // Copyright (C) 2000  Alexander Perry - alex.perry@ieee.org
   4 //
   5 // This program is free software; you can redistribute it and/or
   6 // modify it under the terms of the GNU General Public License as
   7 // published by the Free Software Foundation; either version 2 of the
   8 // License, or (at your option) any later version.
   9 //
  10 // This program is distributed in the hope that it will be useful, but
  11 // WITHOUT ANY WARRANTY; without even the implied warranty of
  12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 // General Public License for more details.
  14 //
  15 // You should have received a copy of the GNU General Public License
  16 // along with this program; if not, write to the Free Software
  17 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18 //
  19 // $Id$
  20
  21
  22 #ifdef HAVE_CONFIG_H
  23 #  include <config.h>
  24 #endif
  25
  26 #include <simgear/compiler.h>
  27
  28 #include STL_IOSTREAM
  29
  30 #include <simgear/constants.h>
  31 #include <simgear/math/sg_types.hxx>
  32 #include <simgear/misc/props.hxx>
  33 #include <Aircraft/aircraft.hxx>
  34 #ifdef FG_WEATHERCM
  35 # include <WeatherCM/FGLocalWeatherDatabase.h>
  36 #else
  37 # include <Environment/environment_mgr.hxx>
  38 # include <Environment/environment.hxx>
  39 #endif
  40
  41 SG_USING_NAMESPACE(std);
  42
  43 #include "radiostack.hxx"
  44 #include "steam.hxx"
  45
  46 static bool isTied = false;
  47
  48
  49 \f
  50 ////////////////////////////////////////////////////////////////////////
  51 // Declare the functions that read the variables
  52 ////////////////////////////////////////////////////////////////////////
  53
  54 double FGSteam::the_STATIC_inhg = 29.92;
  55 double FGSteam::the_ALT_ft = 0.0;  // Indicated altitude
  56 double FGSteam::get_ALT_ft() { _CatchUp(); return the_ALT_ft; }
  57
  58 double FGSteam::the_ALT_datum_mb = 1013.0;
  59 double FGSteam::get_ALT_datum_mb() { return the_ALT_datum_mb; }
  60
  61 void FGSteam::set_ALT_datum_mb ( double datum_mb ) {
  62     the_ALT_datum_mb = datum_mb;
  63 }
  64
  65 double FGSteam::get_ASI_kias() { return fgGetDouble("/velocities/airspeed-kt"); }
  66
  67 double FGSteam::the_VSI_case = 29.92;
  68 double FGSteam::the_VSI_fps = 0.0;
  69 double FGSteam::get_VSI_fps() { _CatchUp(); return the_VSI_fps; }
  70
  71 double FGSteam::the_VACUUM_inhg = 0.0;
  72 double FGSteam::get_VACUUM_inhg() { _CatchUp(); return the_VACUUM_inhg; }
  73
  74 double FGSteam::the_MH_err   = 0.0;
  75 double FGSteam::the_MH_deg   = 0.0;
  76 double FGSteam::the_MH_degps = 0.0;
  77 double FGSteam::get_MH_deg () { _CatchUp(); return the_MH_deg; }
  78
  79 double FGSteam::the_DG_err   = 0.0;
  80 double FGSteam::the_DG_deg   = 0.0;
  81 double FGSteam::the_DG_degps = 0.0;
  82 double FGSteam::the_DG_inhg  = 0.0;
  83 double FGSteam::get_DG_deg () { _CatchUp(); return the_DG_deg; }
  84 double FGSteam::get_DG_err () { _CatchUp(); return the_DG_err; }
  85
  86 void FGSteam::set_DG_err ( double approx_magvar ) {
  87     the_DG_err = approx_magvar;
  88 }
  89
  90 double FGSteam::the_TC_rad   = 0.0;
  91 double FGSteam::the_TC_std   = 0.0;
  92 double FGSteam::get_TC_rad () { _CatchUp(); return the_TC_rad; }
  93 double FGSteam::get_TC_std () { _CatchUp(); return the_TC_std; }
  94
  95 \f
  96 ////////////////////////////////////////////////////////////////////////
  97 // Recording the current time
  98 ////////////////////////////////////////////////////////////////////////
  99
 100
 101 double FGSteam::_UpdateTimePending = 1000000;  /* Forces filters to reset */
 102
 103
 104                                 // FIXME: no need to use static
 105                                 // functions any longer.
 106
 107 void FGSteam::update (double dt)
 108 {
 109     if (!isTied) {
 110         isTied = true;
 111         fgTie("/steam/airspeed-kt", FGSteam::get_ASI_kias);
 112         fgSetArchivable("/steam/airspeed-kt");
 113         fgTie("/steam/altitude-ft", FGSteam::get_ALT_ft);
 114         fgSetArchivable("/steam/altitude-ft");
 115         fgTie("/steam/altimeter-datum-mb",
 116               FGSteam::get_ALT_datum_mb, FGSteam::set_ALT_datum_mb,
 117               false);  /* don't modify the value */
 118         fgSetArchivable("/steam/altimeter-datum-mb");
 119         fgTie("/steam/turn-rate", FGSteam::get_TC_std);
 120         fgSetArchivable("/steam/turn-rate");
 121         fgTie("/steam/slip-skid", FGSteam::get_TC_rad);
 122         fgSetArchivable("/steam/slip-skid");
 123         fgTie("/steam/vertical-speed-fps", FGSteam::get_VSI_fps);
 124         fgSetArchivable("/steam/vertical-speed-fps");
 125         fgTie("/steam/gyro-compass-deg", FGSteam::get_DG_deg);
 126         fgSetArchivable("/steam/gyro-compass-deg");
 127         // fgTie("/steam/adf-deg", FGSteam::get_HackADF_deg);
 128         // fgSetArchivable("/steam/adf-deg");
 129         fgTie("/steam/gyro-compass-error-deg",
 130               FGSteam::get_DG_err, FGSteam::set_DG_err,
 131               false);  /* don't modify the value */
 132         fgSetArchivable("/steam/gyro-compass-error-deg");
 133         fgTie("/steam/mag-compass-deg", FGSteam::get_MH_deg);
 134         fgSetArchivable("/steam/mag-compass-deg");
 135     }
 136     _UpdateTimePending += dt;
 137 }
 138
 139
 140 /* tc should be (elapsed_time_between_updates / desired_smoothing_time) */
 141 void FGSteam::set_lowpass ( double *outthe, double inthe, double tc )
 142 {
 143         if ( tc < 0.0 )
 144         {       if ( tc < -1.0 )
 145                 {       /* time went backwards; kill the filter */
 146                         (*outthe) = inthe;
 147                 } else
 148                 {       /* ignore mildly negative time */
 149                 }
 150         } else
 151         if ( tc < 0.2 )
 152         {       /* Normal mode of operation; fast approximation to exp(-tc) */
 153                 (*outthe) = (*outthe) * ( 1.0 - tc )
 154                           +    inthe  * tc;
 155         } else
 156         if ( tc > 5.0 )
 157         {       /* Huge time step; assume filter has settled */
 158                 (*outthe) = inthe;
 159         } else
 160         {       /* Moderate time step; non linear response */
 161                 double keep = exp ( -tc );
 162                 // printf ( "ARP: Keep is %f\n", keep );
 163                 (*outthe) = (*outthe) * keep
 164                           +    inthe  * ( 1.0 - keep );
 165         }
 166 }
 167
 168
 169 #define INHG_TO_MB 33.86388  /* Inches_of_mercury * INHG_TO_MB == millibars. */
 170
 171 // Convert air pressure to altitude by ICAO Standard Atmosphere
 172 double pressInHgToAltFt(double p_inhg)
 173 {
 174     // Ref. Aviation Formulary, Ed Williams, www.best.com/~williams/avform.htm
 175     const double P_0 = 29.92126;  // Std. MSL pressure, inHg. (=1013.25 mb)
 176     const double p_Tr = 0.2233609 * P_0;  // Pressure at tropopause, same units.
 177     const double h_Tr = 36089.24;  // Alt of tropopause, ft. (=11.0 km)
 178
 179     if (p_inhg > p_Tr)  // 0.0 to 11.0 km
 180         return (1.0 - pow((p_inhg / P_0), 1.0 / 5.2558797)) / 6.8755856e-6;
 181     return h_Tr + log10(p_inhg / p_Tr) / -4.806346e-5;  // 11.0 to 20.0 km
 182     // We could put more code for higher altitudes here.
 183 }
 184
 185
 186 // Convert altitude to air pressure by ICAO Standard Atmosphere
 187 double altFtToPressInHg(double alt_ft)
 188 {
 189     // Ref. Aviation Formulary, Ed Williams, www.best.com/~williams/avform.htm
 190     const double P_0 = 29.92126;  // Std. MSL pressure, inHg. (=1013.25 mb)
 191     const double p_Tr = 0.2233609 * P_0;  // Pressure at tropopause, same units.
 192     const double h_Tr = 36089.24;  // Alt of tropopause, ft. (=11.0 km)
 193
 194     if (alt_ft < h_Tr)  // 0.0 to 11.0 km
 195         return P_0 * pow(1.0 - 6.8755856e-6 * alt_ft, 5.2558797);
 196     return p_Tr * exp(-4.806346e-5 * (alt_ft - h_Tr));  // 11.0 to 20.0 km
 197     // We could put more code for higher altitudes here.
 198 }
 199
 200
 201 \f
 202 ////////////////////////////////////////////////////////////////////////
 203 // Here the fun really begins
 204 ////////////////////////////////////////////////////////////////////////
 205
 206
 207 void FGSteam::_CatchUp()
 208 {
 209   static const SGPropertyNode *heading_deg_node = fgGetNode("/orientation/heading-deg", true);
 210   static const SGPropertyNode *mag_var_deg_node = fgGetNode("/environment/magnetic-variation-deg", true);
 211   static const SGPropertyNode *mag_dip_deg_node = fgGetNode("/environment/magnetic-dip-deg", true);
 212   static const SGPropertyNode *enginge_0_rpm_node = fgGetNode("/engines/engine[0]/rpm", true);
 213
 214   if ( _UpdateTimePending != 0 )
 215   {
 216         double dt = _UpdateTimePending;
 217         double AccN, AccE, AccU;
 218         /* int i, j; */
 219         double d, the_ENGINE_rpm;
 220
 221         /**************************
 222         Someone has called our update function and
 223         it turns out that we are running somewhat behind.
 224         Here, we recalculate everything for a 'dt' second step.
 225         */
 226
 227         /**************************
 228         The ball responds to the acceleration vector in the body
 229         frame, only the components perpendicular to the longitudinal
 230         axis of the aircraft.  This is only related to actual
 231         side slip for a symmetrical aircraft which is not touching
 232         the ground and not changing its attitude.  Math simplifies
 233         by assuming (for small angles) arctan(x)=x in radians.
 234         Obvious failure mode is the absence of liquid in the
 235         tube, which is there to damp the motion, so that instead
 236         the ball will bounce around, hitting the tube ends.
 237         More subtle flaw is having it not move or a travel limit
 238         occasionally due to some dirt in the tube or on the ball.
 239         */
 240         d = - current_aircraft.fdm_state->get_A_Z_pilot();
 241         if ( d < 1 ) d = 1;
 242         set_lowpass ( & the_TC_rad,
 243                 current_aircraft.fdm_state->get_A_Y_pilot () / d,
 244                 dt );
 245
 246         /**************************
 247         The rate of turn indication is from an electric gyro.
 248         We should have it spin up with the master switch.
 249         It is mounted at a funny angle so that it detects
 250         both rate of bank (i.e. rolling into and out of turns)
 251         and the rate of turn (i.e. how fast heading is changing).
 252         */
 253         set_lowpass ( & the_TC_std,
 254                 current_aircraft.fdm_state->get_Phi_dot ()
 255                         * SGD_RADIANS_TO_DEGREES / 20.0 +
 256                 current_aircraft.fdm_state->get_Psi_dot ()
 257                         * SGD_RADIANS_TO_DEGREES / 3.0  , dt );
 258
 259         /**************************
 260         We want to know the pilot accelerations,
 261         to compute the magnetic compass errors.
 262         */
 263         AccN = current_aircraft.fdm_state->get_V_dot_north();
 264         AccE = current_aircraft.fdm_state->get_V_dot_east();
 265         AccU = current_aircraft.fdm_state->get_V_dot_down()
 266              - 9.81 * SG_METER_TO_FEET;
 267         if ( fabs(the_TC_rad) > 0.2 /* 2.0 */ )
 268         {       /* Massive sideslip jams it; it stops turning */
 269                 the_MH_degps = 0.0;
 270                 the_MH_err   = heading_deg_node->getDoubleValue() - the_MH_deg;
 271         } else
 272         {       double MagDip, MagVar, CosDip;
 273                 double FrcN, FrcE, FrcU, AccTot;
 274                 double EdgN, EdgE, EdgU;
 275                 double TrqN, TrqE, TrqU, Torque;
 276                 /* Find a force vector towards exact magnetic north */
 277                 MagVar = mag_var_deg_node->getDoubleValue()
 278                     / SGD_RADIANS_TO_DEGREES;
 279                 MagDip = mag_var_deg_node->getDoubleValue()
 280                     / SGD_RADIANS_TO_DEGREES;
 281                 CosDip = cos ( MagDip );
 282                 FrcN = CosDip * cos ( MagVar );
 283                 FrcE = CosDip * sin ( MagVar );
 284                 FrcU = sin ( MagDip );
 285                 /* Rotation occurs around acceleration axis,
 286                    but axis magnitude is irrelevant.  So compute it. */
 287                 AccTot = AccN*AccN + AccE*AccE + AccU*AccU;
 288                 if ( AccTot > 1.0 )  AccTot = sqrt ( AccTot );
 289                                 else AccTot = 1.0;
 290                 /* Force applies to north marking on compass card */
 291                 EdgN = cos ( the_MH_err / SGD_RADIANS_TO_DEGREES );
 292                 EdgE = sin ( the_MH_err / SGD_RADIANS_TO_DEGREES );
 293                 EdgU = 0.0;
 294                 /* Apply the force to the edge to get torques */
 295                 TrqN = EdgE * FrcU - EdgU * FrcE;
 296                 TrqE = EdgU * FrcN - EdgN * FrcU;
 297                 TrqU = EdgN * FrcE - EdgE * FrcN;
 298                 /* Select the component parallel to the axis */
 299                 Torque = ( TrqN * AccN +
 300                            TrqE * AccE +
 301                            TrqU * AccU ) * 5.0 / AccTot;
 302                 /* The magnetic compass has angular momentum,
 303                    so we apply a torque to it and wait */
 304                 if ( dt < 1.0 )
 305                 {       the_MH_degps= the_MH_degps * (1.0 - dt) - Torque;
 306                         the_MH_err += dt * the_MH_degps;
 307                 }
 308                 if ( the_MH_err >  180.0 ) the_MH_err -= 360.0; else
 309                 if ( the_MH_err < -180.0 ) the_MH_err += 360.0;
 310                 the_MH_deg  = heading_deg_node->getDoubleValue() - the_MH_err;
 311         }
 312
 313         /**************************
 314         This is not actually correct, but provides a
 315         scaling capability for the vacuum pump later on.
 316         When we have a real engine model, we can ask it.
 317         */
 318         the_ENGINE_rpm = enginge_0_rpm_node->getDoubleValue();
 319
 320         /**************************
 321         First, we need to know what the static line is reporting,
 322         which is a whole simulation area in itself.
 323         We filter the actual value by one second to
 324         account for the line impedance of the plumbing.
 325         */
 326 #ifdef FG_WEATHERCM
 327         sgVec3 plane_pos = { cur_fdm_state->get_Latitude(),
 328                              cur_fdm_state->get_Longitude(),
 329                              cur_fdm_state->get_Altitude() * SG_FEET_TO_METER };
 330         double static_inhg = WeatherDatabase->get(plane_pos).AirPressure *
 331             (0.01 / INHG_TO_MB);
 332 #else
 333         double static_inhg = fgGetDouble("/environment/pressure-inhg");
 334 #endif
 335
 336         set_lowpass ( & the_STATIC_inhg, static_inhg, dt );
 337
 338         /*
 339         NO alternate static source error (student feature),
 340         NO possibility of blockage (instructor feature),
 341         NO slip-induced error, important for C172 for example.
 342         */
 343
 344         /**************************
 345         Altimeter.
 346         ICAO standard atmosphere MSL pressure is 1013.25 mb, and pressure
 347         gradient is about 28 ft per mb at MSL increasing to about 32 at
 348         5000 and 38 at 10000 ft.
 349         Standard altimeters apply the subscale offset to the output altitude,
 350         not to the input pressure; I don't know exactly what pressure gradient
 351         they assume for this.  I choose to make it accurate at low altitudes.
 352         Remember, we are trying to simulate a real altimeter, not an ideal one.
 353         */
 354         set_lowpass ( & the_ALT_ft,
 355             pressInHgToAltFt(the_STATIC_inhg) +
 356             (the_ALT_datum_mb - 1013.25) * 28.0, /* accurate at low alt. */
 357             dt * 10 ); /* smoothing time 0.1 s */
 358
 359         /**************************
 360         The VSI case is a low-pass filter of the static line pressure.
 361         The instrument reports the difference, scaled to approx ft.
 362         NO option for student to break glass when static source fails.
 363         NO capability for a fixed non-zero reading when level.
 364         NO capability to have a scaling error of maybe a factor of two.
 365         */
 366         the_VSI_fps = ( the_VSI_case - the_STATIC_inhg )
 367                     * 10000.0; /* manual scaling factor */
 368         set_lowpass ( & the_VSI_case, the_STATIC_inhg, dt/6.0 );
 369
 370         /**************************
 371         The engine driven vacuum pump is directly attached
 372         to the engine shaft, so each engine rotation pumps
 373         a fixed volume.  The amount of air in that volume
 374         is determined by the vacuum line's internal pressure.
 375         The instruments are essentially leaking air like
 376         a fixed source impedance from atmospheric pressure.
 377         The regulator provides a digital limit setting,
 378         which is open circuit unless the pressure drop is big.
 379         Thus, we can compute the vacuum line pressure directly.
 380         We assume that there is negligible reservoir space.
 381         NO failure of the pump supported (yet)
 382         */
 383         the_VACUUM_inhg = the_STATIC_inhg *
 384                 the_ENGINE_rpm / ( the_ENGINE_rpm + 10000.0 );
 385         if ( the_VACUUM_inhg > 5.0 )
 386              the_VACUUM_inhg = 5.0;
 387
 388 /*
 389 > I was merely going to do the engine rpm driven vacuum pump for both
 390 > the AI and DG, have the gyros spin down down in power off descents,
 391 > have it tumble when you exceed the usual pitch or bank limits,
 392 > put in those insidious turning errors ... for now anyway.
 393 */
 394         if ( _UpdateTimePending > 999999 )
 395             the_DG_err = fgGetDouble("/environment/magnetic-variation-deg");
 396         the_DG_degps = 0.01; /* HACK! */
 397         if (dt<1.0) the_DG_err += dt * the_DG_degps;
 398         the_DG_deg = heading_deg_node->getDoubleValue() - the_DG_err;
 399
 400         /**************************
 401         Finished updates, now clear the timer
 402         */
 403         _UpdateTimePending = 0;
 404   } else {
 405       // cout << "0 Updates pending" << endl;
 406   }
 407 }
 408
 409 \f
 410 ////////////////////////////////////////////////////////////////////////
 411 // Everything below is a transient hack; expect it to disappear
 412 ////////////////////////////////////////////////////////////////////////
 413
 414 double FGSteam::get_HackOBS1_deg () {
 415     return current_radiostack->get_nav1_radial();
 416 }
 417
 418 double FGSteam::get_HackOBS2_deg () {
 419     return current_radiostack->get_nav2_radial();
 420 }
 421
 422
 423 // end of steam.cxx