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
  35 SG_USING_NAMESPACE(std);
  36
  37 #include "radiostack.hxx"
  38 #include "steam.hxx"
  39
  40 static bool isTied = false;
  41
  42
  43 \f
  44 ////////////////////////////////////////////////////////////////////////
  45 // Declare the functions that read the variables
  46 ////////////////////////////////////////////////////////////////////////
  47
  48 double FGSteam::the_STATIC_inhg = 29.92;
  49 double FGSteam::the_ALT_ft = 0.0;  // Indicated altitude
  50 double FGSteam::get_ALT_ft() { _CatchUp(); return the_ALT_ft; }
  51
  52 double FGSteam::the_ALT_datum_mb = 1013.0;
  53 double FGSteam::get_ALT_datum_mb() { return the_ALT_datum_mb; }
  54
  55 void FGSteam::set_ALT_datum_mb ( double datum_mb ) {
  56     the_ALT_datum_mb = datum_mb;
  57 }
  58
  59 double FGSteam::get_ASI_kias() { return fgGetDouble("/velocities/airspeed-kt"); }
  60
  61 double FGSteam::the_VSI_case = 29.92;
  62 double FGSteam::the_VSI_fps = 0.0;
  63 double FGSteam::get_VSI_fps() { _CatchUp(); return the_VSI_fps; }
  64
  65 double FGSteam::the_VACUUM_inhg = 0.0;
  66 double FGSteam::get_VACUUM_inhg() { _CatchUp(); return the_VACUUM_inhg; }
  67
  68 double FGSteam::the_MH_err   = 0.0;
  69 double FGSteam::the_MH_deg   = 0.0;
  70 double FGSteam::the_MH_degps = 0.0;
  71 double FGSteam::get_MH_deg () { _CatchUp(); return the_MH_deg; }
  72
  73 double FGSteam::the_DG_err   = 0.0;
  74 double FGSteam::the_DG_deg   = 0.0;
  75 double FGSteam::the_DG_degps = 0.0;
  76 double FGSteam::the_DG_inhg  = 0.0;
  77 double FGSteam::get_DG_deg () { _CatchUp(); return the_DG_deg; }
  78 double FGSteam::get_DG_err () { _CatchUp(); return the_DG_err; }
  79
  80 void FGSteam::set_DG_err ( double approx_magvar ) {
  81     the_DG_err = approx_magvar;
  82 }
  83
  84 double FGSteam::the_TC_rad   = 0.0;
  85 double FGSteam::the_TC_std   = 0.0;
  86 double FGSteam::get_TC_rad () { _CatchUp(); return the_TC_rad; }
  87 double FGSteam::get_TC_std () { _CatchUp(); return the_TC_std; }
  88
  89 \f
  90 ////////////////////////////////////////////////////////////////////////
  91 // Recording the current time
  92 ////////////////////////////////////////////////////////////////////////
  93
  94
  95 int FGSteam::_UpdatesPending = 1000000;  /* Forces filter to reset */
  96
  97
  98                                 // FIXME: no need to use static
  99                                 // functions any longer.
 100
 101 void FGSteam::update ( int timesteps )
 102 {
 103     if (!isTied) {
 104         isTied = true;
 105         fgTie("/steam/airspeed-kt", FGSteam::get_ASI_kias);
 106         fgSetArchivable("/steam/airspeed-kt");
 107         fgTie("/steam/altitude-ft", FGSteam::get_ALT_ft);
 108         fgSetArchivable("/steam/altitude-ft");
 109         fgTie("/steam/altimeter-datum-mb",
 110               FGSteam::get_ALT_datum_mb, FGSteam::set_ALT_datum_mb,
 111               false);  /* don't modify the value */
 112         fgSetArchivable("/steam/altimeter-datum-mb");
 113         fgTie("/steam/turn-rate", FGSteam::get_TC_std);
 114         fgSetArchivable("/steam/turn-rate");
 115         fgTie("/steam/slip-skid", FGSteam::get_TC_rad);
 116         fgSetArchivable("/steam/slip-skid");
 117         fgTie("/steam/vertical-speed-fpm", FGSteam::get_VSI_fps);
 118         fgSetArchivable("/steam/vertical-speed-fpm");
 119         fgTie("/steam/gyro-compass-deg", FGSteam::get_DG_deg);
 120         fgSetArchivable("/steam/gyro-compass-deg");
 121         fgTie("/steam/adf-deg", FGSteam::get_HackADF_deg);
 122         fgSetArchivable("/steam/adf-deg");
 123         fgTie("/steam/gyro-compass-error-deg",
 124               FGSteam::get_DG_err, FGSteam::set_DG_err,
 125               false);  /* don't modify the value */
 126         fgSetArchivable("/steam/gyro-compass-error-deg");
 127         fgTie("/steam/mag-compass-deg", FGSteam::get_MH_deg);
 128         fgSetArchivable("/steam/mag-compass-deg");
 129     }
 130     _UpdatesPending += timesteps;
 131 }
 132
 133
 134 /* tc should be (elapsed_time_between_updates / desired_smoothing_time) */
 135 void FGSteam::set_lowpass ( double *outthe, double inthe, double tc )
 136 {
 137         if ( tc < 0.0 )
 138         {       if ( tc < -1.0 )
 139                 {       /* time went backwards; kill the filter */
 140                         (*outthe) = inthe;
 141                 } else
 142                 {       /* ignore mildly negative time */
 143                 }
 144         } else
 145         if ( tc < 0.2 )
 146         {       /* Normal mode of operation; fast approximation to exp(-tc) */
 147                 (*outthe) = (*outthe) * ( 1.0 - tc )
 148                           +    inthe  * tc;
 149         } else
 150         if ( tc > 5.0 )
 151         {       /* Huge time step; assume filter has settled */
 152                 (*outthe) = inthe;
 153         } else
 154         {       /* Moderate time step; non linear response */
 155                 double keep = exp ( -tc );
 156                 // printf ( "ARP: Keep is %f\n", keep );
 157                 (*outthe) = (*outthe) * keep
 158                           +    inthe  * ( 1.0 - keep );
 159         }
 160 }
 161
 162
 163 #define INHG_TO_MB 33.86388  /* Inches_of_mercury * INHG_TO_MB == millibars. */
 164
 165 // Convert air pressure to altitude by ICAO Standard Atmosphere
 166 double pressInHgToAltFt(double p_inhg)
 167 {
 168     // Ref. Aviation Formulary, Ed Williams, www.best.com/~williams/avform.htm
 169     const double P_0 = 29.92126;  // Std. MSL pressure, inHg. (=1013.25 mb)
 170     const double p_Tr = 0.2233609 * P_0;  // Pressure at tropopause, same units.
 171     const double h_Tr = 36089.24;  // Alt of tropopause, ft. (=11.0 km)
 172
 173     if (p_inhg > p_Tr)  // 0.0 to 11.0 km
 174         return (1.0 - pow((p_inhg / P_0), 1.0 / 5.2558797)) / 6.8755856e-6;
 175     return h_Tr + log10(p_inhg / p_Tr) / -4.806346e-5;  // 11.0 to 20.0 km
 176     // We could put more code for higher altitudes here.
 177 }
 178
 179
 180 // Convert altitude to air pressure by ICAO Standard Atmosphere
 181 double altFtToPressInHg(double alt_ft)
 182 {
 183     // Ref. Aviation Formulary, Ed Williams, www.best.com/~williams/avform.htm
 184     const double P_0 = 29.92126;  // Std. MSL pressure, inHg. (=1013.25 mb)
 185     const double p_Tr = 0.2233609 * P_0;  // Pressure at tropopause, same units.
 186     const double h_Tr = 36089.24;  // Alt of tropopause, ft. (=11.0 km)
 187
 188     if (alt_ft < h_Tr)  // 0.0 to 11.0 km
 189         return P_0 * pow(1.0 - 6.8755856e-6 * alt_ft, 5.2558797);
 190     return p_Tr * exp(-4.806346e-5 * (alt_ft - h_Tr));  // 11.0 to 20.0 km
 191     // We could put more code for higher altitudes here.
 192 }
 193
 194
 195 \f
 196 ////////////////////////////////////////////////////////////////////////
 197 // Here the fun really begins
 198 ////////////////////////////////////////////////////////////////////////
 199
 200
 201 void FGSteam::_CatchUp()
 202 { if ( _UpdatesPending != 0 )
 203   {     double dt = _UpdatesPending * 1.0 /
 204             fgGetInt("/sim/model-hz"); // FIXME: inefficient
 205         double AccN, AccE, AccU;
 206         /* int i, j; */
 207         double d, the_ENGINE_rpm;
 208
 209 #if 0
 210         /**************************
 211         There is the possibility that this is the first call.
 212         If this is the case, we will emit the feature registrations
 213         just to be on the safe side.  Doing it more than once will
 214         waste CPU time but doesn't hurt anything really.
 215         */
 216         if ( _UpdatesPending > 999999 )
 217         { FGFeature::register_int (    "Avionics/NAV1/Localizer", &NAV1_LOC );
 218           FGFeature::register_double ( "Avionics/NAV1/Latitude",  &NAV1_Lat );
 219           FGFeature::register_double ( "Avionics/NAV1/Longitude", &NAV1_Lon );
 220           FGFeature::register_double ( "Avionics/NAV1/Radial",    &NAV1_Rad );
 221           FGFeature::register_double ( "Avionics/NAV1/Altitude",  &NAV1_Alt );
 222           FGFeature::register_int (    "Avionics/NAV2/Localizer", &NAV2_LOC );
 223           FGFeature::register_double ( "Avionics/NAV2/Latitude",  &NAV2_Lat );
 224           FGFeature::register_double ( "Avionics/NAV2/Longitude", &NAV2_Lon );
 225           FGFeature::register_double ( "Avionics/NAV2/Radial",    &NAV2_Rad );
 226           FGFeature::register_double ( "Avionics/NAV2/Altitude",  &NAV2_Alt );
 227           FGFeature::register_double ( "Avionics/ADF/Latitude",   &ADF_Lat );
 228           FGFeature::register_double ( "Avionics/ADF/Longitude",  &ADF_Lon );
 229         }
 230 #endif
 231
 232         /**************************
 233         Someone has called our update function and
 234         it turns out that we are running somewhat behind.
 235         Here, we recalculate everything for a 'dt' second step.
 236         */
 237
 238         /**************************
 239         The ball responds to the acceleration vector in the body
 240         frame, only the components perpendicular to the longitudinal
 241         axis of the aircraft.  This is only related to actual
 242         side slip for a symmetrical aircraft which is not touching
 243         the ground and not changing its attitude.  Math simplifies
 244         by assuming (for small angles) arctan(x)=x in radians.
 245         Obvious failure mode is the absence of liquid in the
 246         tube, which is there to damp the motion, so that instead
 247         the ball will bounce around, hitting the tube ends.
 248         More subtle flaw is having it not move or a travel limit
 249         occasionally due to some dirt in the tube or on the ball.
 250         */
 251         d = - current_aircraft.fdm_state->get_A_Z_pilot();
 252         if ( d < 1 ) d = 1;
 253         set_lowpass ( & the_TC_rad,
 254                 current_aircraft.fdm_state->get_A_Y_pilot () / d,
 255                 dt );
 256
 257         /**************************
 258         The rate of turn indication is from an electric gyro.
 259         We should have it spin up with the master switch.
 260         It is mounted at a funny angle so that it detects
 261         both rate of bank (i.e. rolling into and out of turns)
 262         and the rate of turn (i.e. how fast heading is changing).
 263         */
 264         set_lowpass ( & the_TC_std,
 265                 current_aircraft.fdm_state->get_Phi_dot ()
 266                         * SGD_RADIANS_TO_DEGREES / 20.0 +
 267                 current_aircraft.fdm_state->get_Psi_dot ()
 268                         * SGD_RADIANS_TO_DEGREES / 3.0  , dt );
 269
 270         /**************************
 271         We want to know the pilot accelerations,
 272         to compute the magnetic compass errors.
 273         */
 274         AccN = current_aircraft.fdm_state->get_V_dot_north();
 275         AccE = current_aircraft.fdm_state->get_V_dot_east();
 276         AccU = current_aircraft.fdm_state->get_V_dot_down()
 277              - 9.81 / 0.3;
 278         if ( fabs(the_TC_rad) > 0.2 )
 279         {       /* Massive sideslip jams it; it stops turning */
 280                 the_MH_degps = 0.0;
 281                 the_MH_err   = fgGetDouble("/orientation/heading-deg") - the_MH_deg;
 282         } else
 283         {       double MagDip, MagVar, CosDip;
 284                 double FrcN, FrcE, FrcU, AccTot;
 285                 double EdgN, EdgE, EdgU;
 286                 double TrqN, TrqE, TrqU, Torque;
 287                 /* Find a force vector towards exact magnetic north */
 288                 MagVar = fgGetDouble("/environment/magnetic-variation-deg")
 289                     / SGD_RADIANS_TO_DEGREES;
 290                 MagDip = fgGetDouble("/environment/magnetic-dip-deg")
 291                     / SGD_RADIANS_TO_DEGREES;
 292                 CosDip = cos ( MagDip );
 293                 FrcN = CosDip * cos ( MagVar );
 294                 FrcE = CosDip * sin ( MagVar );
 295                 FrcU = sin ( MagDip );
 296                 /* Rotation occurs around acceleration axis,
 297                    but axis magnitude is irrelevant.  So compute it. */
 298                 AccTot = AccN*AccN + AccE*AccE + AccU*AccU;
 299                 if ( AccTot > 1.0 )  AccTot = sqrt ( AccTot );
 300                                 else AccTot = 1.0;
 301                 /* Force applies to north marking on compass card */
 302                 EdgN = cos ( the_MH_err / SGD_RADIANS_TO_DEGREES );
 303                 EdgE = sin ( the_MH_err / SGD_RADIANS_TO_DEGREES );
 304                 EdgU = 0.0;
 305                 /* Apply the force to the edge to get torques */
 306                 TrqN = EdgE * FrcU - EdgU * FrcE;
 307                 TrqE = EdgU * FrcN - EdgN * FrcU;
 308                 TrqU = EdgN * FrcE - EdgE * FrcN;
 309                 /* Select the component parallel to the axis */
 310                 Torque = ( TrqN * AccN +
 311                            TrqE * AccE +
 312                            TrqU * AccU ) * 5.0 / AccTot;
 313                 /* The magnetic compass has angular momentum,
 314                    so we apply a torque to it and wait */
 315                 if ( dt < 1.0 )
 316                 {       the_MH_degps= the_MH_degps * (1.0 - dt) - Torque;
 317                         the_MH_err += dt * the_MH_degps;
 318                 }
 319                 if ( the_MH_err >  180.0 ) the_MH_err -= 360.0; else
 320                 if ( the_MH_err < -180.0 ) the_MH_err += 360.0;
 321                 the_MH_deg  = fgGetDouble("/orientation/heading-deg") - the_MH_err;
 322         }
 323
 324         /**************************
 325         This is not actually correct, but provides a
 326         scaling capability for the vacuum pump later on.
 327         When we have a real engine model, we can ask it.
 328         */
 329         the_ENGINE_rpm = globals->get_controls()->get_throttle(0) * 26.0;
 330
 331         /**************************
 332         First, we need to know what the static line is reporting,
 333         which is a whole simulation area in itself.  For now, we cheat.
 334         We filter the actual value by one second to
 335         account for the line impedance of the plumbing.
 336         */
 337         double static_inhg
 338             = altFtToPressInHg(fgGetDouble("/position/altitude-ft"));
 339         set_lowpass ( & the_STATIC_inhg, static_inhg, dt );
 340
 341         /*
 342         NO alternate static source error (student feature),
 343         NO possibility of blockage (instructor feature),
 344         NO slip-induced error, important for C172 for example.
 345         */
 346
 347         /**************************
 348         Altimeter.
 349         ICAO standard atmosphere MSL pressure is 1013.25 mb, and pressure
 350         gradient is about 28 ft per mb at MSL increasing to about 32 at
 351         5000 and 38 at 10000 ft.
 352         Standard altimeters apply the subscale offset to the output altitude,
 353         not to the input pressure; I don't know exactly what pressure gradient
 354         they assume for this.  I choose to make it accurate at low altitudes.
 355         Remember, we are trying to simulate a real altimeter, not an ideal one.
 356         */
 357         set_lowpass ( & the_ALT_ft,
 358             pressInHgToAltFt(the_STATIC_inhg) +
 359             (the_ALT_datum_mb - 1013.25) * 28.0, /* accurate at low alt. */
 360             dt * 10 ); /* smoothing time 0.1 s */
 361
 362         /**************************
 363         The VSI case is a low-pass filter of the static line pressure.
 364         The instrument reports the difference, scaled to approx ft.
 365         NO option for student to break glass when static source fails.
 366         NO capability for a fixed non-zero reading when level.
 367         NO capability to have a scaling error of maybe a factor of two.
 368         */
 369         the_VSI_fps = ( the_VSI_case - the_STATIC_inhg )
 370                     * 10000.0; /* manual scaling factor */
 371         set_lowpass ( & the_VSI_case, the_STATIC_inhg, dt/6.0 );
 372
 373         /**************************
 374         The engine driven vacuum pump is directly attached
 375         to the engine shaft, so each engine rotation pumps
 376         a fixed volume.  The amount of air in that volume
 377         is determined by the vacuum line's internal pressure.
 378         The instruments are essentially leaking air like
 379         a fixed source impedance from atmospheric pressure.
 380         The regulator provides a digital limit setting,
 381         which is open circuit unless the pressure drop is big.
 382         Thus, we can compute the vacuum line pressure directly.
 383         We assume that there is negligible reservoir space.
 384         NO failure of the pump supported (yet)
 385         */
 386         the_VACUUM_inhg = the_STATIC_inhg *
 387                 the_ENGINE_rpm / ( the_ENGINE_rpm + 10000.0 );
 388         if ( the_VACUUM_inhg > 5.0 )
 389              the_VACUUM_inhg = 5.0;
 390
 391 /*
 392 > I was merely going to do the engine rpm driven vacuum pump for both
 393 > the AI and DG, have the gyros spin down down in power off descents,
 394 > have it tumble when you exceed the usual pitch or bank limits,
 395 > put in those insidious turning errors ... for now anyway.
 396 */
 397         if ( _UpdatesPending > 999999 )
 398             the_DG_err = fgGetDouble("/environment/magnetic-variation-deg");
 399         the_DG_degps = 0.01; /* HACK! */
 400         if (dt<1.0) the_DG_err += dt * the_DG_degps;
 401         the_DG_deg = fgGetDouble("/orientation/heading-deg") - the_DG_err;
 402
 403         /**************************
 404         Finished updates, now clear the timer
 405         */
 406         _UpdatesPending = 0;
 407   } else {
 408       // cout << "0 Updates pending" << endl;
 409   }
 410 }
 411
 412 \f
 413 ////////////////////////////////////////////////////////////////////////
 414 // Everything below is a transient hack; expect it to disappear
 415 ////////////////////////////////////////////////////////////////////////
 416
 417
 418 #if 0
 419
 420 double FGSteam::get_HackGS_deg () {
 421     if ( current_radiostack->get_nav1_inrange() &&
 422          current_radiostack->get_nav1_has_gs() )
 423     {
 424         double x = current_radiostack->get_nav1_gs_dist();
 425         double y = (fgGetDouble("/position/altitude-ft")
 426                     - current_radiostack->get_nav1_elev())
 427             * SG_FEET_TO_METER;
 428         double angle = atan2( y, x ) * SGD_RADIANS_TO_DEGREES;
 429         return (current_radiostack->get_nav1_target_gs() - angle) * 5.0;
 430     } else {
 431         return 0.0;
 432     }
 433 }
 434
 435
 436 double FGSteam::get_HackVOR1_deg () {
 437     double r;
 438
 439     if ( current_radiostack->get_nav1_inrange() ) {
 440         r = current_radiostack->get_nav1_heading()
 441             - current_radiostack->get_nav1_radial();
 442         // cout << "Radial = " << current_radiostack->get_nav1_radial()
 443         //      << "  Bearing = " << current_radiostack->get_nav1_heading()
 444         //      << endl;
 445
 446         if (r> 180.0) r-=360.0; else
 447             if (r<-180.0) r+=360.0;
 448         if ( fabs(r) > 90.0 )
 449             r = ( r<0.0 ? -r-180.0 : -r+180.0 );
 450         // According to Robin Peel, the ILS is 4x more sensitive than a vor
 451         if ( current_radiostack->get_nav1_loc() ) r *= 4.0;
 452     } else {
 453         r = 0.0;
 454     }
 455
 456     return r;
 457 }
 458
 459
 460 double FGSteam::get_HackVOR2_deg () {
 461     double r;
 462
 463     if ( current_radiostack->get_nav2_inrange() ) {
 464         r = current_radiostack->get_nav2_heading()
 465             - current_radiostack->get_nav2_radial();
 466         // cout << "Radial = " << current_radiostack->get_nav1_radial()
 467         // << "  Bearing = " << current_radiostack->get_nav1_heading() << endl;
 468
 469         if (r> 180.0) r-=360.0; else
 470             if (r<-180.0) r+=360.0;
 471         if ( fabs(r) > 90.0 )
 472             r = ( r<0.0 ? -r-180.0 : -r+180.0 );
 473     } else {
 474         r = 0.0;
 475     }
 476
 477     return r;
 478 }
 479 #endif
 480
 481
 482 double FGSteam::get_HackOBS1_deg () {
 483     return current_radiostack->get_nav1_radial();
 484 }
 485
 486
 487 double FGSteam::get_HackOBS2_deg () {
 488     return current_radiostack->get_nav2_radial();
 489 }
 490
 491
 492 double FGSteam::get_HackADF_deg () {
 493     static double last_r = 0;
 494
 495     if ( current_radiostack->get_adf_inrange() ) {
 496         double r = current_radiostack->get_adf_heading()
 497             - fgGetDouble("/orientation/heading-deg");
 498         last_r = r;
 499         // cout << "Radial = " << current_radiostack->get_adf_heading()
 500         //      << "  Heading = " << FGBFI::getHeading() << endl;
 501         return r;
 502     } else {
 503         return last_r;
 504     }
 505 }
 506
 507
 508 // end of steam.cxx