src/FDM/JSBSim/models/atmosphere/FGMSIS.h

   1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2
   3  Header:       FGMSIS.h
   4  Description:  MSIS-00 Atmosphere
   5  Author:       David Culp
   6  Date started: 12/14/03
   7
   8  ------------- Copyright (C) 2003  David P. Culp (davidculp2@comcast.net) ------
   9
  10  This program is free software; you can redistribute it and/or modify it under
  11  the terms of the GNU Lesser General Public License as published by the Free Software
  12  Foundation; either version 2 of the License, or (at your option) any later
  13  version.
  14
  15  This program is distributed in the hope that it will be useful, but WITHOUT
  16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  17  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  18  details.
  19
  20  You should have received a copy of the GNU Lesser General Public License along with
  21  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
  22  Place - Suite 330, Boston, MA  02111-1307, USA.
  23
  24  Further information about the GNU Lesser General Public License can also be found on
  25  the world wide web at http://www.gnu.org.
  26
  27 HISTORY
  28 --------------------------------------------------------------------------------
  29 12/14/03   DPC   Created
  30 01/11/04   DPC   Derive from FGAtmosphere
  31
  32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  33 SENTRY
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  35
  36 #ifndef FGMSIS_H
  37 #define FGMSIS_H
  38
  39 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  40 INCLUDES
  41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  42
  43 #include "models/FGAtmosphere.h"
  44 #include "FGFDMExec.h"
  45
  46 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  47 DEFINITIONS
  48 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  49
  50 #define ID_MSIS "$Id: FGMSIS.h,v 1.9 2011/05/20 03:18:36 jberndt Exp $"
  51
  52 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  53 FORWARD DECLARATIONS
  54 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  55
  56 namespace JSBSim {
  57
  58 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  59 CLASS DOCUMENTATION
  60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  61
  62 /** Models the MSIS-00 atmosphere.
  63     This is a wrapper for the NRL-MSIS-00 model 2001:
  64
  65     This C++ format model wraps the NRLMSISE-00 C source code package - release
  66     20020503
  67
  68     The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and
  69     Doug Drob. They also wrote a NRLMSISE-00 distribution package in
  70     FORTRAN which is available at
  71     http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm
  72
  73     Dominik Brodowski implemented and maintains this C version. You can
  74     reach him at devel@brodo.de. See the file "DOCUMENTATION" for details,
  75     and check http://www.brodo.de/english/pub/nrlmsise/index.html for
  76     updated releases of this package.
  77     @author David Culp
  78     @version $Id: FGMSIS.h,v 1.9 2011/05/20 03:18:36 jberndt Exp $
  79 */
  80
  81 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  82 STRUCT DECLARATIONS
  83 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
  84
  85 struct nrlmsise_flags {
  86   int switches[24];
  87   double sw[24];
  88   double swc[24];
  89 };
  90 /*
  91  *   Switches: to turn on and off particular variations use these switches.
  92  *   0 is off, 1 is on, and 2 is main effects off but cross terms on.
  93  *
  94  *   Standard values are 0 for switch 0 and 1 for switches 1 to 23. The
  95  *   array "switches" needs to be set accordingly by the calling program.
  96  *   The arrays sw and swc are set internally.
  97  *
  98  *   switches[i]:
  99  *    i - explanation
 100  *   -----------------
 101  *    0 - output in centimeters instead of meters
 102  *    1 - F10.7 effect on mean
 103  *    2 - time independent
 104  *    3 - symmetrical annual
 105  *    4 - symmetrical semiannual
 106  *    5 - asymmetrical annual
 107  *    6 - asymmetrical semiannual
 108  *    7 - diurnal
 109  *    8 - semidiurnal
 110  *    9 - daily ap [when this is set to -1 (!) the pointer
 111  *                  ap_a in struct nrlmsise_input must
 112  *                  point to a struct ap_array]
 113  *   10 - all UT/long effects
 114  *   11 - longitudinal
 115  *   12 - UT and mixed UT/long
 116  *   13 - mixed AP/UT/LONG
 117  *   14 - terdiurnal
 118  *   15 - departures from diffusive equilibrium
 119  *   16 - all TINF var
 120  *   17 - all TLB var
 121  *   18 - all TN1 var
 122  *   19 - all S var
 123  *   20 - all TN2 var
 124  *   21 - all NLB var
 125  *   22 - all TN3 var
 126  *   23 - turbo scale height var
 127  */
 128
 129 struct ap_array {
 130   double a[7];
 131 };
 132 /* Array containing the following magnetic values:
 133  *   0 : daily AP
 134  *   1 : 3 hr AP index for current time
 135  *   2 : 3 hr AP index for 3 hrs before current time
 136  *   3 : 3 hr AP index for 6 hrs before current time
 137  *   4 : 3 hr AP index for 9 hrs before current time
 138  *   5 : Average of eight 3 hr AP indicies from 12 to 33 hrs
 139  *           prior to current time
 140  *   6 : Average of eight 3 hr AP indicies from 36 to 57 hrs
 141  *           prior to current time
 142  */
 143
 144
 145 struct nrlmsise_input {
 146   int year;      /* year, currently ignored                           */
 147   int doy;       /* day of year                                       */
 148   double sec;    /* seconds in day (UT)                               */
 149   double alt;    /* altitude in kilometers                            */
 150   double g_lat;  /* geodetic latitude                                 */
 151   double g_long; /* geodetic longitude                                */
 152   double lst;    /* local apparent solar time (hours), see note below */
 153   double f107A;  /* 81 day average of F10.7 flux (centered on DOY)    */
 154   double f107;   /* daily F10.7 flux for previous day                 */
 155   double ap;     /* magnetic index(daily)                             */
 156   struct ap_array *ap_a; /* see above */
 157 };
 158 /*
 159  *   NOTES ON INPUT VARIABLES:
 160  *      UT, Local Time, and Longitude are used independently in the
 161  *      model and are not of equal importance for every situation.
 162  *      For the most physically realistic calculation these three
 163  *      variables should be consistent (lst=sec/3600 + g_long/15).
 164  *      The Equation of Time departures from the above formula
 165  *      for apparent local time can be included if available but
 166  *      are of minor importance.
 167  *
 168  *      f107 and f107A values used to generate the model correspond
 169  *      to the 10.7 cm radio flux at the actual distance of the Earth
 170  *      from the Sun rather than the radio flux at 1 AU. The following
 171  *      site provides both classes of values:
 172  *      ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/
 173  *
 174  *      f107, f107A, and ap effects are neither large nor well
 175  *      established below 80 km and these parameters should be set to
 176  *      150., 150., and 4. respectively.
 177  */
 178
 179
 180
 181 /* ------------------------------------------------------------------- */
 182 /* ------------------------------ OUTPUT ----------------------------- */
 183 /* ------------------------------------------------------------------- */
 184
 185 struct nrlmsise_output {
 186   double d[9];   /* densities    */
 187   double t[2];   /* temperatures */
 188 };
 189 /*
 190  *   OUTPUT VARIABLES:
 191  *      d[0] - HE NUMBER DENSITY(CM-3)
 192  *      d[1] - O NUMBER DENSITY(CM-3)
 193  *      d[2] - N2 NUMBER DENSITY(CM-3)
 194  *      d[3] - O2 NUMBER DENSITY(CM-3)
 195  *      d[4] - AR NUMBER DENSITY(CM-3)
 196  *      d[5] - TOTAL MASS DENSITY(GM/CM3) [includes d[8] in td7d]
 197  *      d[6] - H NUMBER DENSITY(CM-3)
 198  *      d[7] - N NUMBER DENSITY(CM-3)
 199  *      d[8] - Anomalous oxygen NUMBER DENSITY(CM-3)
 200  *      t[0] - EXOSPHERIC TEMPERATURE
 201  *      t[1] - TEMPERATURE AT ALT
 202  *
 203  *
 204  *      O, H, and N are set to zero below 72.5 km
 205  *
 206  *      t[0], Exospheric temperature, is set to global average for
 207  *      altitudes below 120 km. The 120 km gradient is left at global
 208  *      average value for altitudes below 72 km.
 209  *
 210  *      d[5], TOTAL MASS DENSITY, is NOT the same for subroutines GTD7
 211  *      and GTD7D
 212  *
 213  *        SUBROUTINE GTD7 -- d[5] is the sum of the mass densities of the
 214  *        species labeled by indices 0-4 and 6-7 in output variable d.
 215  *        This includes He, O, N2, O2, Ar, H, and N but does NOT include
 216  *        anomalous oxygen (species index 8).
 217  *
 218  *        SUBROUTINE GTD7D -- d[5] is the "effective total mass density
 219  *        for drag" and is the sum of the mass densities of all species
 220  *        in this model, INCLUDING anomalous oxygen.
 221  */
 222
 223
 224
 225 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 226 CLASS DECLARATION
 227 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
 228
 229 class MSIS : public FGAtmosphere
 230 {
 231 public:
 232
 233   /// Constructor
 234   MSIS(FGFDMExec*);
 235   /// Destructor
 236   ~MSIS();
 237   /** Runs the MSIS-00 atmosphere model; called by the Executive
 238       Can pass in a value indicating if the executive is directing the simulation to Hold.
 239       @param Holding if true, the executive has been directed to hold the sim from
 240                      advancing time. Some models may ignore this flag, such as the Input
 241                      model, which may need to be active to listen on a socket for the
 242                      "Resume" command to be given.
 243       @return false if no error */
 244   bool Run(bool Holding);
 245
 246   bool InitModel(void);
 247
 248   /// Does nothing. External control is not allowed.
 249   void UseExternal(void);
 250
 251 private:
 252
 253   void Calculate(int day,      // day of year (1 to 366)
 254                  double sec,   // seconds in day (0.0 to 86400.0)
 255                  double alt,   // altitude, feet
 256                  double lat,   // geodetic latitude, degrees
 257                  double lon    // geodetic longitude, degrees
 258                 );
 259
 260   void Debug(int from);
 261
 262   nrlmsise_flags flags;
 263   nrlmsise_input input;
 264   nrlmsise_output output;
 265   ap_array aph;
 266
 267   /* PARMB */
 268   double gsurf;
 269   double re;
 270
 271   /* GTS3C */
 272   double dd;
 273
 274   /* DMIX */
 275   double dm04, dm16, dm28, dm32, dm40, dm01, dm14;
 276
 277   /* MESO7 */
 278   double meso_tn1[5];
 279   double meso_tn2[4];
 280   double meso_tn3[5];
 281   double meso_tgn1[2];
 282   double meso_tgn2[2];
 283   double meso_tgn3[2];
 284
 285   /* LPOLY */
 286   double dfa;
 287   double plg[4][9];
 288   double ctloc, stloc;
 289   double c2tloc, s2tloc;
 290   double s3tloc, c3tloc;
 291   double apdf, apt[4];
 292
 293   void tselec(struct nrlmsise_flags *flags);
 294   void glatf(double lat, double *gv, double *reff);
 295   double ccor(double alt, double r, double h1, double zh);
 296   double ccor2(double alt, double r, double h1, double zh, double h2);
 297   double scalh(double alt, double xm, double temp);
 298   double dnet(double dd, double dm, double zhm, double xmm, double xm);
 299   void splini(double *xa, double *ya, double *y2a, int n, double x, double *y);
 300   void splint(double *xa, double *ya, double *y2a, int n, double x, double *y);
 301   void spline(double *x, double *y, int n, double yp1, double ypn, double *y2);
 302   double zeta(double zz, double zl);
 303   double densm(double alt, double d0, double xm, double *tz, int mn3, double *zn3,
 304                double *tn3, double *tgn3, int mn2, double *zn2, double *tn2,
 305                double *tgn2);
 306   double densu(double alt, double dlb, double tinf, double tlb, double xm,
 307                double alpha, double *tz, double zlb, double s2, int mn1,
 308                double *zn1, double *tn1, double *tgn1);
 309   double g0(double a, double *p);
 310   double sumex(double ex);
 311   double sg0(double ex, double *p, double *ap);
 312   double globe7(double *p, nrlmsise_input *input, nrlmsise_flags *flags);
 313   double glob7s(double *p, nrlmsise_input *input, nrlmsise_flags *flags);
 314
 315 // GTD7
 316 // Neutral Atmosphere Empirical Model from the surface to lower exosphere.
 317   void gtd7(nrlmsise_input *input, nrlmsise_flags *flags, nrlmsise_output *output);
 318
 319 // GTD7D
 320 // This subroutine provides Effective Total Mass Density for output
 321 // d[5] which includes contributions from "anomalous oxygen" which can
 322 // affect satellite drag above 500 km. See the section "output" for
 323 // additional details.
 324   void gtd7d(nrlmsise_input *input, nrlmsise_flags *flags, nrlmsise_output *output);
 325
 326 // GHP7
 327 // To specify outputs at a pressure level (press) rather than at
 328 // an altitude.
 329   void ghp7(nrlmsise_input *input, nrlmsise_flags *flags, nrlmsise_output *output, double press);
 330
 331 // GTS7
 332 // Thermospheric portion of NRLMSISE-00
 333   void gts7(nrlmsise_input *input, nrlmsise_flags *flags, nrlmsise_output *output);
 334
 335 };
 336
 337 } // namespace JSBSim
 338
 339 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 340 #endif
 341