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