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