1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 (incorporated into C++ JSBSim class heirarchy, see model authors below)
7 Purpose: Models the MSIS-00 atmosphere
9 ------------- Copyright (C) 2003 David P. Culp (davidculp2@comcast.net) ------
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
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
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.
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.
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.
34 --------------------------------------------------------------------------------
36 01/11/04 DPC Derived from FGAtmosphere
38 --------------------------------------------------------------------
39 --------- N R L M S I S E - 0 0 M O D E L 2 0 0 1 ----------
40 --------------------------------------------------------------------
42 This file is part of the NRLMSISE-00 C source code package - release
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
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.
56 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
62 #include <math.h> /* maths functions */
63 #include <stdlib.h> /* for malloc/free */
64 #include <stdio.h> /* for printf */
65 #include <iostream> // for cout, endl
69 static const char *IdSrc = "$Id$";
70 static const char *IdHdr = ID_MSIS;
72 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
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];
86 extern double ptm[10];
87 extern double pdm[8][10];
88 extern double pavgm[10];
90 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
92 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
95 MSIS::MSIS(FGFDMExec* fdmex) : FGAtmosphere(fdmex)
102 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
112 bool MSIS::InitModel(void)
114 FGModel::InitModel();
118 flags.switches[0] = 0;
119 for (i=1;i<24;i++) flags.switches[i] = 1;
121 for (i=0;i<7;i++) aph.a[i] = 100.0;
123 // set some common magnetic flux values
128 SLtemperature = intTemperature = 518.0;
129 SLpressure = intPressure = 2116.7;
130 SLdensity = intDensity = 0.002378;
131 SLsoundspeed = sqrt(2403.0832 * SLtemperature);
132 rSLtemperature = 1.0/intTemperature;
133 rSLpressure = 1.0/intPressure;
134 rSLdensity = 1.0/intDensity;
135 rSLsoundspeed = 1.0/SLsoundspeed;
136 temperature = &intTemperature;
137 pressure = &intPressure;
138 density = &intDensity;
145 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
149 if (FGModel::Run()) return true;
150 if (FDMExec->Holding()) return false;
152 //do temp, pressure, and density first
154 // get sea-level values
155 Calculate(Auxiliary->GetDayOfYear(),
156 Auxiliary->GetSecondsInDay(),
158 Propagate->GetLocation().GetLatitudeDeg(),
159 Propagate->GetLocation().GetLongitudeDeg());
160 SLtemperature = output.t[1] * 1.8;
161 SLdensity = output.d[5] * 1.940321;
162 SLpressure = 1716.488 * SLdensity * SLtemperature;
163 SLsoundspeed = sqrt(2403.0832 * SLtemperature);
164 rSLtemperature = 1.0/SLtemperature;
165 rSLpressure = 1.0/SLpressure;
166 rSLdensity = 1.0/SLdensity;
167 rSLsoundspeed = 1.0/SLsoundspeed;
169 // get at-altitude values
170 Calculate(Auxiliary->GetDayOfYear(),
171 Auxiliary->GetSecondsInDay(),
173 Propagate->GetLocation().GetLatitudeDeg(),
174 Propagate->GetLocation().GetLongitudeDeg());
175 intTemperature = output.t[1] * 1.8;
176 intDensity = output.d[5] * 1.940321;
177 intPressure = 1716.488 * intDensity * intTemperature;
178 soundspeed = sqrt(2403.0832 * intTemperature);
179 //cout << "T=" << intTemperature << " D=" << intDensity << " P=";
180 //cout << intPressure << " a=" << soundspeed << endl;
183 if (turbType != ttNone) {
185 vWindNED += vTurbulence;
188 if (vWindNED(1) != 0.0) psiw = atan2( vWindNED(2), vWindNED(1) );
190 if (psiw < 0) psiw += 2*M_PI;
197 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199 void MSIS::Calculate(int day, double sec, double alt, double lat, double lon)
204 input.alt = alt / 3281; //feet to kilometers
208 input.lst = (sec/3600) + (lon/15);
209 if (input.lst > 24.0) input.lst -= 24.0;
210 if (input.lst < 0.0) input.lst = 24 - input.lst;
212 gtd7d(&input, &flags, &output);
215 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
218 void MSIS::UseExternal(void){
219 // do nothing, external control not allowed
223 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
226 void MSIS::tselec(struct nrlmsise_flags *flags)
231 if (flags->switches[i]==1)
235 if (flags->switches[i]>0)
240 flags->sw[i]=flags->switches[i];
241 flags->swc[i]=flags->switches[i];
247 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
249 void MSIS::glatf(double lat, double *gv, double *reff)
251 double dgtr = 1.74533E-2;
253 c2 = cos(2.0*dgtr*lat);
254 *gv = 980.616 * (1.0 - 0.0026373 * c2);
255 *reff = 2.0 * (*gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5;
258 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
260 double MSIS::ccor(double alt, double r, double h1, double zh)
262 /* CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS
265 * H1 - transition scale length
266 * ZH - altitude of 1/2 R
280 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
282 double MSIS::ccor2(double alt, double r, double h1, double zh, double h2)
284 /* CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS
287 * H1 - transition scale length
288 * ZH - altitude of 1/2 R
289 * H2 - transition scale length #2 ?
294 e1 = (alt - zh) / h1;
295 e2 = (alt - zh) / h2;
296 if ((e1 > 70) || (e2 > 70))
298 if ((e1 < -70) && (e2 < -70))
302 ccor2v = r / (1.0 + 0.5 * (ex1 + ex2));
306 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308 double MSIS::scalh(double alt, double xm, double temp)
312 g = gsurf / (pow((1.0 + alt/re),2.0));
313 g = rgas * temp / (g * xm);
317 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
319 double MSIS::dnet (double dd, double dm, double zhm, double xmm, double xm)
321 /* TURBOPAUSE CORRECTION FOR MSIS MODELS
323 * DD - diffusive density
324 * DM - full mixed density
325 * ZHM - transition scale length
326 * XMM - full mixed molecular weight
327 * XM - species molecular weight
328 * DNET - combined density
333 if (!((dm>0) && (dd>0))) {
334 printf("dnet log error %e %e %e\n",dm,dd,xm);
335 if ((dd==0) && (dm==0))
342 ylog = a * log(dm/dd);
347 a = dd*pow((1.0 + exp(ylog)),(1.0/a));
351 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
353 void MSIS::splini (double *xa, double *ya, double *y2a, int n, double x, double *y)
355 /* INTEGRATE CUBIC SPLINE FUNCTION FROM XA(1) TO X
356 * XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
357 * Y2A: ARRAY OF SECOND DERIVATIVES
358 * N: SIZE OF ARRAYS XA,YA,Y2A
359 * X: ABSCISSA ENDPOINT FOR INTEGRATION
365 double xx, h, a, b, a2, b2;
366 while ((x>xa[klo]) && (khi<n)) {
374 h = xa[khi] - xa[klo];
375 a = (xa[khi] - xx)/h;
376 b = (xx - xa[klo])/h;
379 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;
386 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
388 void MSIS::splint (double *xa, double *ya, double *y2a, int n, double x, double *y)
390 /* CALCULATE CUBIC SPLINE INTERP VALUE
391 * ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL.
392 * XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
393 * Y2A: ARRAY OF SECOND DERIVATIVES
394 * N: SIZE OF ARRAYS XA,YA,Y2A
395 * X: ABSCISSA FOR INTERPOLATION
403 while ((khi-klo)>1) {
410 h = xa[khi] - xa[klo];
412 printf("bad XA input to splint");
415 yi = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0;
419 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
421 void MSIS::spline (double *x, double *y, int n, double yp1, double ypn, double *y2)
423 /* CALCULATE 2ND DERIVATIVES OF CUBIC SPLINE INTERP FUNCTION
424 * ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL
425 * X,Y: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
426 * N: SIZE OF ARRAYS X,Y
427 * YP1,YPN: SPECIFIED DERIVATIVES AT X[0] AND X[N-1]; VALUES
428 * >= 1E30 SIGNAL SIGNAL SECOND DERIVATIVE ZERO
429 * Y2: OUTPUT ARRAY OF SECOND DERIVATIVES
432 double sig, p, qn, un;
434 u=(double*)malloc(sizeof(double)*n);
436 printf("Out Of Memory in spline - ERROR");
444 u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
446 for (i=1;i<(n-1);i++) {
447 sig = (x[i]-x[i-1])/(x[i+1] - x[i-1]);
448 p = sig * y2[i-1] + 2.0;
449 y2[i] = (sig - 1.0) / p;
450 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;
457 un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])/(x[n-1] - x[n-2]));
459 y2[n-1] = (un - qn * u[n-2]) / (qn * y2[n-2] + 1.0);
461 y2[k] = y2[k] * y2[k+1] + u[k];
466 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
468 double MSIS::zeta(double zz, double zl)
470 return ((zz-zl)*(re+zl)/(re+zz));
473 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
475 double MSIS::densm(double alt, double d0, double xm, double *tz, int mn3,
476 double *zn3, double *tn3, double *tgn3, int mn2, double *zn2,
477 double *tn2, double *tgn2)
479 /* Calculate Temperature and Density Profiles for lower atmos. */
480 double xs[10], ys[10], y2out[10];
482 double z, z1, z2, t1, t2, zg, zgdif;
485 double expl, gamm, glb;
497 /* STRATOSPHERE/MESOSPHERE TEMPERATURE */
508 zgdif = zeta(z2, z1);
510 /* set up spline nodes */
512 xs[k]=zeta(zn2[k],z1)/zgdif;
515 yd1=-tgn2[0] / (t1*t1) * zgdif;
516 yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
518 /* calculate spline coefficients */
519 spline (xs, ys, mn, yd1, yd2, y2out);
521 splint (xs, ys, y2out, mn, x, &y);
523 /* temperature at altitude */
526 /* calaculate stratosphere / mesospehere density */
527 glb = gsurf / (pow((1.0 + z1/re),2.0));
528 gamm = xm * glb * zgdif / rgas;
530 /* Integrate temperature profile */
531 splini(xs, ys, y2out, mn, x, &yi);
536 /* Density at altitude */
537 densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
547 /* troposhere / stratosphere temperature */
557 /* set up spline nodes */
559 xs[k] = zeta(zn3[k],z1) / zgdif;
560 ys[k] = 1.0 / tn3[k];
562 yd1=-tgn3[0] / (t1*t1) * zgdif;
563 yd2=-tgn3[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
565 /* calculate spline coefficients */
566 spline (xs, ys, mn, yd1, yd2, y2out);
568 splint (xs, ys, y2out, mn, x, &y);
570 /* temperature at altitude */
573 /* calaculate tropospheric / stratosphere density */
574 glb = gsurf / (pow((1.0 + z1/re),2.0));
575 gamm = xm * glb * zgdif / rgas;
577 /* Integrate temperature profile */
578 splini(xs, ys, y2out, mn, x, &yi);
583 /* Density at altitude */
584 densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
592 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
594 double MSIS::densu(double alt, double dlb, double tinf, double tlb, double xm,
595 double alpha, double *tz, double zlb, double s2, int mn1,
596 double *zn1, double *tn1, double *tgn1)
598 /* Calculate Temperature and Density Profiles for MSIS models
599 * New lower thermo polynomial
601 double yd2, yd1, x=0.0, y=0.0;
603 double densu_temp=1.0;
604 double za, z, zg2, tt, ta=0.0;
605 double dta, z1=0.0, z2, t1=0.0, t2, zg, zgdif=0.0;
613 double xs[5], ys[5], y2out[5];
614 /* joining altitudes of Bates and spline */
621 /* geopotential altitude difference from ZLB */
624 /* Bates temperature */
625 tt = tinf - (tinf - tlb) * exp(-s2*zg2);
631 /* calculate temperature below ZA
632 * temperature gradient at ZA from Bates profile */
633 dta = (tinf - ta) * s2 * pow(((re+zlb)/(re+za)),2.0);
645 /* geopotental difference from z1 */
647 zgdif = zeta(z2, z1);
648 /* set up spline nodes */
650 xs[k] = zeta(zn1[k], z1) / zgdif;
651 ys[k] = 1.0 / tn1[k];
653 /* end node derivatives */
654 yd1 = -tgn1[0] / (t1*t1) * zgdif;
655 yd2 = -tgn1[1] / (t2*t2) * zgdif * pow(((re+z2)/(re+z1)),2.0);
656 /* calculate spline coefficients */
657 spline (xs, ys, mn, yd1, yd2, y2out);
659 splint (xs, ys, y2out, mn, x, &y);
660 /* temperature at altitude */
667 /* calculate density above za */
668 glb = gsurf / pow((1.0 + zlb/re),2.0);
669 gamma = xm * glb / (s2 * rgas * tinf);
670 expl = exp(-s2 * gamma * zg2);
676 /* density at altitude */
677 densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl;
682 /* calculate density below za */
683 glb = gsurf / pow((1.0 + z1/re),2.0);
684 gamm = xm * glb * zgdif / rgas;
686 /* integrate spline temperatures */
687 splini (xs, ys, y2out, mn, x, &yi);
694 /* density at altitude */
695 densu_temp = densu_temp * pow ((t1 / *tz),(1.0 + alpha)) * exp(-expl);
699 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
701 /* 3hr Magnetic activity functions */
703 double MSIS::g0(double a, double *p)
705 return (a - 4.0 + (p[25] - 1.0) * (a - 4.0 + (exp(-sqrt(p[24]*p[24]) *
706 (a - 4.0)) - 1.0) / sqrt(p[24]*p[24])));
709 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
712 double MSIS::sumex(double ex)
714 return (1.0 + (1.0 - pow(ex,19.0)) / (1.0 - ex) * pow(ex,0.5));
717 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
720 double MSIS::sg0(double ex, double *p, double *ap)
722 return (g0(ap[1],p) + (g0(ap[2],p)*ex + g0(ap[3],p)*ex*ex +
723 g0(ap[4],p)*pow(ex,3.0) + (g0(ap[5],p)*pow(ex,4.0) +
724 g0(ap[6],p)*pow(ex,12.0))*(1.0-pow(ex,8.0))/(1.0-ex)))/sumex(ex);
727 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
729 double MSIS::globe7(double *p, struct nrlmsise_input *input,
730 struct nrlmsise_flags *flags)
732 /* CALCULATE G(L) FUNCTION
733 * Upper Thermosphere Parameters */
740 double c, s, c2, c4, s2;
741 double sr = 7.2722E-5;
742 double dgtr = 1.74533E-2;
743 double dr = 1.72142E-2;
745 double cd32, cd18, cd14, cd39;
746 double p32, p18, p14, p39;
757 else if (flags->sw[9]<0)
759 xlong = input->g_long;
761 /* calculate legendre polynomials */
762 c = sin(input->g_lat * dgtr);
763 s = cos(input->g_lat * dgtr);
769 plg[0][2] = 0.5*(3.0*c2 -1.0);
770 plg[0][3] = 0.5*(5.0*c*c2-3.0*c);
771 plg[0][4] = (35.0*c4 - 30.0*c2 + 3.0)/8.0;
772 plg[0][5] = (63.0*c2*c2*c - 70.0*c2*c + 15.0*c)/8.0;
773 plg[0][6] = (11.0*c*plg[0][5] - 5.0*plg[0][4])/6.0;
774 /* plg[0][7] = (13.0*c*plg[0][6] - 6.0*plg[0][5])/7.0; */
777 plg[1][3] = 1.5*(5.0*c2-1.0)*s;
778 plg[1][4] = 2.5*(7.0*c2*c-3.0*c)*s;
779 plg[1][5] = 1.875*(21.0*c4 - 14.0*c2 +1.0)*s;
780 plg[1][6] = (11.0*c*plg[1][5]-6.0*plg[1][4])/5.0;
781 /* plg[1][7] = (13.0*c*plg[1][6]-7.0*plg[1][5])/6.0; */
782 /* plg[1][8] = (15.0*c*plg[1][7]-8.0*plg[1][6])/7.0; */
784 plg[2][3] = 15.0*s2*c;
785 plg[2][4] = 7.5*(7.0*c2 -1.0)*s2;
786 plg[2][5] = 3.0*c*plg[2][4]-2.0*plg[2][3];
787 plg[2][6] =(11.0*c*plg[2][5]-7.0*plg[2][4])/4.0;
788 plg[2][7] =(13.0*c*plg[2][6]-8.0*plg[2][5])/5.0;
789 plg[3][3] = 15.0*s2*s;
790 plg[3][4] = 105.0*s2*s*c;
791 plg[3][5] =(9.0*c*plg[3][4]-7.*plg[3][3])/2.0;
792 plg[3][6] =(11.0*c*plg[3][5]-8.*plg[3][4])/3.0;
794 if (!(((flags->sw[7]==0)&&(flags->sw[8]==0))&&(flags->sw[14]==0))) {
795 stloc = sin(hr*tloc);
796 ctloc = cos(hr*tloc);
797 s2tloc = sin(2.0*hr*tloc);
798 c2tloc = cos(2.0*hr*tloc);
799 s3tloc = sin(3.0*hr*tloc);
800 c3tloc = cos(3.0*hr*tloc);
803 cd32 = cos(dr*(input->doy-p[31]));
804 cd18 = cos(2.0*dr*(input->doy-p[17]));
805 cd14 = cos(dr*(input->doy-p[13]));
806 cd39 = cos(2.0*dr*(input->doy-p[38]));
813 df = input->f107 - input->f107A;
814 dfa = input->f107A - 150.0;
815 t[0] = p[19]*df*(1.0+p[59]*dfa) + p[20]*df*df + p[21]*dfa + p[29]*pow(dfa,2.0);
816 f1 = 1.0 + (p[47]*dfa +p[19]*df+p[20]*df*df)*flags->swc[1];
817 f2 = 1.0 + (p[49]*dfa+p[19]*df+p[20]*df*df)*flags->swc[1];
819 /* TIME INDEPENDENT */
820 t[1] = (p[1]*plg[0][2]+ p[2]*plg[0][4]+p[22]*plg[0][6]) +
821 (p[14]*plg[0][2])*dfa*flags->swc[1] +p[26]*plg[0][1];
823 /* SYMMETRICAL ANNUAL */
826 /* SYMMETRICAL SEMIANNUAL */
827 t[3] = (p[15]+p[16]*plg[0][2])*cd18;
829 /* ASYMMETRICAL ANNUAL */
830 t[4] = f1*(p[9]*plg[0][1]+p[10]*plg[0][3])*cd14;
832 /* ASYMMETRICAL SEMIANNUAL */
833 t[5] = p[37]*plg[0][1]*cd39;
838 t71 = (p[11]*plg[1][2])*cd14*flags->swc[5];
839 t72 = (p[12]*plg[1][2])*cd14*flags->swc[5];
840 t[6] = f2*((p[3]*plg[1][1] + p[4]*plg[1][3] + p[27]*plg[1][5] + t71) * \
841 ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + p[28]*plg[1][5] \
848 t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
849 t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
850 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);
855 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);
858 /* magnetic activity based on daily ap */
859 if (flags->sw[9]==-1) {
863 exp1 = exp(-10800.0*sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-sqrt(input->g_lat*input->g_lat))));
868 apt[0]=sg0(exp1,p,ap->a);
869 /* apt[1]=sg2(exp1,p,ap->a);
870 apt[2]=sg0(exp2,p,ap->a);
871 apt[3]=sg2(exp2,p,ap->a);
874 t[8] = apt[0]*(p[50]+p[96]*plg[0][2]+p[54]*plg[0][4]+ \
875 (p[125]*plg[0][1]+p[126]*plg[0][3]+p[127]*plg[0][5])*cd14*flags->swc[5]+ \
876 (p[128]*plg[1][1]+p[129]*plg[1][3]+p[130]*plg[1][5])*flags->swc[7]* \
877 cos(hr*(tloc-p[131])));
887 apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44);
889 t[8]=apdf*(p[32]+p[45]*plg[0][2]+p[34]*plg[0][4]+ \
890 (p[100]*plg[0][1]+p[101]*plg[0][3]+p[102]*plg[0][5])*cd14*flags->swc[5]+
891 (p[121]*plg[1][1]+p[122]*plg[1][3]+p[123]*plg[1][5])*flags->swc[7]*
892 cos(hr*(tloc-p[124])));
896 if ((flags->sw[10])&&(input->g_long>-1000.0)) {
900 t[10] = (1.0 + p[80]*dfa*flags->swc[1])* \
901 ((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
902 +p[103]*plg[1][1]+p[104]*plg[1][3]+p[105]*plg[1][5]\
903 +flags->swc[5]*(p[109]*plg[1][1]+p[110]*plg[1][3]+p[111]*plg[1][5])*cd14)* \
904 cos(dgtr*input->g_long) \
905 +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
906 +p[106]*plg[1][1]+p[107]*plg[1][3]+p[108]*plg[1][5]\
907 +flags->swc[5]*(p[112]*plg[1][1]+p[113]*plg[1][3]+p[114]*plg[1][5])*cd14)* \
908 sin(dgtr*input->g_long));
911 /* ut and mixed ut, longitude */
913 t[11]=(1.0+p[95]*plg[0][1])*(1.0+p[81]*dfa*flags->swc[1])*\
914 (1.0+p[119]*plg[0][1]*flags->swc[5]*cd14)*\
915 ((p[68]*plg[0][1]+p[69]*plg[0][3]+p[70]*plg[0][5])*\
916 cos(sr*(input->sec-p[71])));
917 t[11]+=flags->swc[11]*\
918 (p[76]*plg[2][3]+p[77]*plg[2][5]+p[78]*plg[2][7])*\
919 cos(sr*(input->sec-p[79])+2.0*dgtr*input->g_long)*(1.0+p[137]*dfa*flags->swc[1]);
922 /* ut, longitude magnetic activity */
924 if (flags->sw[9]==-1) {
926 t[12]=apt[0]*flags->swc[11]*(1.+p[132]*plg[0][1])*\
927 ((p[52]*plg[1][2]+p[98]*plg[1][4]+p[67]*plg[1][6])*\
928 cos(dgtr*(input->g_long-p[97])))\
929 +apt[0]*flags->swc[11]*flags->swc[5]*\
930 (p[133]*plg[1][1]+p[134]*plg[1][3]+p[135]*plg[1][5])*\
931 cd14*cos(dgtr*(input->g_long-p[136])) \
932 +apt[0]*flags->swc[12]* \
933 (p[55]*plg[0][1]+p[56]*plg[0][3]+p[57]*plg[0][5])*\
934 cos(sr*(input->sec-p[58]));
937 t[12] = apdf*flags->swc[11]*(1.0+p[120]*plg[0][1])*\
938 ((p[60]*plg[1][2]+p[61]*plg[1][4]+p[62]*plg[1][6])*\
939 cos(dgtr*(input->g_long-p[63])))\
940 +apdf*flags->swc[11]*flags->swc[5]* \
941 (p[115]*plg[1][1]+p[116]*plg[1][3]+p[117]*plg[1][5])* \
942 cd14*cos(dgtr*(input->g_long-p[118])) \
943 + apdf*flags->swc[12]* \
944 (p[83]*plg[0][1]+p[84]*plg[0][3]+p[85]*plg[0][5])* \
945 cos(sr*(input->sec-p[75]));
950 /* parms not used: 82, 89, 99, 139-149 */
953 tinf = tinf + fabs(flags->sw[i+1])*t[i];
957 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
959 double MSIS::glob7s(double *p, struct nrlmsise_input *input,
960 struct nrlmsise_flags *flags)
962 /* VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99
967 double cd32, cd18, cd14, cd39;
968 double p32, p18, p14, p39;
970 double dr=1.72142E-2;
971 double dgtr=1.74533E-2;
972 /* confirm parameter set */
976 printf("Wrong parameter set for glob7s\n");
981 cd32 = cos(dr*(input->doy-p[31]));
982 cd18 = cos(2.0*dr*(input->doy-p[17]));
983 cd14 = cos(dr*(input->doy-p[13]));
984 cd39 = cos(2.0*dr*(input->doy-p[38]));
993 /* time independent */
994 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];
996 /* SYMMETRICAL ANNUAL */
997 t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32;
999 /* SYMMETRICAL SEMIANNUAL */
1000 t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18;
1002 /* ASYMMETRICAL ANNUAL */
1003 t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14;
1005 /* ASYMMETRICAL SEMIANNUAL */
1006 t[5]=(p[37]*plg[0][1])*cd39;
1011 t71 = p[11]*plg[1][2]*cd14*flags->swc[5];
1012 t72 = p[12]*plg[1][2]*cd14*flags->swc[5];
1013 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) ;
1019 t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
1020 t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
1021 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);
1025 if (flags->sw[14]) {
1026 t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
1029 /* MAGNETIC ACTIVITY */
1031 if (flags->sw[9]==1)
1032 t[8] = apdf * (p[32] + p[45] * plg[0][2] * flags->swc[2]);
1033 if (flags->sw[9]==-1)
1034 t[8]=(p[50]*apt[0] + p[96]*plg[0][2] * apt[0]*flags->swc[2]);
1038 if (!((flags->sw[10]==0) || (flags->sw[11]==0) || (input->g_long<=-1000.0))) {
1039 t[10] = (1.0 + plg[0][1]*(p[80]*flags->swc[5]*cos(dr*(input->doy-p[81]))\
1040 +p[85]*flags->swc[6]*cos(2.0*dr*(input->doy-p[86])))\
1041 +p[83]*flags->swc[3]*cos(dr*(input->doy-p[84]))\
1042 +p[87]*flags->swc[4]*cos(2.0*dr*(input->doy-p[88])))\
1043 *((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
1044 +p[74]*plg[1][1]+p[75]*plg[1][3]+p[76]*plg[1][5]\
1045 )*cos(dgtr*input->g_long)\
1046 +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
1047 +p[77]*plg[1][1]+p[78]*plg[1][3]+p[79]*plg[1][5]\
1048 )*sin(dgtr*input->g_long));
1052 tt+=fabs(flags->sw[i+1])*t[i];
1056 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1058 void MSIS::gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1059 struct nrlmsise_output *output)
1064 double zn3[5]={32.5,20.0,15.0,10.0,0.0};
1066 double zn2[4]={72.5,55.0,45.0,32.5};
1075 struct nrlmsise_output soutput;
1080 /* Latitude variation of gravity (none for sw[2]=0) */
1082 if (flags->sw[2]==0)
1084 glatf(xlat, &gsurf, &re);
1088 /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */
1089 if (input->alt>zn2[0])
1096 gts7(input, flags, &soutput);
1099 if (flags->sw[0]) /* metric adjustment */
1103 output->t[0]=soutput.t[0];
1104 output->t[1]=soutput.t[1];
1105 if (input->alt>=zn2[0]) {
1107 output->d[i]=soutput.d[i];
1111 /* LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])
1112 * Temperature at nodes and gradients at end nodes
1113 * Inverse temperature a linear function of spherical harmonics
1115 meso_tgn2[0]=meso_tgn1[1];
1116 meso_tn2[0]=meso_tn1[4];
1117 meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags));
1118 meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags));
1119 meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags));
1120 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));
1121 meso_tn3[0]=meso_tn2[3];
1123 if (input->alt<zn3[0]) {
1124 /* LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])
1125 * Temperature at nodes and gradients at end nodes
1126 * Inverse temperature a linear function of spherical harmonics
1128 meso_tgn3[0]=meso_tgn2[1];
1129 meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags->sw[22]*glob7s(pma[3], input, flags));
1130 meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags));
1131 meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags));
1132 meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags));
1133 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));
1136 /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */
1139 if (input->alt>zmix)
1140 dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix);
1143 /**** N2 density ****/
1144 dmr=soutput.d[2] / dm28m - 1.0;
1145 output->d[2]=densm(input->alt,dm28m,xmm, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
1146 output->d[2]=output->d[2] * (1.0 + dmr*dmc);
1148 /**** HE density ****/
1149 dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0;
1150 output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc);
1152 /**** O density ****/
1156 /**** O2 density ****/
1157 dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0;
1158 output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc);
1160 /**** AR density ***/
1161 dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0;
1162 output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc);
1164 /**** Hydrogen density ****/
1167 /**** Atomic nitrogen density ****/
1170 /**** Total mass density */
1171 output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] +
1172 28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4]
1173 + output->d[6] + 14.0 * output->d[7]);
1176 output->d[5]=output->d[5]/1000;
1178 /**** temperature at altitude ****/
1179 dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3,
1180 mn2, zn2, meso_tn2, meso_tgn2);
1185 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1187 void MSIS::gtd7d(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1188 struct nrlmsise_output *output)
1190 gtd7(input, flags, output);
1191 output->d[5] = 1.66E-24 * (4.0 * output->d[0] + 16.0 * output->d[1] +
1192 28.0 * output->d[2] + 32.0 * output->d[3] + 40.0 * output->d[4]
1193 + output->d[6] + 14.0 * output->d[7] + 16.0 * output->d[8]);
1196 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1198 void MSIS::ghp7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1199 struct nrlmsise_output *output, double press)
1201 double bm = 1.3806E-19;
1202 double rgas = 831.4;
1203 double test = 0.00043;
1210 double xn, xm, diff;
1216 zi = 18.06 * (3.00 - pl);
1217 else if ((pl>0.075) && (pl<=2.5))
1218 zi = 14.98 * (3.08 - pl);
1219 else if ((pl>-1) && (pl<=0.075))
1220 zi = 17.80 * (2.72 - pl);
1221 else if ((pl>-2) && (pl<=-1))
1222 zi = 14.28 * (3.64 - pl);
1223 else if ((pl>-4) && (pl<=-2))
1224 zi = 12.72 * (4.32 -pl);
1226 zi = 25.3 * (0.11 - pl);
1227 cl = input->g_lat/90.0;
1230 cd = (1.0 - (double) input->doy) / 91.25;
1232 cd = ((double) input->doy) / 91.25 - 3.0;
1234 if ((pl > -1.11) && (pl<=-0.23))
1237 ca = (2.79 - pl) / (2.79 + 0.23);
1238 if ((pl <= -1.11) && (pl>-3))
1239 ca = (-2.93 - pl)/(-2.93 + 1.11);
1240 z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl;
1242 z = 22.0 * pow((pl + 4.0),2.0) + 110.0;
1244 /* iteration loop */
1249 gtd7(input, flags, output);
1251 xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7];
1252 p = bm * xn * output->t[1];
1255 diff = pl - log10(p);
1256 if (sqrt(diff*diff)<test)
1259 printf("ERROR: ghp7 not converging for press %e, diff %e",press,diff);
1262 xm = output->d[5] / xn / 1.66E-24;
1265 g = gsurf / (pow((1.0 + z/re),2.0));
1266 sh = rgas * output->t[1] / (xm * g);
1268 /* new altitude estimate using scale height */
1270 z = z - sh * diff * 2.302;
1276 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1278 void MSIS::gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags,
1279 struct nrlmsise_output *output)
1281 /* Thermospheric portion of NRLMSISE-00
1282 * See GTD7 for more extensive comments
1288 double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
1293 double s, z0, t0, tr12;
1294 double db01, db04, db14, db16, db28, db32, db40, db48;
1295 double zh28, zh04, zh16, zh32, zh40, zh01, zh14;
1296 double zhm28, zhm04, zhm16, zhm32, zhm40, zhm01, zhm14;
1298 double b28, b04, b16, b32, b40, b01, b14;
1300 double g28, g4, g16, g32, g40, g1, g14;
1302 double zc04, zc16, zc32, zc40, zc01, zc14;
1303 double hc04, hc16, hc32, hc40, hc01, hc14;
1304 double hcc16, hcc32, hcc01, hcc14;
1305 double zcc16, zcc32, zcc01, zcc14;
1306 double rc16, rc32, rc01, rc14;
1308 double g16h, db16h, tho, zsht, zmho, zsho;
1309 double dgtr=1.74533E-2;
1310 double dr=1.72142E-2;
1311 double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
1312 double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
1314 double hc216, hcc232;
1320 /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */
1321 if (input->alt>zn1[0])
1322 tinf = ptm[0]*pt[0] * \
1323 (1.0+flags->sw[16]*globe7(pt,input,flags));
1325 tinf = ptm[0]*pt[0];
1328 /* GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */
1329 if (input->alt>zn1[4])
1330 g0 = ptm[3]*ps[0] * \
1331 (1.0+flags->sw[19]*globe7(ps,input,flags));
1334 tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0];
1335 s = g0 / (tinf - tlb);
1337 /* Lower thermosphere temp variations not significant for
1338 * density above 300 km */
1339 if (input->alt<300.0) {
1340 meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags));
1341 meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags));
1342 meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags));
1343 meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags));
1344 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));
1346 meso_tn1[1]=ptm[6]*ptl[0][0];
1347 meso_tn1[2]=ptm[2]*ptl[1][0];
1348 meso_tn1[3]=ptm[7]*ptl[2][0];
1349 meso_tn1[4]=ptm[4]*ptl[3][0];
1350 meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0));
1357 /* N2 variation factor at Zlb */
1358 g28=flags->sw[21]*globe7(pd[2], input, flags);
1360 /* VARIATION OF TURBOPAUSE HEIGHT */
1361 zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13])));
1367 /**** N2 DENSITY ****/
1369 /* Diffusive density at Zlb */
1370 db28 = pdm[2][0]*exp(g28)*pd[2][0];
1371 /* Diffusive density at Alt */
1372 output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1376 zhm28=pdm[2][3]*pdl[1][5];
1378 /* Mixed density at Zlb */
1379 b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1380 if ((flags->sw[15])&&(z<=altl[2])) {
1381 /* Mixed density at Alt */
1382 dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1383 /* Net density at Alt */
1384 output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0);
1388 /**** HE DENSITY ****/
1390 /* Density variation factor at Zlb */
1391 g4 = flags->sw[21]*globe7(pd[0], input, flags);
1392 /* Diffusive density at Zlb */
1393 db04 = pdm[0][0]*exp(g4)*pd[0][0];
1394 /* Diffusive density at Alt */
1395 output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1397 if ((flags->sw[15]) && (z<altl[0])) {
1400 /* Mixed density at Zlb */
1401 b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1402 /* Mixed density at Alt */
1403 dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1405 /* Net density at Alt */
1406 output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.);
1407 /* Correction to specified mixing ratio at ground */
1408 rl=log(b28*pdm[0][1]/b04);
1409 zc04=pdm[0][4]*pdl[1][0];
1410 hc04=pdm[0][5]*pdl[1][1];
1411 /* Net density corrected at Alt */
1412 output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
1416 /**** O DENSITY ****/
1418 /* Density variation factor at Zlb */
1419 g16= flags->sw[21]*globe7(pd[1],input,flags);
1420 /* Diffusive density at Zlb */
1421 db16 = pdm[1][0]*exp(g16)*pd[1][0];
1422 /* Diffusive density at Alt */
1423 output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1425 if ((flags->sw[15]) && (z<=altl[1])) {
1428 /* Mixed density at Zlb */
1429 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);
1430 /* Mixed density at Alt */
1431 dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1433 /* Net density at Alt */
1434 output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.);
1435 rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0));
1436 hc16=pdm[1][5]*pdl[1][3];
1437 zc16=pdm[1][4]*pdl[1][2];
1438 hc216=pdm[1][5]*pdl[1][4];
1439 output->d[1]=output->d[1]*ccor2(z,rl,hc16,zc16,hc216);
1440 /* Chemistry correction */
1441 hcc16=pdm[1][7]*pdl[1][13];
1442 zcc16=pdm[1][6]*pdl[1][12];
1443 rc16=pdm[1][3]*pdl[1][14];
1444 /* Net density corrected at Alt */
1445 output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16);
1449 /**** O2 DENSITY ****/
1451 /* Density variation factor at Zlb */
1452 g32= flags->sw[21]*globe7(pd[4], input, flags);
1453 /* Diffusive density at Zlb */
1454 db32 = pdm[3][0]*exp(g32)*pd[4][0];
1455 /* Diffusive density at Alt */
1456 output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1458 if (flags->sw[15]) {
1462 /* Mixed density at Zlb */
1463 b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1464 /* Mixed density at Alt */
1465 dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1467 /* Net density at Alt */
1468 output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.);
1469 /* Correction to specified mixing ratio at ground */
1470 rl=log(b28*pdm[3][1]/b32);
1471 hc32=pdm[3][5]*pdl[1][7];
1472 zc32=pdm[3][4]*pdl[1][6];
1473 output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32);
1475 /* Correction for general departure from diffusive equilibrium above Zlb */
1476 hcc32=pdm[3][7]*pdl[1][22];
1477 hcc232=pdm[3][7]*pdl[0][22];
1478 zcc32=pdm[3][6]*pdl[1][21];
1479 rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.));
1480 /* Net density corrected at Alt */
1481 output->d[3]=output->d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
1485 /**** AR DENSITY ****/
1487 /* Density variation factor at Zlb */
1488 g40= flags->sw[20]*globe7(pd[5],input,flags);
1489 /* Diffusive density at Zlb */
1490 db40 = pdm[4][0]*exp(g40)*pd[5][0];
1491 /* Diffusive density at Alt */
1492 output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1494 if ((flags->sw[15]) && (z<=altl[4])) {
1497 /* Mixed density at Zlb */
1498 b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1499 /* Mixed density at Alt */
1500 dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1502 /* Net density at Alt */
1503 output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.);
1504 /* Correction to specified mixing ratio at ground */
1505 rl=log(b28*pdm[4][1]/b40);
1506 hc40=pdm[4][5]*pdl[1][9];
1507 zc40=pdm[4][4]*pdl[1][8];
1508 /* Net density corrected at Alt */
1509 output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40);
1513 /**** HYDROGEN DENSITY ****/
1515 /* Density variation factor at Zlb */
1516 g1 = flags->sw[21]*globe7(pd[6], input, flags);
1517 /* Diffusive density at Zlb */
1518 db01 = pdm[5][0]*exp(g1)*pd[6][0];
1519 /* Diffusive density at Alt */
1520 output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1522 if ((flags->sw[15]) && (z<=altl[6])) {
1525 /* Mixed density at Zlb */
1526 b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1527 /* Mixed density at Alt */
1528 dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1530 /* Net density at Alt */
1531 output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.);
1532 /* Correction to specified mixing ratio at ground */
1533 rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01);
1534 hc01=pdm[5][5]*pdl[1][11];
1535 zc01=pdm[5][4]*pdl[1][10];
1536 output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01);
1537 /* Chemistry correction */
1538 hcc01=pdm[5][7]*pdl[1][19];
1539 zcc01=pdm[5][6]*pdl[1][18];
1540 rc01=pdm[5][3]*pdl[1][20];
1541 /* Net density corrected at Alt */
1542 output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01);
1546 /**** ATOMIC NITROGEN DENSITY ****/
1548 /* Density variation factor at Zlb */
1549 g14 = flags->sw[21]*globe7(pd[7],input,flags);
1550 /* Diffusive density at Zlb */
1551 db14 = pdm[6][0]*exp(g14)*pd[7][0];
1552 /* Diffusive density at Alt */
1553 output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1555 if ((flags->sw[15]) && (z<=altl[7])) {
1558 /* Mixed density at Zlb */
1559 b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1560 /* Mixed density at Alt */
1561 dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1563 /* Net density at Alt */
1564 output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.);
1565 /* Correction to specified mixing ratio at ground */
1566 rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14);
1567 hc14=pdm[6][5]*pdl[0][1];
1568 zc14=pdm[6][4]*pdl[0][0];
1569 output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14);
1570 /* Chemistry correction */
1571 hcc14=pdm[6][7]*pdl[0][4];
1572 zcc14=pdm[6][6]*pdl[0][3];
1573 rc14=pdm[6][3]*pdl[0][5];
1574 /* Net density corrected at Alt */
1575 output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14);
1579 /**** Anomalous OXYGEN DENSITY ****/
1581 g16h = flags->sw[21]*globe7(pd[8],input,flags);
1582 db16h = pdm[7][0]*exp(g16h)*pd[8][0];
1583 tho = pdm[7][9]*pdl[0][6];
1584 dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1587 zsho=scalh(zmho,16.0,tho);
1588 output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.));
1591 /* total mass density */
1592 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]);
1593 db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14);
1598 z = sqrt(input->alt*input->alt);
1599 ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1);
1602 output->d[i]=output->d[i]*1.0E6;
1603 output->d[5]=output->d[5]/1000;
1608 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1609 // The bitmasked value choices are as follows:
1610 // unset: In this case (the default) JSBSim would only print
1611 // out the normally expected messages, essentially echoing
1612 // the config files as they are read. If the environment
1613 // variable is not set, debug_lvl is set to 1 internally
1614 // 0: This requests JSBSim not to output any messages
1616 // 1: This value explicity requests the normal JSBSim
1618 // 2: This value asks for a message to be printed out when
1619 // a class is instantiated
1620 // 4: When this value is set, a message is displayed when a
1621 // FGModel object executes its Run() method
1622 // 8: When this value is set, various runtime state variables
1623 // are printed out periodically
1624 // 16: When set various parameters are sanity checked and
1625 // a message is printed out when they go out of bounds
1627 void MSIS::Debug(int from)
1629 if (debug_lvl <= 0) return;
1631 if (debug_lvl & 1) { // Standard console startup message output
1632 if (from == 0) { // Constructor
1635 if (debug_lvl & 2 ) { // Instantiation/Destruction notification
1636 if (from == 0) cout << "Instantiated: MSIS" << endl;
1637 if (from == 1) cout << "Destroyed: MSIS" << endl;
1639 if (debug_lvl & 4 ) { // Run() method entry print for FGModel-derived objects
1641 if (debug_lvl & 8 ) { // Runtime state variables
1643 if (debug_lvl & 16) { // Sanity checking
1645 if (debug_lvl & 32) { // Turbulence
1646 if (first_pass && from == 2) {
1647 cout << "vTurbulence(X), vTurbulence(Y), vTurbulence(Z), "
1648 << "vTurbulenceGrad(X), vTurbulenceGrad(Y), vTurbulenceGrad(Z), "
1649 << "vDirection(X), vDirection(Y), vDirection(Z), "
1651 << "vTurbPQR(P), vTurbPQR(Q), vTurbPQR(R), " << endl;
1654 cout << vTurbulence << ", " << vTurbulenceGrad << ", " << vDirection << ", " << Magnitude << ", " << vTurbPQR << endl;
1657 if (debug_lvl & 64) {
1658 if (from == 0) { // Constructor
1659 cout << IdSrc << endl;
1660 cout << IdHdr << endl;
1667 } // namespace JSBSim