]> git.mxchange.org Git - flightgear.git/blob - src/FDM/JSBSim/models/atmosphere/FGMSIS.cpp
Merge branch 'next' into durk-atc
[flightgear.git] / src / FDM / JSBSim / models / atmosphere / FGMSIS.cpp
1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2
3  Module:       FGMSIS.cpp
4  Author:       David Culp
5                (incorporated into C++ JSBSim class heirarchy, see model authors below)
6  Date started: 12/14/03
7  Purpose:      Models the MSIS-00 atmosphere
8
9  ------------- Copyright (C) 2003  David P. Culp (davidculp2@comcast.net) ------
10
11  This program is free software; you can redistribute it and/or modify it under
12  the terms of the GNU Lesser General Public License as published by the Free Software
13  Foundation; either version 2 of the License, or (at your option) any later
14  version.
15
16  This program is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18  FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
19  details.
20
21  You should have received a copy of the GNU Lesser General Public License along with
22  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
23  Place - Suite 330, Boston, MA  02111-1307, USA.
24
25  Further information about the GNU Lesser General Public License can also be found on
26  the world wide web at http://www.gnu.org.
27
28 FUNCTIONAL DESCRIPTION
29 --------------------------------------------------------------------------------
30 Models the MSIS-00 atmosphere. Provides temperature and density to FGAtmosphere,
31 given day-of-year, time-of-day, altitude, latitude, longitude and local time.
32
33 HISTORY
34 --------------------------------------------------------------------------------
35 12/14/03   DPC   Created
36 01/11/04   DPC   Derived from FGAtmosphere
37
38  --------------------------------------------------------------------
39  ---------  N R L M S I S E - 0 0    M O D E L    2 0 0 1  ----------
40  --------------------------------------------------------------------
41
42  This file is part of the NRLMSISE-00  C source code package - release
43  20020503
44
45  The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and
46  Doug Drob. They also wrote a NRLMSISE-00 distribution package in
47  FORTRAN which is available at
48  http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm
49
50  Dominik Brodowski implemented and maintains this C version. You can
51  reach him at devel@brodo.de. See the file "DOCUMENTATION" for details,
52  and check http://www.brodo.de/english/pub/nrlmsise/index.html for
53  updated releases of this package.
54 */
55
56 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 INCLUDES
58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
59
60 #include "FGMSIS.h"
61 #include "models/FGAuxiliary.h"
62 #include <cmath>          /* maths functions */
63 #include <iostream>        // for cout, endl
64
65 using namespace std;
66
67 namespace JSBSim {
68
69 static const char *IdSrc = "$Id: FGMSIS.cpp,v 1.17 2011/05/20 03:18:36 jberndt Exp $";
70 static const char *IdHdr = ID_MSIS;
71
72 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73 EXTERNAL GLOBAL DATA
74 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
75
76   /* POWER7 */
77   extern double pt[150];
78   extern double pd[9][150];
79   extern double ps[150];
80   extern double pdl[2][25];
81   extern double ptl[4][100];
82   extern double pma[10][100];
83   extern double sam[100];
84
85   /* LOWER7 */
86   extern double ptm[10];
87   extern double pdm[8][10];
88   extern double pavgm[10];
89
90 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 CLASS IMPLEMENTATION
92 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
93
94
95 MSIS::MSIS(FGFDMExec* fdmex) : FGAtmosphere(fdmex)
96 {
97   Name = "MSIS";
98
99   for (int i=0; i<9; i++) output.d[i] = 0.0;
100   for (int i=0; i<2; i++) output.t[i] = 0.0;
101
102   dm04 = dm16 = dm28 = dm32 = dm40 = dm01 = dm14 = dfa = 0.0;
103
104   for (int i=0; i<5; i++) meso_tn1[i] = 0.0;
105   for (int i=0; i<4; i++) meso_tn2[i] = 0.0;
106   for (int i=0; i<5; i++) meso_tn3[i] = 0.0;
107   for (int i=0; i<2; i++) meso_tgn1[i] = 0.0;
108   for (int i=0; i<2; i++) meso_tgn2[i] = 0.0;
109   for (int i=0; i<2; i++) meso_tgn3[i] = 0.0;
110
111   Debug(0);
112 }
113
114 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115
116 MSIS::~MSIS()
117 {
118   Debug(1);
119 }
120
121 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
122
123 bool MSIS::InitModel(void)
124 {
125   unsigned int i;
126
127   flags.switches[0] = 0;
128   flags.sw[0] = 0;
129   flags.swc[0] = 0;
130   for (i=1;i<24;i++) {
131     flags.switches[i] = 1;
132     flags.sw[i] = 1;
133     flags.swc[i] = 1;
134   }
135
136   for (i=0;i<7;i++) aph.a[i] = 100.0;
137
138   // set some common magnetic flux values
139   input.f107A = 150.0;
140   input.f107 = 150.0;
141   input.ap = 4.0;
142
143   UseInternal();
144
145   SLtemperature = intTemperature = 518.0;
146   SLpressure    = intPressure = 2116.7;
147   SLdensity     = intDensity = 0.002378;
148   SLsoundspeed  = sqrt(2403.0832 * SLtemperature);
149   rSLtemperature = 1.0/intTemperature;
150   rSLpressure    = 1.0/intPressure;
151   rSLdensity     = 1.0/intDensity;
152   rSLsoundspeed  = 1.0/SLsoundspeed;
153
154   return true;
155 }
156
157 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
158
159 bool MSIS::Run(bool Holding)
160 {
161   if (FGModel::Run(Holding)) return true;
162   if (Holding) return false;
163
164   RunPreFunctions();
165
166   h = FDMExec->GetPropagate()->GetAltitudeASL();
167
168   //do temp, pressure, and density first
169   if (!useExternal) {
170     // get sea-level values
171     Calculate(FDMExec->GetAuxiliary()->GetDayOfYear(),
172               FDMExec->GetAuxiliary()->GetSecondsInDay(),
173               0.0,
174               FDMExec->GetPropagate()->GetLocation().GetLatitudeDeg(),
175               FDMExec->GetPropagate()->GetLocation().GetLongitudeDeg());
176     SLtemperature = output.t[1] * 1.8;
177     SLdensity     = output.d[5] * 1.940321;
178     SLpressure    = 1716.488 * SLdensity * SLtemperature;
179     SLsoundspeed  = sqrt(2403.0832 * SLtemperature);
180     rSLtemperature = 1.0/SLtemperature;
181     rSLpressure    = 1.0/SLpressure;
182     rSLdensity     = 1.0/SLdensity;
183     rSLsoundspeed  = 1.0/SLsoundspeed;
184
185     // get at-altitude values
186     Calculate(FDMExec->GetAuxiliary()->GetDayOfYear(),
187               FDMExec->GetAuxiliary()->GetSecondsInDay(),
188               h,
189               FDMExec->GetPropagate()->GetLocation().GetLatitudeDeg(),
190               FDMExec->GetPropagate()->GetLocation().GetLongitudeDeg());
191     intTemperature = output.t[1] * 1.8;
192     intDensity     = output.d[5] * 1.940321;
193     intPressure    = 1716.488 * intDensity * intTemperature;
194     //cout << "T=" << intTemperature << " D=" << intDensity << " P=";
195     //cout << intPressure << " a=" << soundspeed << endl;
196   }
197
198   CalculateDerived();
199
200   RunPostFunctions();
201
202   Debug(2);
203
204   return false;
205 }
206
207 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
208
209 void MSIS::Calculate(int day, double sec, double alt, double lat, double lon)
210 {
211   input.year = 2000;
212   input.doy = day;
213   input.sec = sec;
214   input.alt = alt / 3281;  //feet to kilometers
215   input.g_lat = lat;
216   input.g_long = lon;
217
218   input.lst = (sec/3600) + (lon/15);
219   if (input.lst > 24.0) input.lst -= 24.0;
220   if (input.lst < 0.0) input.lst = 24 - input.lst;
221
222   gtd7d(&input, &flags, &output);
223 }
224
225 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
226
227
228 void MSIS::UseExternal(void){
229   // do nothing, external control not allowed
230 }
231
232
233 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
234
235
236 void MSIS::tselec(struct nrlmsise_flags *flags)
237 {
238   int i;
239   for (i=0;i<24;i++) {
240     if (i!=9) {
241       if (flags->switches[i]==1)
242         flags->sw[i]=1;
243       else
244         flags->sw[i]=0;
245       if (flags->switches[i]>0)
246         flags->swc[i]=1;
247       else
248         flags->swc[i]=0;
249     } else {
250       flags->sw[i]=flags->switches[i];
251       flags->swc[i]=flags->switches[i];
252     }
253   }
254 }
255
256
257 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
258
259 void MSIS::glatf(double lat, double *gv, double *reff)
260 {
261   double dgtr = 1.74533E-2;
262   double c2;
263   c2 = cos(2.0*dgtr*lat);
264   *gv = 980.616 * (1.0 - 0.0026373 * c2);
265   *reff = 2.0 * (*gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5;
266 }
267
268 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
269
270 double MSIS::ccor(double alt, double r, double h1, double zh)
271 {
272 /*        CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS
273  *         ALT - altitude
274  *         R - target ratio
275  *         H1 - transition scale length
276  *         ZH - altitude of 1/2 R
277  */
278   double e;
279   double ex;
280   e = (alt - zh) / h1;
281   if (e>70)
282     return exp(0.0);
283   if (e<-70)
284     return exp(r);
285   ex = exp(e);
286   e = r / (1.0 + ex);
287   return exp(e);
288 }
289
290 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
291
292 double MSIS::ccor2(double alt, double r, double h1, double zh, double h2)
293 {
294 /*        CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS
295  *         ALT - altitude
296  *         R - target ratio
297  *         H1 - transition scale length
298  *         ZH - altitude of 1/2 R
299  *         H2 - transition scale length #2 ?
300  */
301   double e1, e2;
302   double ex1, ex2;
303   double ccor2v;
304   e1 = (alt - zh) / h1;
305   e2 = (alt - zh) / h2;
306   if ((e1 > 70) || (e2 > 70))
307     return exp(0.0);
308   if ((e1 < -70) && (e2 < -70))
309     return exp(r);
310   ex1 = exp(e1);
311   ex2 = exp(e2);
312   ccor2v = r / (1.0 + 0.5 * (ex1 + ex2));
313   return exp(ccor2v);
314 }
315
316 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
317
318 double MSIS::scalh(double alt, double xm, double temp)
319 {
320   double g;
321   double rgas=831.4;
322   g = gsurf / (pow((1.0 + alt/re),2.0));
323   g = rgas * temp / (g * xm);
324   return g;
325 }
326
327 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
328
329 double MSIS::dnet (double dd, double dm, double zhm, double xmm, double xm)
330 {
331 /*       TURBOPAUSE CORRECTION FOR MSIS MODELS
332  *        Root mean density
333  *         DD - diffusive density
334  *         DM - full mixed density
335  *         ZHM - transition scale length
336  *         XMM - full mixed molecular weight
337  *         XM  - species molecular weight
338  *         DNET - combined density
339  */
340   double a;
341   double ylog;
342   a  = zhm / (xmm-xm);
343   if (!((dm>0) && (dd>0))) {
344     cerr << "dnet log error " << dm << ' ' << dd << ' ' << xm << ' ' << endl;
345     if ((dd==0) && (dm==0))
346       dd=1;
347     if (dm==0)
348       return dd;
349     if (dd==0)
350       return dm;
351   }
352   ylog = a * log(dm/dd);
353   if (ylog<-10)
354     return dd;
355   if (ylog>10)
356     return dm;
357   a = dd*pow((1.0 + exp(ylog)),(1.0/a));
358   return a;
359 }
360
361 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
362
363 void MSIS::splini (double *xa, double *ya, double *y2a, int n, double x, double *y)
364 {
365 /*      INTEGRATE CUBIC SPLINE FUNCTION FROM XA(1) TO X
366  *       XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
367  *       Y2A: ARRAY OF SECOND DERIVATIVES
368  *       N: SIZE OF ARRAYS XA,YA,Y2A
369  *       X: ABSCISSA ENDPOINT FOR INTEGRATION
370  *       Y: OUTPUT VALUE
371  */
372   double yi=0;
373   int klo=0;
374   int khi=1;
375   double xx=0.0, h=0.0, a=0.0, b=0.0, a2=0.0, b2=0.0;
376   while ((x>xa[klo]) && (khi<n)) {
377     xx=x;
378     if (khi<(n-1)) {
379       if (x<xa[khi])
380         xx=x;
381       else
382         xx=xa[khi];
383     }
384     h = xa[khi] - xa[klo];
385     a = (xa[khi] - xx)/h;
386     b = (xx - xa[klo])/h;
387     a2 = a*a;
388     b2 = b*b;
389     yi += ((1.0 - a2) * ya[klo] / 2.0 + b2 * ya[khi] / 2.0 + ((-(1.0+a2*a2)/4.0 + a2/2.0) * y2a[klo] + (b2*b2/4.0 - b2/2.0) * y2a[khi]) * h * h / 6.0) * h;
390     klo++;
391     khi++;
392   }
393   *y = yi;
394 }
395
396 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
397
398 void MSIS::splint (double *xa, double *ya, double *y2a, int n, double x, double *y)
399 {
400 /*      CALCULATE CUBIC SPLINE INTERP VALUE
401  *       ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL.
402  *       XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
403  *       Y2A: ARRAY OF SECOND DERIVATIVES
404  *       N: SIZE OF ARRAYS XA,YA,Y2A
405  *       X: ABSCISSA FOR INTERPOLATION
406  *       Y: OUTPUT VALUE
407  */
408   int klo=0;
409   int khi=n-1;
410   int k;
411   double h;
412   double a, b, yi;
413   while ((khi-klo)>1) {
414     k=(khi+klo)/2;
415     if (xa[k]>x)
416       khi=k;
417     else
418       klo=k;
419   }
420   h = xa[khi] - xa[klo];
421   if (h==0.0)
422     cerr << "bad XA input to splint" << endl;
423   a = (xa[khi] - x)/h;
424   b = (x - xa[klo])/h;
425   yi = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0;
426   *y = yi;
427 }
428
429 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
430
431 void MSIS::spline (double *x, double *y, int n, double yp1, double ypn, double *y2)
432 {
433 /*       CALCULATE 2ND DERIVATIVES OF CUBIC SPLINE INTERP FUNCTION
434  *       ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL
435  *       X,Y: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
436  *       N: SIZE OF ARRAYS X,Y
437  *       YP1,YPN: SPECIFIED DERIVATIVES AT X[0] AND X[N-1]; VALUES
438  *                >= 1E30 SIGNAL SIGNAL SECOND DERIVATIVE ZERO
439  *       Y2: OUTPUT ARRAY OF SECOND DERIVATIVES
440  */
441   double *u;
442   double sig, p, qn, un;
443   int i, k;
444   u=new double[n];
445   if (u==NULL) {
446     cerr << "Out Of Memory in spline - ERROR" << endl;
447     return;
448   }
449   if (yp1>0.99E30) {
450     y2[0]=0;
451     u[0]=0;
452   } else {
453     y2[0]=-0.5;
454     u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
455   }
456   for (i=1;i<(n-1);i++) {
457     sig = (x[i]-x[i-1])/(x[i+1] - x[i-1]);
458     p = sig * y2[i-1] + 2.0;
459     y2[i] = (sig - 1.0) / p;
460     u[i] = (6.0 * ((y[i+1] - y[i])/(x[i+1] - x[i]) -(y[i] - y[i-1]) / (x[i] - x[i-1]))/(x[i+1] - x[i-1]) - sig * u[i-1])/p;
461   }
462   if (ypn>0.99E30) {
463     qn = 0;
464     un = 0;
465   } else {
466     qn = 0.5;
467     un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])/(x[n-1] - x[n-2]));
468   }
469   y2[n-1] = (un - qn * u[n-2]) / (qn * y2[n-2] + 1.0);
470   for (k=n-2;k>=0;k--)
471     y2[k] = y2[k] * y2[k+1] + u[k];
472
473   delete u;
474 }
475
476 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
477
478 double MSIS::zeta(double zz, double zl)
479 {
480   return ((zz-zl)*(re+zl)/(re+zz));
481 }
482
483 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
484
485 double MSIS::densm(double alt, double d0, double xm, double *tz, int mn3,
486                      double *zn3, double *tn3, double *tgn3, int mn2, double *zn2,
487                      double *tn2, double *tgn2)
488 {
489 /*      Calculate Temperature and Density Profiles for lower atmos.  */
490   double xs[10] = {0,0,0,0,0,0,0,0,0,0};
491   double ys[10] = {0,0,0,0,0,0,0,0,0,0};
492   double y2out[10] = {0,0,0,0,0,0,0,0,0,0};
493   double rgas = 831.4;
494   double z=0, z1=0, z2=0, t1=0, t2=0, zg=0, zgdif=0;
495   double yd1=0, yd2=0;
496   double x=0, y=0, yi=0;
497   double expl=0, gamm=0, glb=0;
498   double densm_tmp=0;
499   int mn=0;
500   int k=0;
501   densm_tmp=d0;
502   if (alt>zn2[0]) {
503     if (xm==0.0)
504       return *tz;
505     else
506       return d0;
507   }
508
509   /* STRATOSPHERE/MESOSPHERE TEMPERATURE */
510   if (alt>zn2[mn2-1])
511     z=alt;
512   else
513     z=zn2[mn2-1];
514   mn=mn2;
515   z1=zn2[0];
516   z2=zn2[mn-1];
517   t1=tn2[0];
518   t2=tn2[mn-1];
519   zg = zeta(z, z1);
520   zgdif = zeta(z2, z1);
521
522   /* set up spline nodes */
523   for (k=0;k<mn;k++) {
524     xs[k]=zeta(zn2[k],z1)/zgdif;
525     ys[k]=1.0 / tn2[k];
526   }
527   yd1=-tgn2[0] / (t1*t1) * zgdif;
528   yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
529
530   /* calculate spline coefficients */
531   spline (xs, ys, mn, yd1, yd2, y2out);
532   x = zg/zgdif;
533   splint (xs, ys, y2out, mn, x, &y);
534
535   /* temperature at altitude */
536   *tz = 1.0 / y;
537   if (xm!=0.0) {
538     /* calaculate stratosphere / mesospehere density */
539     glb = gsurf / (pow((1.0 + z1/re),2.0));
540     gamm = xm * glb * zgdif / rgas;
541
542     /* Integrate temperature profile */
543     splini(xs, ys, y2out, mn, x, &yi);
544     expl=gamm*yi;
545     if (expl>50.0)
546       expl=50.0;
547
548     /* Density at altitude */
549     densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
550   }
551
552   if (alt>zn3[0]) {
553     if (xm==0.0)
554       return *tz;
555     else
556       return densm_tmp;
557   }
558
559   /* troposhere / stratosphere temperature */
560   z = alt;
561   mn = mn3;
562   z1=zn3[0];
563   z2=zn3[mn-1];
564   t1=tn3[0];
565   t2=tn3[mn-1];
566   zg=zeta(z,z1);
567   zgdif=zeta(z2,z1);
568
569   /* set up spline nodes */
570   for (k=0;k<mn;k++) {
571     xs[k] = zeta(zn3[k],z1) / zgdif;
572     ys[k] = 1.0 / tn3[k];
573   }
574   yd1=-tgn3[0] / (t1*t1) * zgdif;
575   yd2=-tgn3[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
576
577   /* calculate spline coefficients */
578   spline (xs, ys, mn, yd1, yd2, y2out);
579   x = zg/zgdif;
580   splint (xs, ys, y2out, mn, x, &y);
581
582   /* temperature at altitude */
583   *tz = 1.0 / y;
584   if (xm!=0.0) {
585     /* calaculate tropospheric / stratosphere density */
586     glb = gsurf / (pow((1.0 + z1/re),2.0));
587     gamm = xm * glb * zgdif / rgas;
588
589     /* Integrate temperature profile */
590     splini(xs, ys, y2out, mn, x, &yi);
591     expl=gamm*yi;
592     if (expl>50.0)
593       expl=50.0;
594
595     /* Density at altitude */
596     densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
597   }
598   if (xm==0.0)
599     return *tz;
600   else
601     return densm_tmp;
602 }
603
604 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
605
606 double MSIS::densu(double alt, double dlb, double tinf, double tlb, double xm,
607                      double alpha, double *tz, double zlb, double s2, int mn1,
608                      double *zn1, double *tn1, double *tgn1)
609 {
610 /*      Calculate Temperature and Density Profiles for MSIS models
611  *      New lower thermo polynomial
612  */
613   double yd2=0.0, yd1=0.0, x=0.0, y=0.0;
614   double rgas=831.4;
615   double densu_temp=1.0;
616   double za=0.0, z=0.0, zg2=0.0, tt=0.0, ta=0.0;
617   double dta=0.0, z1=0.0, z2=0.0, t1=0.0, t2=0.0, zg=0.0, zgdif=0.0;
618   int mn=0;
619   int k=0;
620   double glb=0.0;
621   double expl=0.0;
622   double yi=0.0;
623   double densa=0.0;
624   double gamma=0.0, gamm=0.0;
625   double xs[5]={0.0,0.0,0.0,0.0,0.0}, ys[5]={0.0,0.0,0.0,0.0,0.0}, y2out[5]={0.0,0.0,0.0,0.0,0.0};
626   /* joining altitudes of Bates and spline */
627   za=zn1[0];
628   if (alt>za)
629     z=alt;
630   else
631     z=za;
632
633   /* geopotential altitude difference from ZLB */
634   zg2 = zeta(z, zlb);
635
636   /* Bates temperature */
637   tt = tinf - (tinf - tlb) * exp(-s2*zg2);
638   ta = tt;
639   *tz = tt;
640   densu_temp = *tz;
641
642   if (alt<za) {
643     /* calculate temperature below ZA
644      * temperature gradient at ZA from Bates profile */
645     dta = (tinf - ta) * s2 * pow(((re+zlb)/(re+za)),2.0);
646     tgn1[0]=dta;
647     tn1[0]=ta;
648     if (alt>zn1[mn1-1])
649       z=alt;
650     else
651       z=zn1[mn1-1];
652     mn=mn1;
653     z1=zn1[0];
654     z2=zn1[mn-1];
655     t1=tn1[0];
656     t2=tn1[mn-1];
657     /* geopotental difference from z1 */
658     zg = zeta (z, z1);
659     zgdif = zeta(z2, z1);
660     /* set up spline nodes */
661     for (k=0;k<mn;k++) {
662       xs[k] = zeta(zn1[k], z1) / zgdif;
663       ys[k] = 1.0 / tn1[k];
664     }
665     /* end node derivatives */
666     yd1 = -tgn1[0] / (t1*t1) * zgdif;
667     yd2 = -tgn1[1] / (t2*t2) * zgdif * pow(((re+z2)/(re+z1)),2.0);
668     /* calculate spline coefficients */
669     spline (xs, ys, mn, yd1, yd2, y2out);
670     x = zg / zgdif;
671     splint (xs, ys, y2out, mn, x, &y);
672     /* temperature at altitude */
673     *tz = 1.0 / y;
674     densu_temp = *tz;
675   }
676   if (xm==0)
677     return densu_temp;
678
679   /* calculate density above za */
680   glb = gsurf / pow((1.0 + zlb/re),2.0);
681   gamma = xm * glb / (s2 * rgas * tinf);
682   expl = exp(-s2 * gamma * zg2);
683   if (expl>50.0)
684       expl=50.0;
685   if (tt<=0)
686     expl=50.0;
687
688   /* density at altitude */
689   densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl;
690   densu_temp=densa;
691   if (alt>=za)
692     return densu_temp;
693
694   /* calculate density below za */
695   glb = gsurf / pow((1.0 + z1/re),2.0);
696   gamm = xm * glb * zgdif / rgas;
697
698   /* integrate spline temperatures */
699   splini (xs, ys, y2out, mn, x, &yi);
700   expl = gamm * yi;
701   if (expl>50.0)
702     expl=50.0;
703   if (*tz<=0)
704     expl=50.0;
705
706   /* density at altitude */
707   densu_temp = densu_temp * pow ((t1 / *tz),(1.0 + alpha)) * exp(-expl);
708   return densu_temp;
709 }
710
711 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
712
713 /*    3hr Magnetic activity functions */
714 /*    Eq. A24d */
715 double MSIS::g0(double a, double *p)
716 {
717   return (a - 4.0 + (p[25] - 1.0) * (a - 4.0 + (exp(-sqrt(p[24]*p[24]) *
718                 (a - 4.0)) - 1.0) / sqrt(p[24]*p[24])));
719 }
720
721 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
722
723 /*    Eq. A24c */
724 double MSIS::sumex(double ex)
725 {
726   return (1.0 + (1.0 - pow(ex,19.0)) / (1.0 - ex) * pow(ex,0.5));
727 }
728
729 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
730
731 /*    Eq. A24a */
732 double MSIS::sg0(double ex, double *p, double *ap)
733 {
734   return (g0(ap[1],p) + (g0(ap[2],p)*ex + g0(ap[3],p)*ex*ex +
735                 g0(ap[4],p)*pow(ex,3.0)  + (g0(ap[5],p)*pow(ex,4.0) +
736                 g0(ap[6],p)*pow(ex,12.0))*(1.0-pow(ex,8.0))/(1.0-ex)))/sumex(ex);
737 }
738
739 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
740
741 double MSIS::globe7(double *p, struct nrlmsise_input *input,
742                       struct nrlmsise_flags *flags)
743 {
744 /*       CALCULATE G(L) FUNCTION
745  *       Upper Thermosphere Parameters */
746   double t[15];
747   int i,j;
748   int sw9=1;
749   double apd;
750   double xlong;
751   double tloc;
752   double c, s, c2, c4, s2;
753   double sr = 7.2722E-5;
754   double dgtr = 1.74533E-2;
755   double dr = 1.72142E-2;
756   double hr = 0.2618;
757   double cd32, cd18, cd14, cd39;
758   double p32, p18, p14, p39;
759   double df;
760   double f1, f2;
761   double tinf;
762   struct ap_array *ap;
763
764   tloc=input->lst;
765   for (j=0;j<14;j++)
766     t[j]=0;
767   if (flags->sw[9]>0)
768     sw9=1;
769   else if (flags->sw[9]<0)
770     sw9=-1;
771   xlong = input->g_long;
772
773   /* calculate legendre polynomials */
774   c = sin(input->g_lat * dgtr);
775   s = cos(input->g_lat * dgtr);
776   c2 = c*c;
777   c4 = c2*c2;
778   s2 = s*s;
779
780   plg[0][1] = c;
781   plg[0][2] = 0.5*(3.0*c2 -1.0);
782   plg[0][3] = 0.5*(5.0*c*c2-3.0*c);
783   plg[0][4] = (35.0*c4 - 30.0*c2 + 3.0)/8.0;
784   plg[0][5] = (63.0*c2*c2*c - 70.0*c2*c + 15.0*c)/8.0;
785   plg[0][6] = (11.0*c*plg[0][5] - 5.0*plg[0][4])/6.0;
786 /*      plg[0][7] = (13.0*c*plg[0][6] - 6.0*plg[0][5])/7.0; */
787   plg[1][1] = s;
788   plg[1][2] = 3.0*c*s;
789   plg[1][3] = 1.5*(5.0*c2-1.0)*s;
790   plg[1][4] = 2.5*(7.0*c2*c-3.0*c)*s;
791   plg[1][5] = 1.875*(21.0*c4 - 14.0*c2 +1.0)*s;
792   plg[1][6] = (11.0*c*plg[1][5]-6.0*plg[1][4])/5.0;
793 /*      plg[1][7] = (13.0*c*plg[1][6]-7.0*plg[1][5])/6.0; */
794 /*      plg[1][8] = (15.0*c*plg[1][7]-8.0*plg[1][6])/7.0; */
795   plg[2][2] = 3.0*s2;
796   plg[2][3] = 15.0*s2*c;
797   plg[2][4] = 7.5*(7.0*c2 -1.0)*s2;
798   plg[2][5] = 3.0*c*plg[2][4]-2.0*plg[2][3];
799   plg[2][6] =(11.0*c*plg[2][5]-7.0*plg[2][4])/4.0;
800   plg[2][7] =(13.0*c*plg[2][6]-8.0*plg[2][5])/5.0;
801   plg[3][3] = 15.0*s2*s;
802   plg[3][4] = 105.0*s2*s*c;
803   plg[3][5] =(9.0*c*plg[3][4]-7.*plg[3][3])/2.0;
804   plg[3][6] =(11.0*c*plg[3][5]-8.*plg[3][4])/3.0;
805
806   if (!(((flags->sw[7]==0)&&(flags->sw[8]==0))&&(flags->sw[14]==0))) {
807     stloc = sin(hr*tloc);
808     ctloc = cos(hr*tloc);
809     s2tloc = sin(2.0*hr*tloc);
810     c2tloc = cos(2.0*hr*tloc);
811     s3tloc = sin(3.0*hr*tloc);
812     c3tloc = cos(3.0*hr*tloc);
813   }
814
815   cd32 = cos(dr*(input->doy-p[31]));
816   cd18 = cos(2.0*dr*(input->doy-p[17]));
817   cd14 = cos(dr*(input->doy-p[13]));
818   cd39 = cos(2.0*dr*(input->doy-p[38]));
819   p32=p[31];
820   p18=p[17];
821   p14=p[13];
822   p39=p[38];
823
824   /* F10.7 EFFECT */
825   df = input->f107 - input->f107A;
826   dfa = input->f107A - 150.0;
827   t[0] =  p[19]*df*(1.0+p[59]*dfa) + p[20]*df*df + p[21]*dfa + p[29]*pow(dfa,2.0);
828   f1 = 1.0 + (p[47]*dfa +p[19]*df+p[20]*df*df)*flags->swc[1];
829   f2 = 1.0 + (p[49]*dfa+p[19]*df+p[20]*df*df)*flags->swc[1];
830
831   /*  TIME INDEPENDENT */
832   t[1] = (p[1]*plg[0][2]+ p[2]*plg[0][4]+p[22]*plg[0][6]) +
833         (p[14]*plg[0][2])*dfa*flags->swc[1] +p[26]*plg[0][1];
834
835   /*  SYMMETRICAL ANNUAL */
836   t[2] = p[18]*cd32;
837
838   /*  SYMMETRICAL SEMIANNUAL */
839   t[3] = (p[15]+p[16]*plg[0][2])*cd18;
840
841   /*  ASYMMETRICAL ANNUAL */
842   t[4] =  f1*(p[9]*plg[0][1]+p[10]*plg[0][3])*cd14;
843
844   /*  ASYMMETRICAL SEMIANNUAL */
845   t[5] =    p[37]*plg[0][1]*cd39;
846
847         /* DIURNAL */
848   if (flags->sw[7]) {
849     double t71, t72;
850     t71 = (p[11]*plg[1][2])*cd14*flags->swc[5];
851     t72 = (p[12]*plg[1][2])*cd14*flags->swc[5];
852     t[6] = f2*((p[3]*plg[1][1] + p[4]*plg[1][3] + p[27]*plg[1][5] + t71) * \
853          ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + p[28]*plg[1][5] \
854             + t72)*stloc);
855 }
856
857   /* SEMIDIURNAL */
858   if (flags->sw[8]) {
859     double t81, t82;
860     t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
861     t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
862     t[7] = f2*((p[5]*plg[2][2]+ p[41]*plg[2][4] + t81)*c2tloc +(p[8]*plg[2][2] + p[42]*plg[2][4] + t82)*s2tloc);
863   }
864
865   /* TERDIURNAL */
866   if (flags->sw[14]) {
867     t[13] = f2 * ((p[39]*plg[3][3]+(p[93]*plg[3][4]+p[46]*plg[3][6])*cd14*flags->swc[5])* s3tloc +(p[40]*plg[3][3]+(p[94]*plg[3][4]+p[48]*plg[3][6])*cd14*flags->swc[5])* c3tloc);
868 }
869
870   /* magnetic activity based on daily ap */
871   if (flags->sw[9]==-1) {
872     ap = input->ap_a;
873     if (p[51]!=0) {
874       double exp1;
875       exp1 = exp(-10800.0*sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-sqrt(input->g_lat*input->g_lat))));
876       if (exp1>0.99999)
877         exp1=0.99999;
878       if (p[24]<1.0E-4)
879         p[24]=1.0E-4;
880       apt[0]=sg0(exp1,p,ap->a);
881       /* apt[1]=sg2(exp1,p,ap->a);
882          apt[2]=sg0(exp2,p,ap->a);
883          apt[3]=sg2(exp2,p,ap->a);
884       */
885       if (flags->sw[9]) {
886         t[8] = apt[0]*(p[50]+p[96]*plg[0][2]+p[54]*plg[0][4]+ \
887      (p[125]*plg[0][1]+p[126]*plg[0][3]+p[127]*plg[0][5])*cd14*flags->swc[5]+ \
888      (p[128]*plg[1][1]+p[129]*plg[1][3]+p[130]*plg[1][5])*flags->swc[7]* \
889                  cos(hr*(tloc-p[131])));
890       }
891     }
892   } else {
893     double p44, p45;
894     apd=input->ap-4.0;
895     p44=p[43];
896     p45=p[44];
897     if (p44<0)
898       p44 = 1.0E-5;
899     apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44);
900     if (flags->sw[9]) {
901       t[8]=apdf*(p[32]+p[45]*plg[0][2]+p[34]*plg[0][4]+ \
902      (p[100]*plg[0][1]+p[101]*plg[0][3]+p[102]*plg[0][5])*cd14*flags->swc[5]+
903      (p[121]*plg[1][1]+p[122]*plg[1][3]+p[123]*plg[1][5])*flags->swc[7]*
904             cos(hr*(tloc-p[124])));
905     }
906   }
907
908   if ((flags->sw[10])&&(input->g_long>-1000.0)) {
909
910     /* longitudinal */
911     if (flags->sw[11]) {
912       t[10] = (1.0 + p[80]*dfa*flags->swc[1])* \
913      ((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
914       +p[103]*plg[1][1]+p[104]*plg[1][3]+p[105]*plg[1][5]\
915       +flags->swc[5]*(p[109]*plg[1][1]+p[110]*plg[1][3]+p[111]*plg[1][5])*cd14)* \
916           cos(dgtr*input->g_long) \
917       +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
918       +p[106]*plg[1][1]+p[107]*plg[1][3]+p[108]*plg[1][5]\
919       +flags->swc[5]*(p[112]*plg[1][1]+p[113]*plg[1][3]+p[114]*plg[1][5])*cd14)* \
920       sin(dgtr*input->g_long));
921     }
922
923     /* ut and mixed ut, longitude */
924     if (flags->sw[12]){
925       t[11]=(1.0+p[95]*plg[0][1])*(1.0+p[81]*dfa*flags->swc[1])*\
926         (1.0+p[119]*plg[0][1]*flags->swc[5]*cd14)*\
927         ((p[68]*plg[0][1]+p[69]*plg[0][3]+p[70]*plg[0][5])*\
928         cos(sr*(input->sec-p[71])));
929       t[11]+=flags->swc[11]*\
930         (p[76]*plg[2][3]+p[77]*plg[2][5]+p[78]*plg[2][7])*\
931         cos(sr*(input->sec-p[79])+2.0*dgtr*input->g_long)*(1.0+p[137]*dfa*flags->swc[1]);
932     }
933
934     /* ut, longitude magnetic activity */
935     if (flags->sw[13]) {
936       if (flags->sw[9]==-1) {
937         if (p[51]) {
938           t[12]=apt[0]*flags->swc[11]*(1.+p[132]*plg[0][1])*\
939             ((p[52]*plg[1][2]+p[98]*plg[1][4]+p[67]*plg[1][6])*\
940              cos(dgtr*(input->g_long-p[97])))\
941             +apt[0]*flags->swc[11]*flags->swc[5]*\
942             (p[133]*plg[1][1]+p[134]*plg[1][3]+p[135]*plg[1][5])*\
943             cd14*cos(dgtr*(input->g_long-p[136])) \
944             +apt[0]*flags->swc[12]* \
945             (p[55]*plg[0][1]+p[56]*plg[0][3]+p[57]*plg[0][5])*\
946             cos(sr*(input->sec-p[58]));
947         }
948       } else {
949         t[12] = apdf*flags->swc[11]*(1.0+p[120]*plg[0][1])*\
950           ((p[60]*plg[1][2]+p[61]*plg[1][4]+p[62]*plg[1][6])*\
951           cos(dgtr*(input->g_long-p[63])))\
952           +apdf*flags->swc[11]*flags->swc[5]* \
953           (p[115]*plg[1][1]+p[116]*plg[1][3]+p[117]*plg[1][5])* \
954           cd14*cos(dgtr*(input->g_long-p[118])) \
955           + apdf*flags->swc[12]* \
956           (p[83]*plg[0][1]+p[84]*plg[0][3]+p[85]*plg[0][5])* \
957           cos(sr*(input->sec-p[75]));
958       }
959     }
960   }
961
962   /* parms not used: 82, 89, 99, 139-149 */
963   tinf = p[30];
964   for (i=0;i<14;i++)
965     tinf = tinf + fabs(flags->sw[i+1])*t[i];
966   return tinf;
967 }
968
969 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
970
971 double MSIS::glob7s(double *p, struct nrlmsise_input *input,
972                       struct nrlmsise_flags *flags)
973 {
974 /*    VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99
975  */
976   double pset=2.0;
977   double t[14] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,};
978   double tt=0.0;
979   double cd32=0.0, cd18=0.0, cd14=0.0, cd39=0.0;
980   double p32=0.0, p18=0.0, p14=0.0, p39=0.0;
981   int i=0,j=0;
982   double dr=1.72142E-2;
983   double dgtr=1.74533E-2;
984   /* confirm parameter set */
985   if (p[99]==0)
986     p[99]=pset;
987   if (p[99]!=pset) {
988     cerr << "Wrong parameter set for glob7s" << endl;
989     return -1;
990   }
991   for (j=0;j<14;j++)
992     t[j]=0.0;
993   cd32 = cos(dr*(input->doy-p[31]));
994   cd18 = cos(2.0*dr*(input->doy-p[17]));
995   cd14 = cos(dr*(input->doy-p[13]));
996   cd39 = cos(2.0*dr*(input->doy-p[38]));
997   p32=p[31];
998   p18=p[17];
999   p14=p[13];
1000   p39=p[38];
1001
1002   /* F10.7 */
1003   t[0] = p[21]*dfa;
1004
1005   /* time independent */
1006   t[1]=p[1]*plg[0][2] + p[2]*plg[0][4] + p[22]*plg[0][6] + p[26]*plg[0][1] + p[14]*plg[0][3] + p[59]*plg[0][5];
1007
1008         /* SYMMETRICAL ANNUAL */
1009   t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32;
1010
1011         /* SYMMETRICAL SEMIANNUAL */
1012   t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18;
1013
1014         /* ASYMMETRICAL ANNUAL */
1015   t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14;
1016
1017   /* ASYMMETRICAL SEMIANNUAL */
1018   t[5]=(p[37]*plg[0][1])*cd39;
1019
1020         /* DIURNAL */
1021   if (flags->sw[7]) {
1022     double t71, t72;
1023     t71 = p[11]*plg[1][2]*cd14*flags->swc[5];
1024     t72 = p[12]*plg[1][2]*cd14*flags->swc[5];
1025     t[6] = ((p[3]*plg[1][1] + p[4]*plg[1][3] + t71) * ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + t72) * stloc) ;
1026   }
1027
1028   /* SEMIDIURNAL */
1029   if (flags->sw[8]) {
1030     double t81, t82;
1031     t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
1032     t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
1033     t[7] = ((p[5]*plg[2][2] + p[41]*plg[2][4] + t81) * c2tloc + (p[8]*plg[2][2] + p[42]*plg[2][4] + t82) * s2tloc);
1034   }
1035
1036   /* TERDIURNAL */
1037   if (flags->sw[14]) {
1038     t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
1039   }
1040
1041   /* MAGNETIC ACTIVITY */
1042   if (flags->sw[9]) {
1043     if (flags->sw[9]==1)
1044       t[8] = apdf * (p[32] + p[45] * plg[0][2] * flags->swc[2]);
1045     if (flags->sw[9]==-1)
1046       t[8]=(p[50]*apt[0] + p[96]*plg[0][2] * apt[0]*flags->swc[2]);
1047   }
1048
1049   /* LONGITUDINAL */
1050   if (!((flags->sw[10]==0) || (flags->sw[11]==0) || (input->g_long<=-1000.0))) {
1051     t[10] = (1.0 + plg[0][1]*(p[80]*flags->swc[5]*cos(dr*(input->doy-p[81]))\
1052             +p[85]*flags->swc[6]*cos(2.0*dr*(input->doy-p[86])))\
1053       +p[83]*flags->swc[3]*cos(dr*(input->doy-p[84]))\
1054       +p[87]*flags->swc[4]*cos(2.0*dr*(input->doy-p[88])))\
1055       *((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
1056       +p[74]*plg[1][1]+p[75]*plg[1][3]+p[76]*plg[1][5]\
1057       )*cos(dgtr*input->g_long)\
1058       +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
1059       +p[77]*plg[1][1]+p[78]*plg[1][3]+p[79]*plg[1][5]\
1060       )*sin(dgtr*input->g_long));
1061   }
1062   tt=0;
1063   for (i=0;i<14;i++)
1064     tt+=fabs(flags->sw[i+1])*t[i];
1065   return tt;
1066 }
1067
1068 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1069
1070 void MSIS::gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1071                   struct nrlmsise_output *output)
1072 {
1073   double xlat=0.0;
1074   double xmm=0.0;
1075   int mn3 = 5;
1076   double zn3[5]={32.5,20.0,15.0,10.0,0.0};
1077   int mn2 = 4;
1078   double zn2[4]={72.5,55.0,45.0,32.5};
1079   double altt=0.0;
1080   double zmix=62.5;
1081   double tmp=0.0;
1082   double dm28m=0.0;
1083   double tz=0.0;
1084   double dmc=0.0;
1085   double dmr=0.0;
1086   double dz28=0.0;
1087   struct nrlmsise_output soutput;
1088   int i;
1089
1090   for (int i=0; i<9; i++) soutput.d[i] = 0.0;
1091   for (int i=0; i<2; i++) soutput.t[i] = 0.0;
1092
1093   tselec(flags);
1094
1095   /* Latitude variation of gravity (none for sw[2]=0) */
1096   xlat=input->g_lat;
1097   if (flags->sw[2]==0)
1098     xlat=45.0;
1099   glatf(xlat, &gsurf, &re);
1100
1101   xmm = pdm[2][4];
1102
1103   /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */
1104   if (input->alt>zn2[0])
1105     altt=input->alt;
1106   else
1107     altt=zn2[0];
1108
1109   tmp=input->alt;
1110   input->alt=altt;
1111   gts7(input, flags, &soutput);
1112   altt=input->alt;
1113   input->alt=tmp;
1114   if (flags->sw[0])   /* metric adjustment */
1115     dm28m=dm28*1.0E6;
1116   else
1117     dm28m=dm28;
1118   output->t[0]=soutput.t[0];
1119   output->t[1]=soutput.t[1];
1120   if (input->alt>=zn2[0]) {
1121     for (i=0;i<9;i++)
1122       output->d[i]=soutput.d[i];
1123     return;
1124   }
1125
1126 /*       LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])
1127  *         Temperature at nodes and gradients at end nodes
1128  *         Inverse temperature a linear function of spherical harmonics
1129  */
1130   meso_tgn2[0]=meso_tgn1[1];
1131   meso_tn2[0]=meso_tn1[4];
1132         meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags));
1133         meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags));
1134         meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags));
1135   meso_tgn2[1]=pavgm[8]*pma[9][0]*(1.0+flags->sw[20]*flags->sw[22]*glob7s(pma[9], input, flags))*meso_tn2[3]*meso_tn2[3]/(pow((pma[2][0]*pavgm[2]),2.0));
1136   meso_tn3[0]=meso_tn2[3];
1137
1138   if (input->alt<zn3[0]) {
1139 /*       LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])
1140  *         Temperature at nodes and gradients at end nodes
1141  *         Inverse temperature a linear function of spherical harmonics
1142  */
1143     meso_tgn3[0]=meso_tgn2[1];
1144     meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags->sw[22]*glob7s(pma[3], input, flags));
1145     meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags));
1146     meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags));
1147     meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags));
1148     meso_tgn3[1]=pma[7][0]*pavgm[7]*(1.0+flags->sw[22]*glob7s(pma[7], input, flags)) *meso_tn3[4]*meso_tn3[4]/(pow((pma[6][0]*pavgm[6]),2.0));
1149   }
1150
1151         /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */
1152
1153   dmc=0;
1154   if (input->alt>zmix)
1155     dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix);
1156   dz28=soutput.d[2];
1157
1158   /**** N2 density ****/
1159   dmr=soutput.d[2] / dm28m - 1.0;
1160   output->d[2]=densm(input->alt,dm28m,xmm, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
1161   output->d[2]=output->d[2] * (1.0 + dmr*dmc);
1162
1163   /**** HE density ****/
1164   dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0;
1165   output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc);
1166
1167   /**** O density ****/
1168   output->d[1] = 0;
1169   output->d[8] = 0;
1170
1171   /**** O2 density ****/
1172   dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0;
1173   output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc);
1174
1175   /**** AR density ***/
1176   dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0;
1177   output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc);
1178
1179   /**** Hydrogen density ****/
1180   output->d[6] = 0;
1181
1182   /**** Atomic nitrogen density ****/
1183   output->d[7] = 0;
1184
1185   /**** Total mass density */
1186   output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] +
1187                      28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4]
1188                      + output->d[6] + 14.0 * output->d[7]);
1189
1190   if (flags->sw[0])
1191     output->d[5]=output->d[5]/1000;
1192
1193   /**** temperature at altitude ****/
1194   dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3,
1195                    mn2, zn2, meso_tn2, meso_tgn2);
1196   output->t[1]=tz;
1197
1198 }
1199
1200 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1201
1202 void MSIS::gtd7d(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1203                    struct nrlmsise_output *output)
1204 {
1205   gtd7(input, flags, output);
1206   output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] +
1207                    28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4]
1208                    + output->d[6] + 14.0 * output->d[7] + 16.0 * output->d[8]);
1209   if (flags->sw[0])
1210     output->d[5]=output->d[5]/1000;
1211 }
1212
1213 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1214
1215 void MSIS::ghp7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1216                   struct nrlmsise_output *output, double press)
1217 {
1218   double bm = 1.3806E-19;
1219   double rgas = 831.4;
1220   double test = 0.00043;
1221   double ltest = 12;
1222   double pl, p;
1223   double zi = 0.0;
1224   double z;
1225   double cl, cl2;
1226   double ca, cd;
1227   double xn, xm, diff;
1228   double g, sh;
1229   int l;
1230   pl = log10(press);
1231   if (pl >= -5.0) {
1232     if (pl>2.5)
1233       zi = 18.06 * (3.00 - pl);
1234     else if ((pl>0.075) && (pl<=2.5))
1235       zi = 14.98 * (3.08 - pl);
1236     else if ((pl>-1) && (pl<=0.075))
1237       zi = 17.80 * (2.72 - pl);
1238     else if ((pl>-2) && (pl<=-1))
1239       zi = 14.28 * (3.64 - pl);
1240     else if ((pl>-4) && (pl<=-2))
1241       zi = 12.72 * (4.32 -pl);
1242     else if (pl<=-4)
1243       zi = 25.3 * (0.11 - pl);
1244     cl = input->g_lat/90.0;
1245     cl2 = cl*cl;
1246     if (input->doy<182)
1247       cd = (1.0 - (double) input->doy) / 91.25;
1248     else
1249       cd = ((double) input->doy) / 91.25 - 3.0;
1250     ca = 0;
1251     if ((pl > -1.11) && (pl<=-0.23))
1252       ca = 1.0;
1253     if (pl > -0.23)
1254       ca = (2.79 - pl) / (2.79 + 0.23);
1255     if ((pl <= -1.11) && (pl>-3))
1256       ca = (-2.93 - pl)/(-2.93 + 1.11);
1257     z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl;
1258   } else
1259     z = 22.0 * pow((pl + 4.0),2.0) + 110.0;
1260
1261   /* iteration  loop */
1262   l = 0;
1263   do {
1264     l++;
1265     input->alt = z;
1266     gtd7(input, flags, output);
1267     z = input->alt;
1268     xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7];
1269     p = bm * xn * output->t[1];
1270     if (flags->sw[0])
1271       p = p*1.0E-6;
1272     diff = pl - log10(p);
1273     if (sqrt(diff*diff)<test)
1274       return;
1275     if (l==ltest) {
1276       cerr << "ERROR: ghp7 not converging for press " << press << ", diff " << diff << endl;
1277       return;
1278     }
1279     xm = output->d[5] / xn / 1.66E-24;
1280     if (flags->sw[0])
1281       xm = xm * 1.0E3;
1282     g = gsurf / (pow((1.0 + z/re),2.0));
1283     sh = rgas * output->t[1] / (xm * g);
1284
1285     /* new altitude estimate using scale height */
1286     if (l <  6)
1287       z = z - sh * diff * 2.302;
1288     else
1289       z = z - sh * diff;
1290   } while (1==1);
1291 }
1292
1293 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1294
1295 void MSIS::gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1296                   struct nrlmsise_output *output)
1297 {
1298 /*     Thermospheric portion of NRLMSISE-00
1299  *     See GTD7 for more extensive comments
1300  *     alt > 72.5 km!
1301  */
1302   double za=0.0;
1303   int i, j;
1304   double ddum=0.0, z=0.0;
1305   double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
1306   double tinf=0.0;
1307   int mn1 = 5;
1308   double g0=0.0;
1309   double tlb=0.0;
1310   double s=0.0, z0=0.0, t0=0.0, tr12=0.0;
1311   double db01=0.0, db04=0.0, db14=0.0, db16=0.0, db28=0.0, db32=0.0, db40=0.0, db48=0.0;
1312   double zh28=0.0, zh04=0.0, zh16=0.0, zh32=0.0, zh40=0.0, zh01=0.0, zh14=0.0;
1313   double zhm28=0.0, zhm04=0.0, zhm16=0.0, zhm32=0.0, zhm40=0.0, zhm01=0.0, zhm14=0.0;
1314   double xmd=0.0;
1315   double b28=0.0, b04=0.0, b16=0.0, b32=0.0, b40=0.0, b01=0.0, b14=0.0;
1316   double tz=0.0;
1317   double g28=0.0, g4=0.0, g16=0.0, g32=0.0, g40=0.0, g1=0.0, g14=0.0;
1318   double zhf=0.0, xmm=0.0;
1319   double zc04=0.0, zc16=0.0, zc32=0.0, zc40=0.0, zc01=0.0, zc14=0.0;
1320   double hc04=0.0, hc16=0.0, hc32=0.0, hc40=0.0, hc01=0.0, hc14=0.0;
1321   double hcc16=0.0, hcc32=0.0, hcc01=0.0, hcc14=0.0;
1322   double zcc16=0.0, zcc32=0.0, zcc01=0.0, zcc14=0.0;
1323   double rc16=0.0, rc32=0.0, rc01=0.0, rc14=0.0;
1324   double rl=0.0;
1325   double g16h=0.0, db16h=0.0, tho=0.0, zsht=0.0, zmho=0.0, zsho=0.0;
1326   double dgtr=1.74533E-2;
1327   double dr=1.72142E-2;
1328   double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
1329   double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
1330   double dd=0.0;
1331   double hc216=0.0, hcc232=0.0;
1332   za = pdl[1][15];
1333   zn1[0] = za;
1334   for (j=0;j<9;j++)
1335     output->d[j]=0;
1336
1337   /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */
1338   if (input->alt>zn1[0])
1339     tinf = ptm[0]*pt[0] * \
1340       (1.0+flags->sw[16]*globe7(pt,input,flags));
1341   else
1342     tinf = ptm[0]*pt[0];
1343   output->t[0]=tinf;
1344
1345   /*  GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */
1346   if (input->alt>zn1[4])
1347     g0 = ptm[3]*ps[0] * \
1348       (1.0+flags->sw[19]*globe7(ps,input,flags));
1349   else
1350     g0 = ptm[3]*ps[0];
1351   tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0];
1352   s = g0 / (tinf - tlb);
1353
1354 /*      Lower thermosphere temp variations not significant for
1355  *       density above 300 km */
1356   if (input->alt<300.0) {
1357     meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags));
1358     meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags));
1359     meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags));
1360     meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags));
1361     meso_tgn1[1]=ptm[8]*pma[8][0]*(1.0+flags->sw[18]*flags->sw[20]*glob7s(pma[8], input, flags))*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
1362   } else {
1363     meso_tn1[1]=ptm[6]*ptl[0][0];
1364     meso_tn1[2]=ptm[2]*ptl[1][0];
1365     meso_tn1[3]=ptm[7]*ptl[2][0];
1366     meso_tn1[4]=ptm[4]*ptl[3][0];
1367     meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
1368   }
1369
1370   z0 = zn1[3];
1371   t0 = meso_tn1[3];
1372   tr12 = 1.0;
1373
1374   /* N2 variation factor at Zlb */
1375   g28=flags->sw[21]*globe7(pd[2], input, flags);
1376
1377   /* VARIATION OF TURBOPAUSE HEIGHT */
1378   zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13])));
1379   output->t[0]=tinf;
1380   xmm = pdm[2][4];
1381   z = input->alt;
1382
1383
1384         /**** N2 DENSITY ****/
1385
1386   /* Diffusive density at Zlb */
1387   db28 = pdm[2][0]*exp(g28)*pd[2][0];
1388   /* Diffusive density at Alt */
1389   output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1390   dd=output->d[2];
1391   /* Turbopause */
1392   zh28=pdm[2][2]*zhf;
1393   zhm28=pdm[2][3]*pdl[1][5];
1394   xmd=28.0-xmm;
1395   /* Mixed density at Zlb */
1396   b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1397   if ((flags->sw[15])&&(z<=altl[2])) {
1398     /*  Mixed density at Alt */
1399     dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1400     /*  Net density at Alt */
1401     output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0);
1402   }
1403
1404
1405         /**** HE DENSITY ****/
1406
1407   /*   Density variation factor at Zlb */
1408   g4 = flags->sw[21]*globe7(pd[0], input, flags);
1409   /*  Diffusive density at Zlb */
1410   db04 = pdm[0][0]*exp(g4)*pd[0][0];
1411         /*  Diffusive density at Alt */
1412   output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1413   dd=output->d[0];
1414   if ((flags->sw[15]) && (z<altl[0])) {
1415     /*  Turbopause */
1416     zh04=pdm[0][2];
1417     /*  Mixed density at Zlb */
1418     b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1419     /*  Mixed density at Alt */
1420     dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1421     zhm04=zhm28;
1422     /*  Net density at Alt */
1423     output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.);
1424     /*  Correction to specified mixing ratio at ground */
1425     rl=log(b28*pdm[0][1]/b04);
1426     zc04=pdm[0][4]*pdl[1][0];
1427     hc04=pdm[0][5]*pdl[1][1];
1428     /*  Net density corrected at Alt */
1429     output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
1430   }
1431
1432
1433         /**** O DENSITY ****/
1434
1435   /*  Density variation factor at Zlb */
1436   g16= flags->sw[21]*globe7(pd[1],input,flags);
1437   /*  Diffusive density at Zlb */
1438   db16 =  pdm[1][0]*exp(g16)*pd[1][0];
1439         /*   Diffusive density at Alt */
1440   output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1441   dd=output->d[1];
1442   if ((flags->sw[15]) && (z<=altl[1])) {
1443     /*   Turbopause */
1444     zh16=pdm[1][2];
1445     /*  Mixed density at Zlb */
1446     b16=densu(zh16,db16,tinf,tlb,16.0-xmm,(alpha[1]-1.0), &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1447     /*  Mixed density at Alt */
1448     dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1449     zhm16=zhm28;
1450     /*  Net density at Alt */
1451     output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.);
1452     rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0));
1453     hc16=pdm[1][5]*pdl[1][3];
1454     zc16=pdm[1][4]*pdl[1][2];
1455     hc216=pdm[1][5]*pdl[1][4];
1456     output->d[1]=output->d[1]*ccor2(z,rl,hc16,zc16,hc216);
1457     /*   Chemistry correction */
1458     hcc16=pdm[1][7]*pdl[1][13];
1459     zcc16=pdm[1][6]*pdl[1][12];
1460     rc16=pdm[1][3]*pdl[1][14];
1461     /*  Net density corrected at Alt */
1462     output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16);
1463   }
1464
1465
1466         /**** O2 DENSITY ****/
1467
1468         /*   Density variation factor at Zlb */
1469   g32= flags->sw[21]*globe7(pd[4], input, flags);
1470         /*  Diffusive density at Zlb */
1471   db32 = pdm[3][0]*exp(g32)*pd[4][0];
1472         /*   Diffusive density at Alt */
1473   output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1474   dd=output->d[3];
1475   if (flags->sw[15]) {
1476     if (z<=altl[3]) {
1477       /*   Turbopause */
1478       zh32=pdm[3][2];
1479       /*  Mixed density at Zlb */
1480       b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1481       /*  Mixed density at Alt */
1482       dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1483       zhm32=zhm28;
1484       /*  Net density at Alt */
1485       output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.);
1486       /*   Correction to specified mixing ratio at ground */
1487       rl=log(b28*pdm[3][1]/b32);
1488       hc32=pdm[3][5]*pdl[1][7];
1489       zc32=pdm[3][4]*pdl[1][6];
1490       output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32);
1491     }
1492     /*  Correction for general departure from diffusive equilibrium above Zlb */
1493     hcc32=pdm[3][7]*pdl[1][22];
1494     hcc232=pdm[3][7]*pdl[0][22];
1495     zcc32=pdm[3][6]*pdl[1][21];
1496     rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.));
1497     /*  Net density corrected at Alt */
1498     output->d[3]=output->d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
1499   }
1500
1501
1502         /**** AR DENSITY ****/
1503
1504         /*   Density variation factor at Zlb */
1505   g40= flags->sw[21]*globe7(pd[5],input,flags);
1506         /*  Diffusive density at Zlb */
1507   db40 = pdm[4][0]*exp(g40)*pd[5][0];
1508   /*   Diffusive density at Alt */
1509   output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1510   dd=output->d[4];
1511   if ((flags->sw[15]) && (z<=altl[4])) {
1512     /*   Turbopause */
1513     zh40=pdm[4][2];
1514     /*  Mixed density at Zlb */
1515     b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1516     /*  Mixed density at Alt */
1517     dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1518     zhm40=zhm28;
1519     /*  Net density at Alt */
1520     output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.);
1521     /*   Correction to specified mixing ratio at ground */
1522     rl=log(b28*pdm[4][1]/b40);
1523     hc40=pdm[4][5]*pdl[1][9];
1524     zc40=pdm[4][4]*pdl[1][8];
1525     /*  Net density corrected at Alt */
1526     output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40);
1527     }
1528
1529
1530         /**** HYDROGEN DENSITY ****/
1531
1532         /*   Density variation factor at Zlb */
1533   g1 = flags->sw[21]*globe7(pd[6], input, flags);
1534         /*  Diffusive density at Zlb */
1535   db01 = pdm[5][0]*exp(g1)*pd[6][0];
1536         /*   Diffusive density at Alt */
1537   output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1538   dd=output->d[6];
1539   if ((flags->sw[15]) && (z<=altl[6])) {
1540     /*   Turbopause */
1541     zh01=pdm[5][2];
1542     /*  Mixed density at Zlb */
1543     b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1544     /*  Mixed density at Alt */
1545     dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1546     zhm01=zhm28;
1547     /*  Net density at Alt */
1548     output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.);
1549     /*   Correction to specified mixing ratio at ground */
1550     rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01);
1551     hc01=pdm[5][5]*pdl[1][11];
1552     zc01=pdm[5][4]*pdl[1][10];
1553     output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01);
1554     /*   Chemistry correction */
1555     hcc01=pdm[5][7]*pdl[1][19];
1556     zcc01=pdm[5][6]*pdl[1][18];
1557     rc01=pdm[5][3]*pdl[1][20];
1558     /*  Net density corrected at Alt */
1559     output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01);
1560 }
1561
1562
1563         /**** ATOMIC NITROGEN DENSITY ****/
1564
1565   /*   Density variation factor at Zlb */
1566   g14 = flags->sw[21]*globe7(pd[7],input,flags);
1567         /*  Diffusive density at Zlb */
1568   db14 = pdm[6][0]*exp(g14)*pd[7][0];
1569         /*   Diffusive density at Alt */
1570   output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1571   dd=output->d[7];
1572   if ((flags->sw[15]) && (z<=altl[7])) {
1573     /*   Turbopause */
1574     zh14=pdm[6][2];
1575     /*  Mixed density at Zlb */
1576     b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1577     /*  Mixed density at Alt */
1578     dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1579     zhm14=zhm28;
1580     /*  Net density at Alt */
1581     output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.);
1582     /*   Correction to specified mixing ratio at ground */
1583     rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14);
1584     hc14=pdm[6][5]*pdl[0][1];
1585     zc14=pdm[6][4]*pdl[0][0];
1586     output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14);
1587     /*   Chemistry correction */
1588     hcc14=pdm[6][7]*pdl[0][4];
1589     zcc14=pdm[6][6]*pdl[0][3];
1590     rc14=pdm[6][3]*pdl[0][5];
1591     /*  Net density corrected at Alt */
1592     output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14);
1593   }
1594
1595
1596         /**** Anomalous OXYGEN DENSITY ****/
1597
1598   g16h = flags->sw[21]*globe7(pd[8],input,flags);
1599   db16h = pdm[7][0]*exp(g16h)*pd[8][0];
1600   tho = pdm[7][9]*pdl[0][6];
1601   dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1602   zsht=pdm[7][5];
1603   zmho=pdm[7][4];
1604   zsho=scalh(zmho,16.0,tho);
1605   output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.));
1606
1607
1608   /* total mass density */
1609   output->d[5] = 1.66E-24*(4.0*output->d[0]+16.0*output->d[1]+28.0*output->d[2]+32.0*output->d[3]+40.0*output->d[4]+ output->d[6]+14.0*output->d[7]);
1610   db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14);
1611
1612
1613
1614   /* temperature */
1615   z = sqrt(input->alt*input->alt);
1616   ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1);
1617   if (flags->sw[0]) {
1618     for(i=0;i<9;i++)
1619       output->d[i]=output->d[i]*1.0E6;
1620     output->d[5]=output->d[5]/1000;
1621   }
1622 }
1623
1624
1625 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1626 //    The bitmasked value choices are as follows:
1627 //    unset: In this case (the default) JSBSim would only print
1628 //       out the normally expected messages, essentially echoing
1629 //       the config files as they are read. If the environment
1630 //       variable is not set, debug_lvl is set to 1 internally
1631 //    0: This requests JSBSim not to output any messages
1632 //       whatsoever.
1633 //    1: This value explicity requests the normal JSBSim
1634 //       startup messages
1635 //    2: This value asks for a message to be printed out when
1636 //       a class is instantiated
1637 //    4: When this value is set, a message is displayed when a
1638 //       FGModel object executes its Run() method
1639 //    8: When this value is set, various runtime state variables
1640 //       are printed out periodically
1641 //    16: When set various parameters are sanity checked and
1642 //       a message is printed out when they go out of bounds
1643
1644 void MSIS::Debug(int from)
1645 {
1646   if (debug_lvl <= 0) return;
1647
1648   if (debug_lvl & 1) { // Standard console startup message output
1649     if (from == 0) { // Constructor
1650     }
1651   }
1652   if (debug_lvl & 2 ) { // Instantiation/Destruction notification
1653     if (from == 0) cout << "Instantiated: MSIS" << endl;
1654     if (from == 1) cout << "Destroyed:    MSIS" << endl;
1655   }
1656   if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
1657   }
1658   if (debug_lvl & 8 ) { // Runtime state variables
1659   }
1660   if (debug_lvl & 16) { // Sanity checking
1661   }
1662   if (debug_lvl & 32) { // Turbulence
1663     if (first_pass && from == 2) {
1664       cout << "vTurbulenceNED(X), vTurbulenceNED(Y), vTurbulenceNED(Z), "
1665            << "vTurbulenceGrad(X), vTurbulenceGrad(Y), vTurbulenceGrad(Z), "
1666            << "vDirection(X), vDirection(Y), vDirection(Z), "
1667            << "Magnitude, "
1668            << "vTurbPQR(P), vTurbPQR(Q), vTurbPQR(R), " << endl;
1669     }
1670     if (from == 2) {
1671       cout << vTurbulenceNED << ", " << vTurbulenceGrad << ", " << vDirection << ", " << Magnitude << ", " << vTurbPQR << endl;
1672     }
1673   }
1674   if (debug_lvl & 64) {
1675     if (from == 0) { // Constructor
1676       cout << IdSrc << endl;
1677       cout << IdHdr << endl;
1678     }
1679   }
1680 }
1681
1682
1683
1684 } // namespace JSBSim
1685