COSMOS core  1.0.2 (beta)
Comprehensive Open-architecture Solution for Mission Operations Systems
Functions | Variables
nrlmsise-00.cpp File Reference
#include "physics/nrlmsise-00.h"
#include <cmath>
#include <stdio.h>
#include <stdlib.h>
Include dependency graph for nrlmsise-00.cpp:

Functions

void tselec (struct nrlmsise_flags *flags)
 
void glatf (double lat, double *gv, double *reff)
 
double ccor (double alt, double r, double h1, double zh)
 
double ccor2 (double alt, double r, double h1, double zh, double h2)
 
double scalh (double alt, double xm, double temp)
 
double dnet (double dd, double dm, double zhm, double xmm, double xm)
 
void splini (double *xa, double *ya, double *y2a, int n, double x, double *y)
 
void splint (double *xa, double *ya, double *y2a, int n, double x, double *y)
 
void spline (double *x, double *y, int n, double yp1, double ypn, double *y2)
 
__inline_double zeta (double zz, double zl)
 
double densm (double alt, double d0, double xm, double &tz, int mn3, double *zn3, double *tn3, double *tgn3, int mn2, double *zn2, double *tn2, double *tgn2)
 
double densu (double alt, double dlb, double tinf, double tlb, double xm, double alpha, double *tz, double zlb, double s2, int mn1, double *zn1, double *tn1, double *tgn1)
 
__inline_double g0 (double a, double *p)
 
__inline_double sumex (double ex)
 
__inline_double sg0 (double ex, double *p, double *ap)
 
double globe7 (double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags)
 
double glob7s (double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags)
 
void gtd7 (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
 
void gtd7d (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
 
void ghp7 (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output, double press)
 
void gts7 (struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
 

Variables

static double gsurf
 
static double re
 
static double dd
 
static double dm04
 
static double dm16
 
static double dm28
 
static double dm32
 
static double dm40
 
static double dm01
 
static double dm14
 
static double meso_tn1 [5]
 
static double meso_tn2 [4]
 
static double meso_tn3 [5]
 
static double meso_tgn1 [2]
 
static double meso_tgn2 [2]
 
static double meso_tgn3 [2]
 
double pt [150]
 
double pd [9][150]
 
double ps [150]
 
double pdl [2][25]
 
double ptl [4][100]
 
double pma [10][100]
 
double sam [100]
 
double ptm [10]
 
double pdm [8][10]
 
double pavgm [10]
 
static double dfa
 
static double plg [4][9]
 
static double ctloc
 
static double stloc
 
static double c2tloc
 
static double s2tloc
 
static double s3tloc
 
static double c3tloc
 
static double apdf
 
static double apt [4]
 

Function Documentation

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

Variable Documentation

double gsurf
static
double re
static
double dd
static
double dm04
static
double dm16
static
double dm28
static
double dm32
static
double dm40
static
double dm01
static
double dm14
static
double meso_tn1[5]
static
double meso_tn2[4]
static
double meso_tn3[5]
static
double meso_tgn1[2]
static
double meso_tgn2[2]
static
double meso_tgn3[2]
static
double pt[150]
double pd[9][150]
double ps[150]
double pdl[2][25]
double ptl[4][100]
double pma[10][100]
double sam[100]
double ptm[10]
double pdm[8][10]
double pavgm[10]
double dfa
static
double plg[4][9]
static
double ctloc
static
double stloc
static
double c2tloc
static
double s2tloc
static
double s3tloc
static
double c3tloc
static
double apdf
static
double apt[4]
static