attachment:nrlmsise-00.c of HitchhikerReentry

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