]> git.mxchange.org Git - flightgear.git/blob - src/FDM/JSBSim/models/atmosphere/FGMSIS.h
Merge branch 'next' into durk-atc
[flightgear.git] / 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