61 #include "models/FGAuxiliary.h" 69 IDENT(IdSrc,
"$Id: FGMSIS.cpp,v 1.24 2014/01/13 10:46:07 ehofman Exp $");
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];
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;
102 dm04 = dm16 = dm28 = dm32 = dm40 = dm01 = dm14 = dfa = 0.0;
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;
123 bool MSIS::InitModel(
void)
127 flags.switches[0] = 0;
131 flags.switches[i] = 1;
136 for (i=0;i<7;i++) aph.a[i] = 100.0;
162 if (Holding)
return false;
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;
203 void MSIS::Calculate(
int day,
double sec,
double alt,
double lat,
double lon)
208 input.alt = alt / 3281;
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;
216 gtd7d(&input, &flags, &output);
235 if (flags->switches[i]==1)
239 if (flags->switches[i]>0)
244 flags->sw[i]=flags->switches[i];
245 flags->swc[i]=flags->switches[i];
253 void MSIS::glatf(
double lat,
double *gv,
double *reff)
255 double dgtr = 1.74533E-2;
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;
264 double MSIS::ccor(
double alt,
double r,
double h1,
double zh)
286 double MSIS::ccor2(
double alt,
double r,
double h1,
double zh,
double h2)
298 e1 = (alt - zh) / h1;
299 e2 = (alt - zh) / h2;
300 if ((e1 > 70) || (e2 > 70))
302 if ((e1 < -70) && (e2 < -70))
306 ccor2v = r / (1.0 + 0.5 * (ex1 + ex2));
312 double MSIS::scalh(
double alt,
double xm,
double temp)
316 g = gsurf / (pow((1.0 + alt/re),2.0));
317 g = rgas * temp / (g * xm);
323 double MSIS::dnet (
double dd,
double dm,
double zhm,
double xmm,
double xm)
337 if (!((dm>0) && (dd>0))) {
338 cerr <<
"dnet log error " << dm <<
' ' << dd <<
' ' << xm <<
' ' << endl;
339 if ((dd==0) && (dm==0))
346 ylog = a * log(dm/dd);
351 a = dd*pow((1.0 + exp(ylog)),(1.0/a));
357 void MSIS::splini (
double *xa,
double *ya,
double *y2a,
int n,
double x,
double *y)
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)) {
378 h = xa[khi] - xa[klo];
379 a = (xa[khi] - xx)/h;
380 b = (xx - xa[klo])/h;
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;
392 void MSIS::splint (
double *xa,
double *ya,
double *y2a,
int n,
double x,
double *y)
407 while ((khi-klo)>1) {
414 h = xa[khi] - xa[klo];
416 cerr <<
"bad XA input to splint" << endl;
419 yi = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0;
425 void MSIS::spline (
double *x,
double *y,
int n,
double yp1,
double ypn,
double *y2)
436 double sig, p, qn, un;
440 cerr <<
"Out Of Memory in spline - ERROR" << endl;
448 u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
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;
461 un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])/(x[n-1] - x[n-2]));
463 y2[n-1] = (un - qn * u[n-2]) / (qn * y2[n-2] + 1.0);
465 y2[k] = y2[k] * y2[k+1] + u[k];
472 double MSIS::zeta(
double zz,
double zl)
474 return ((zz-zl)*(re+zl)/(re+zz));
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)
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};
488 double z=0, z1=0, z2=0, t1=0, t2=0, zg=0, zgdif=0;
490 double x=0, y=0, yi=0;
491 double expl=0, gamm=0, glb=0;
514 zgdif = zeta(z2, z1);
518 xs[k]=zeta(zn2[k],z1)/zgdif;
521 yd1=-tgn2[0] / (t1*t1) * zgdif;
522 yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
525 spline (xs, ys, mn, yd1, yd2, y2out);
527 splint (xs, ys, y2out, mn, x, &y);
533 glb = gsurf / (pow((1.0 + z1/re),2.0));
534 gamm = xm * glb * zgdif / rgas;
537 splini(xs, ys, y2out, mn, x, &yi);
543 densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
565 xs[k] = zeta(zn3[k],z1) / zgdif;
566 ys[k] = 1.0 / tn3[k];
568 yd1=-tgn3[0] / (t1*t1) * zgdif;
569 yd2=-tgn3[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0));
572 spline (xs, ys, mn, yd1, yd2, y2out);
574 splint (xs, ys, y2out, mn, x, &y);
580 glb = gsurf / (pow((1.0 + z1/re),2.0));
581 gamm = xm * glb * zgdif / rgas;
584 splini(xs, ys, y2out, mn, x, &yi);
590 densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl);
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)
607 double yd2=0.0, yd1=0.0, x=0.0, y=0.0;
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;
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};
631 tt = tinf - (tinf - tlb) * exp(-s2*zg2);
639 dta = (tinf - ta) * s2 * pow(((re+zlb)/(re+za)),2.0);
653 zgdif = zeta(z2, z1);
656 xs[k] = zeta(zn1[k], z1) / zgdif;
657 ys[k] = 1.0 / tn1[k];
660 yd1 = -tgn1[0] / (t1*t1) * zgdif;
661 yd2 = -tgn1[1] / (t2*t2) * zgdif * pow(((re+z2)/(re+z1)),2.0);
663 spline (xs, ys, mn, yd1, yd2, y2out);
665 splint (xs, ys, y2out, mn, x, &y);
674 glb = gsurf / pow((1.0 + zlb/re),2.0);
675 gamma = xm * glb / (s2 * rgas * tinf);
676 expl = exp(-s2 * gamma * zg2);
683 densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl;
689 glb = gsurf / pow((1.0 + z1/re),2.0);
690 gamm = xm * glb * zgdif / rgas;
693 splini (xs, ys, y2out, mn, x, &yi);
701 densu_temp = densu_temp * pow ((t1 / *tz),(1.0 + alpha)) * exp(-expl);
709 double MSIS::g0(
double a,
double *p)
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])));
718 double MSIS::sumex(
double ex)
720 return (1.0 + (1.0 - pow(ex,19.0)) / (1.0 - ex) * pow(ex,0.5));
726 double MSIS::sg0(
double ex,
double *p,
double *ap)
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);
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;
749 double cd32, cd18, cd14, cd39;
760 c = sin(input->g_lat * dgtr);
761 s = cos(input->g_lat * dgtr);
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;
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;
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;
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);
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]));
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];
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];
821 t[3] = (p[15]+p[16]*plg[0][2])*cd18;
824 t[4] = f1*(p[9]*plg[0][1]+p[10]*plg[0][3])*cd14;
827 t[5] = p[37]*plg[0][1]*cd39;
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] \
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);
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);
853 if (flags->sw[9]==-1) {
857 exp1 = exp(-10800.0*sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-sqrt(input->g_lat*input->g_lat))));
862 apt[0]=sg0(exp1,p,ap->a);
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])));
881 apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44);
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])));
890 if ((flags->sw[10])&&(input->g_long>-1000.0)) {
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));
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]);
918 if (flags->sw[9]==-1) {
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]));
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]));
947 tinf = tinf + fabs(flags->sw[i+1])*t[i];
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,};
961 double cd32=0.0, cd18=0.0, cd14=0.0, cd39=0.0;
963 double dr=1.72142E-2;
964 double dgtr=1.74533E-2;
969 cerr <<
"Wrong parameter set for glob7s" << endl;
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]));
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];
986 t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32;
989 t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18;
992 t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14;
995 t[5]=(p[37]*plg[0][1])*cd39;
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) ;
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);
1014 if (flags->sw[14]) {
1015 t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
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]);
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));
1041 tt+=fabs(flags->sw[i+1])*t[i];
1053 double zn3[5]={32.5,20.0,15.0,10.0,0.0};
1055 double zn2[4]={72.5,55.0,45.0,32.5};
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;
1074 if (flags->sw[2]==0)
1076 glatf(xlat, &gsurf, &re);
1081 if (input->alt>zn2[0])
1088 gts7(input, flags, &soutput);
1095 output->t[0]=soutput.t[0];
1096 output->t[1]=soutput.t[1];
1097 if (input->alt>=zn2[0]) {
1099 output->d[i]=soutput.d[i];
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];
1115 if (input->alt<zn3[0]) {
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));
1131 if (input->alt>zmix)
1132 dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix);
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);
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);
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);
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);
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]);
1168 output->d[5]=output->d[5]/1000;
1171 dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3,
1172 mn2, zn2, meso_tn2, meso_tgn2);
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]);
1187 output->d[5]=output->d[5]/1000;
1195 double bm = 1.3806E-19;
1196 double rgas = 831.4;
1197 double test = 0.00043;
1204 double xn, xm, diff;
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);
1220 zi = 25.3 * (0.11 - pl);
1221 cl = input->g_lat/90.0;
1224 cd = (1.0 - (double) input->doy) / 91.25;
1226 cd = ((double) input->doy) / 91.25 - 3.0;
1228 if ((pl > -1.11) && (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;
1236 z = 22.0 * pow((pl + 4.0),2.0) + 110.0;
1243 gtd7(input, flags, output);
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];
1249 diff = pl - log10(p);
1250 if (sqrt(diff*diff)<test)
1253 cerr <<
"ERROR: ghp7 not converging for press " << press <<
", diff " << diff << endl;
1256 xm = output->d[5] / xn / 1.66E-24;
1259 g = gsurf / (pow((1.0 + z/re),2.0));
1260 sh = rgas * output->t[1] / (xm * g);
1264 z = z - sh * diff * 2.302;
1282 double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
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;
1292 double b28=0.0, b04=0.0, b16=0.0, b32=0.0, b40=0.0, b01=0.0, b14=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;
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};
1308 double hc216=0.0, hcc232=0.0;
1315 if (input->alt>zn1[0])
1316 tinf = ptm[0]*pt[0] * \
1317 (1.0+flags->sw[16]*globe7(pt,input,flags));
1319 tinf = ptm[0]*pt[0];
1323 if (input->alt>zn1[4])
1324 g0 = ptm[3]*ps[0] * \
1325 (1.0+flags->sw[19]*globe7(ps,input,flags));
1328 tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0];
1329 s = g0 / (tinf - tlb);
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));
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));
1348 g28=flags->sw[21]*globe7(pd[2], input, flags);
1351 zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13])));
1360 db28 = pdm[2][0]*exp(g28)*pd[2][0];
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);
1366 zhm28=pdm[2][3]*pdl[1][5];
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])) {
1372 dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1374 output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0);
1381 g4 = flags->sw[21]*globe7(pd[0], input, flags);
1383 db04 = pdm[0][0]*exp(g4)*pd[0][0];
1385 output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1387 if ((flags->sw[15]) && (z<altl[0])) {
1391 b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1393 dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1396 output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.);
1398 rl=log(b28*pdm[0][1]/b04);
1399 zc04=pdm[0][4]*pdl[1][0];
1400 hc04=pdm[0][5]*pdl[1][1];
1402 output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
1409 g16= flags->sw[21]*globe7(pd[1],input,flags);
1411 db16 = pdm[1][0]*exp(g16)*pd[1][0];
1413 output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1415 if ((flags->sw[15]) && (z<=altl[1])) {
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);
1421 dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
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);
1431 hcc16=pdm[1][7]*pdl[1][13];
1432 zcc16=pdm[1][6]*pdl[1][12];
1433 rc16=pdm[1][3]*pdl[1][14];
1435 output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16);
1442 g32= flags->sw[21]*globe7(pd[4], input, flags);
1444 db32 = pdm[3][0]*exp(g32)*pd[4][0];
1446 output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
1448 if (flags->sw[15]) {
1453 b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1455 dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1458 output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.);
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);
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.));
1471 output->d[3]=output->d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
1478 g40= flags->sw[21]*globe7(pd[5],input,flags);
1480 db40 = pdm[4][0]*exp(g40)*pd[5][0];
1482 output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1484 if ((flags->sw[15]) && (z<=altl[4])) {
1488 b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1490 dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1493 output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.);
1495 rl=log(b28*pdm[4][1]/b40);
1496 hc40=pdm[4][5]*pdl[1][9];
1497 zc40=pdm[4][4]*pdl[1][8];
1499 output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40);
1506 g1 = flags->sw[21]*globe7(pd[6], input, flags);
1508 db01 = pdm[5][0]*exp(g1)*pd[6][0];
1510 output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1512 if ((flags->sw[15]) && (z<=altl[6])) {
1516 b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1518 dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1521 output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.);
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);
1528 hcc01=pdm[5][7]*pdl[1][19];
1529 zcc01=pdm[5][6]*pdl[1][18];
1530 rc01=pdm[5][3]*pdl[1][20];
1532 output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01);
1539 g14 = flags->sw[21]*globe7(pd[7],input,flags);
1541 db14 = pdm[6][0]*exp(g14)*pd[7][0];
1543 output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1545 if ((flags->sw[15]) && (z<=altl[7])) {
1549 b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1551 dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
1554 output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.);
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);
1561 hcc14=pdm[6][7]*pdl[0][4];
1562 zcc14=pdm[6][6]*pdl[0][3];
1563 rc14=pdm[6][3]*pdl[0][5];
1565 output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14);
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);
1577 zsho=scalh(zmho,16.0,tho);
1578 output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.));
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]);
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);
1591 output->d[i]=output->d[i]*1.0E6;
1592 output->d[5]=output->d[5]/1000;
1616 void MSIS::Debug(
int from)
1618 if (debug_lvl <= 0)
return;
1620 if (debug_lvl & 1) {
1624 if (debug_lvl & 2 ) {
1625 if (from == 0) cout <<
"Instantiated: MSIS" << endl;
1626 if (from == 1) cout <<
"Destroyed: MSIS" << endl;
1628 if (debug_lvl & 4 ) {
1630 if (debug_lvl & 8 ) {
1632 if (debug_lvl & 16) {
1634 if (debug_lvl & 32) {
1636 if (debug_lvl & 64) {
1638 cout << IdSrc << endl;
1639 cout << IdHdr << endl;
FGAuxiliary * GetAuxiliary(void)
Returns the FGAuxiliary pointer.
double GetAltitudeASL(void) const
Returns the current altitude above sea level.
virtual bool Run(bool Holding)
Runs the model; called by the Executive.
void UseExternal(void)
Does nothing. External control is not allowed.
bool Run(bool Holding)
Runs the MSIS-00 atmosphere model; called by the Executive Can pass in a value indicating if the exec...
Models an empty, abstract base atmosphere class.
Models the MSIS-00 atmosphere.
double GetLongitudeDeg() const
Get the longitude.
double GetLatitudeDeg() const
Get the latitude.
Encapsulates the JSBSim simulation executive.
FGPropagate * GetPropagate(void)
Returns the FGPropagate pointer.