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