diff options
| author | Dominik Brodowski <linux@dominikbrodowski.net> | 2002-03-05 12:00:00 +0100 |
|---|---|---|
| committer | Dominik Brodowski <linux@dominikbrodowski.net> | 2017-09-24 10:21:32 +0200 |
| commit | 1daf607ecfd464fe121da3feb6aed020eed210e6 (patch) | |
| tree | 8489abd28e985372c58f802dc13fcc89252f27dd | |
| download | nrlmsise-00-1daf607ecfd464fe121da3feb6aed020eed210e6.tar.gz | |
NRLMSISE-00 initial C release
C source code for the NRLMSISE-00 empirical atmosphere model
The NRLMSIS-00 (sic!) empirical atmosphere model was developed
by Mike Picone, Alan Hedin, and Doug Drob. It describes the
neutral temperature and densities in Earth's atmosphere from
ground to thermospheric heights.
(quoted from http://modelweb.gsfc.nasa.gov/ ) - you can find
a longer explanation of this model there, too.
The authors of NRLMSISE-00 have released a FORTRAN version
which is available at
http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm.
Based on the Official Beta Release 1.0 (NRLMSISE-00.DIST12.TXT)
Dominik Brodowski wrote an implementation in C.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. Please
inform the maintainer of the C release (Dominik Brodowski -
mail@brodo.de) of any patches and bug-fixes you implement for
NRLMSISE-00 so that this C package can be updated with these
improvements.
Signed-off-by: Dominik Brodowski <linux@dominikbrodwski.net>
| -rw-r--r-- | DOCUMENTATION | 250 | ||||
| -rw-r--r-- | makefile | 21 | ||||
| -rw-r--r-- | nrlmsise-00.c | 1448 | ||||
| -rw-r--r-- | nrlmsise-00.h | 222 | ||||
| -rw-r--r-- | nrlmsise-00_data.c | 740 | ||||
| -rw-r--r-- | nrlmsise-00_test.c | 168 |
6 files changed, 2849 insertions, 0 deletions
diff --git a/DOCUMENTATION b/DOCUMENTATION new file mode 100644 index 0000000..5562a2c --- /dev/null +++ b/DOCUMENTATION @@ -0,0 +1,250 @@ +-------------------------------------------------------------------- +--------- N R L M S I S E - 0 0 M O D E L 2 0 0 1 ---------- +-------------------------------------------------------------------- + + +Table of Contents +----------------- + +1. Legal Information +2. Brief Description +3. Source Code Availability +4. This C Release + 4.1 Files + 4.2 Differences between FORTRAN and C version +5. Interface +6. Release Notes +7. Testing Output + + + +1. LEGAL INFORMATION +==================== + +This package is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. Please inform the +maintainer of the C release (Dominik Brodowski - devel@brodo.de) of +any patches and bug-fixes you implement for NRLMSISE-00 so that this C +package can be updated with these improvements. + + + +2. BRIEF DESCRIPTION +==================== + +The NRLMSIS-00 empirical atmosphere model was developed by Mike +Picone, Alan Hedin, and Doug Drob based on the MSISE90 model. + +The MSISE90 model describes the neutral temperature and densities in +Earth's atmosphere from ground to thermospheric heights. Below 72.5 km +the model is primarily based on the MAP Handbook (Labitzke et al., +1985) tabulation of zonal average temperature and pressure by Barnett +and Corney, which was also used for the CIRA-86. Below 20 km these +data were supplemented with averages from the National Meteorological +Center (NMC). In addition, pitot tube, falling sphere, and grenade +sounder rocket measurements from 1947 to 1972 were taken into +consideration. Above 72.5 km MSISE-90 is essentially a revised MSIS-86 +model taking into account data derived from space shuttle flights and +newer incoherent scatter results. For someone interested only in the +thermosphere (above 120 km), the author recommends the MSIS-86 +model. MSISE is also not the model of preference for specialized +tropospheric work. It is rather for studies that reach across several +atmospheric boundaries. +(quoted from http://nssdc.gsfc.nasa.gov/space/model/atmos/nrlmsise00.html) + + + +3. SOURCE CODE AVAILABILITY +=========================== + +The authors of the NRLMSISE-00 model have released a FORTRAN version +which is available at +http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm + +Based on the Official Beta Release 1.0 (NRLMSISE-00.DIST12.TXT) +Dominik Brodowski wrote an implementation in C which is available on +http://www.brodo.de/english/pub/nrlmsise/ +Please check for updates often during the next weeks, since the C +package probably has some bugs. Additionally, some modifications +to the NRLMSISE-00 FORTRAN package will be published soon, and I'll +try to update the C package accordingly within a few days. + + + +4. THIS C RELEASE +================= + +When "INLINE" is not externally defined the source code should be +clean, plain ANSI which any C compiler should be able to work with. If +you find any problems, please report them to the maintainer Dominik +Brodowski at devel@brodo.de. Thank you! + + +4.1 Files +--------- + +DOCUMENTATION - this Documentation +nrlmsise-00.h - header file for nrlmsise-00 +nrlmsise-00.c - source code for nrlmsise-00 +nrlmsise-00_data.c - block data for nrlmsise-00 +nrlmsise-00_test.c - test-driver for nrlmsise-00 +makefile - makefile for nrlmsise-test (gnu-make and gcc) + + +4.2 Differences between FORTRAN and C version +--------------------------------------------- + +The C package does not save the last results internally to speed up +program execution as the FORTRAN version does. + +The "switches" have to be specified before _every_ call. + +The "DL" array is not printed in the testing routine since it's not an +output value. + +The C version probably contains some bugs and still has to be regarded +more unstable than the FORTRAN relase. Please report any bugs or +incorrect values to the maintainer at devel@brodo.de + + + +5. INTERFACE +============ + +To access the NRLMSISE-00 functions you need to include the header +file nrlmsise-00.h. In this file you can find comments which explain +the in- and output, and the differences between the various functions. + + + +6. RELEASE NOTES +================ + +20020302 - first release + + + +7. TESTING OUTPUT +================= + +nrlmsise-test should generate the following output: + + +6.665177E+05 1.138806E+08 1.998211E+07 3.964051E+05 3.557465E+03 4.074402E-15 3.475312E+04 4.095913E+06 2.667273E+04 1.250540E+03 1.241416E+03 + +3.407293E+06 1.586333E+08 1.391117E+07 3.214942E+05 1.559618E+03 5.001593E-15 4.854208E+04 4.380967E+06 6.956682E+03 1.166754E+03 1.161710E+03 + +1.123767E+05 6.934130E+04 4.247105E+01 1.303616E-01 2.618848E-05 2.756772E-18 2.016750E+04 5.741256E+03 2.374394E+04 1.239892E+03 1.239891E+03 + +7.078194E+07 1.946309E+11 6.115826E+12 1.183637E+12 6.023212E+10 3.563077E-10 1.059880E+07 2.615737E+05 2.819879E-42 1.027318E+03 2.068878E+02 + +1.851122E+06 1.476555E+08 1.579356E+07 2.595354E+05 1.588781E+03 4.809426E-15 5.816167E+04 5.478984E+06 1.264446E+03 1.212396E+03 1.208135E+03 + +8.673095E+05 1.278862E+08 1.822577E+07 2.879564E+05 2.402962E+03 4.355639E-15 3.686389E+04 3.897276E+06 2.667273E+04 1.220146E+03 1.212712E+03 + +5.776251E+05 6.979139E+07 1.236814E+07 2.456484E+05 1.405739E+03 2.470458E-15 5.291986E+04 1.069814E+06 2.667273E+04 1.116385E+03 1.112999E+03 + +3.740304E+05 4.782720E+07 5.240380E+06 1.711538E+05 5.501649E+02 1.571632E-15 8.896776E+04 1.979741E+06 9.121815E+03 1.031247E+03 1.024848E+03 + +6.748339E+05 1.245315E+08 2.369010E+07 4.839898E+05 4.578781E+03 4.564039E-15 3.244595E+04 5.370833E+06 2.667273E+04 1.306052E+03 1.293374E+03 + +5.528601E+05 1.198041E+08 3.495798E+07 9.203305E+05 1.096255E+04 4.973819E-15 2.686428E+04 4.889974E+06 2.805445E+04 1.361868E+03 1.347389E+03 + +1.375488E+14 0.000000E+00 2.049687E+19 5.498695E+18 2.451733E+17 1.261066E-03 0.000000E+00 0.000000E+00 0.000000E+00 1.027318E+03 2.814648E+02 + +4.427443E+13 0.000000E+00 6.597567E+18 1.769929E+18 7.891680E+16 4.059139E-04 0.000000E+00 0.000000E+00 0.000000E+00 1.027318E+03 2.274180E+02 + +2.127829E+12 0.000000E+00 3.170791E+17 8.506280E+16 3.792741E+15 1.950822E-05 0.000000E+00 0.000000E+00 0.000000E+00 1.027318E+03 2.374389E+02 + +1.412184E+11 0.000000E+00 2.104370E+16 5.645392E+15 2.517142E+14 1.294709E-06 0.000000E+00 0.000000E+00 0.000000E+00 1.027318E+03 2.795551E+02 + +1.261027E+10 0.000000E+00 1.874533E+15 4.809049E+14 2.239685E+13 1.141612E-07 0.000000E+00 0.000000E+00 0.000000E+00 1.027318E+03 2.190732E+02 + +5.196477E+05 1.274494E+08 4.850450E+07 1.695722E+06 2.354487E+04 5.880606E-15 2.500078E+04 6.279210E+06 2.667273E+04 1.426412E+03 1.408608E+03 + +5.673929E+07 1.259077E+11 4.929562E+12 1.021646E+12 4.993465E+10 2.900560E-10 8.831229E+06 2.252516E+05 2.415246E-42 1.027318E+03 1.934071E+02 + + +DAY 172 81 172 172 172 +UT 29000 29000 75000 29000 29000 +ALT 400 400 1000 100 400 +LAT 60 60 60 60 0 +LONG -70 -70 -70 -70 -70 +LST 16 16 16 16 16 +F107A 150 150 150 150 150 +F107 150 150 150 150 150 + + +TINF 1250.54 1166.75 1239.89 1027.32 1212.40 +TG 1241.42 1161.71 1239.89 206.89 1208.14 +HE 6.665e+05 3.407e+06 1.124e+05 7.078e+07 1.851e+06 +O 1.139e+08 1.586e+08 6.934e+04 1.946e+11 1.477e+08 +N2 1.998e+07 1.391e+07 4.247e+01 6.116e+12 1.579e+07 +O2 3.964e+05 3.215e+05 1.304e-01 1.184e+12 2.595e+05 +AR 3.557e+03 1.560e+03 2.619e-05 6.023e+10 1.589e+03 +H 3.475e+04 4.854e+04 2.017e+04 1.060e+07 5.816e+04 +N 4.096e+06 4.381e+06 5.741e+03 2.616e+05 5.479e+06 +ANM 0 2.667e+04 6.957e+03 2.374e+04 2.820e-42 1.264e+03 +RHO 4.074e-15 5.002e-15 2.757e-18 3.563e-10 4.809e-15 + + +DAY 172 172 172 172 172 +UT 29000 29000 29000 29000 29000 +ALT 400 400 400 400 400 +LAT 60 60 60 60 60 +LONG 0 -70 -70 -70 -70 +LST 16 4 16 16 16 +F107A 150 150 70 150 150 +F107 150 150 150 180 150 + + +TINF 1220.15 1116.39 1031.25 1306.05 1361.87 +TG 1212.71 1113.00 1024.85 1293.37 1347.39 +HE 8.673e+05 5.776e+05 3.740e+05 6.748e+05 5.529e+05 +O 1.279e+08 6.979e+07 4.783e+07 1.245e+08 1.198e+08 +N2 1.823e+07 1.237e+07 5.240e+06 2.369e+07 3.496e+07 +O2 2.880e+05 2.456e+05 1.712e+05 4.840e+05 9.203e+05 +AR 2.403e+03 1.406e+03 5.502e+02 4.579e+03 1.096e+04 +H 3.686e+04 5.292e+04 8.897e+04 3.245e+04 2.686e+04 +N 3.897e+06 1.070e+06 1.980e+06 5.371e+06 4.890e+06 +ANM 0 2.667e+04 2.667e+04 9.122e+03 2.667e+04 2.805e+04 +RHO 4.356e-15 2.470e-15 1.572e-15 4.564e-15 4.974e-15 + + +DAY 172 172 172 172 172 +UT 29000 29000 29000 29000 29000 +ALT 0 10 30 50 70 +LAT 60 60 60 60 60 +LONG -70 -70 -70 -70 -70 +LST 16 16 16 16 16 +F107A 150 150 150 150 150 +F107 150 150 150 150 150 + + +TINF 1027.32 1027.32 1027.32 1027.32 1027.32 +TG 281.46 227.42 237.44 279.56 219.07 +HE 1.375e+14 4.427e+13 2.128e+12 1.412e+11 1.261e+10 +O 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 +N2 2.050e+19 6.598e+18 3.171e+17 2.104e+16 1.875e+15 +O2 5.499e+18 1.770e+18 8.506e+16 5.645e+15 4.809e+14 +AR 2.452e+17 7.892e+16 3.793e+15 2.517e+14 2.240e+13 +H 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 +N 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 +ANM 0 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 +RHO 1.261e-03 4.059e-04 1.951e-05 1.295e-06 1.142e-07 + + +Note: These values equal those of the official FORTRAN package with +one notable exception: the FORTRAN version reports for "anomalous +oxygen" in test-runs 4 and 17 exactly 0.000E-00, while my C compiler +generates code which calculates 2.820e-42. When only 16-bit wide +double variables are used, this value reduces to 0.000E-00 as well. + + + +======================================================== + +Rosenheim, Germany, on March 5th, 2002 + +Dominik Brodowski diff --git a/makefile b/makefile new file mode 100644 index 0000000..fa2a07e --- /dev/null +++ b/makefile @@ -0,0 +1,21 @@ +CFLAGS = -Wall -g -DINLINE +MFLAGS = $(CFLAGS) -lm +CC = gcc $(CFLAGS) +MCC = gcc $(MFLAGS) + +nrlmsise-test : nrlmsise-00.o nrlmsise-00_test.o nrlmsise-00_data.o + $(MCC) -o nrlmsise-test nrlmsise-00.o nrlmsise-00_test.o \ + nrlmsise-00_data.o + +nrlmsise-00.o : nrlmsise-00.c nrlmsise-00.h + $(CC) -c nrlmsise-00.c + +nrlmsise-00_test.o : nrlmsise-00_test.c nrlmsise-00.h + $(CC) -c nrlmsise-00_test.c + +nrlmsise-00_data.o : nrlmsise-00_data.c nrlmsise-00.h + $(CC) -c nrlmsise-00_data.c + +clean : + rm -rf nrlmsise-test nrlmsise-00.o nrlmsise-00_test.o \ + nrlmsise-00_data.o diff --git a/nrlmsise-00.c b/nrlmsise-00.c new file mode 100644 index 0000000..73a8062 --- /dev/null +++ b/nrlmsise-00.c @@ -0,0 +1,1448 @@ +/* -------------------------------------------------------------------- */ +/* --------- N R L M S I S E - 0 0 M O D E L 2 0 0 1 ---------- */ +/* -------------------------------------------------------------------- */ + +/* This file is part of the NRLMSISE-00 C source code package - release + * 20020305 + * + * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and + * Doug Drob. They also wrote a NRLMSISE-00 distribution package in + * FORTRAN which is available at + * http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm + * + * Dominik Brodowski implemented and maintains this C version. You can + * reach him at devel@brodo.de. See the file "DOCUMENTATION" for details, + * and check http://www.brodo.de/english/pub/nrlmsise/index.html for + * updated releases of this package. + */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ INCLUDES --------------------------- */ +/* ------------------------------------------------------------------- */ + +#include "nrlmsise-00.h" /* header for nrlmsise-00.h */ +#include <math.h> /* maths functions */ +#include <stdio.h> /* for error messages. TBD: remove this */ +#include <stdlib.h> /* for malloc/free */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------- SHARED VARIABLES ------------------------ */ +/* ------------------------------------------------------------------- */ + +/* PARMB */ +static double gsurf; +static double re; + +/* GTS3C */ +static double dd; + +/* DMIX */ +static double dm04, dm16, dm28, dm32, dm40, dm01, dm14; + +/* MESO7 */ +static double meso_tn1[5]; +static double meso_tn2[4]; +static double meso_tn3[3]; +static double meso_tgn1[2]; +static double meso_tgn2[2]; +static double meso_tgn3[2]; + +/* POWER7 */ +extern double pt[150]; +extern double pd[9][150]; +extern double ps[150]; +extern double pdl[2][25]; +extern double ptl[4][100]; +extern double pma[10][100]; +extern double sam[100]; + +/* LOWER7 */ +extern double ptm[10]; +extern double pdm[8][10]; +extern double pavgm[10]; + +/* LPOLY */ +static double dfa; +static double plg[4][9]; +static double ctloc, stloc; +static double c2tloc, s2tloc; +static double s3tloc, c3tloc; +static double apdf, apt[4]; + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ TSELEC ----------------------------- */ +/* ------------------------------------------------------------------- */ + +void tselec(struct nrlmsise_flags *flags) { + int i; + for (i=0;i<24;i++) { + if (i!=9) { + if (flags->switches[i]==1) + flags->sw[i]=1; + else + flags->sw[i]=0; + if (flags->switches[i]>0) + flags->swc[i]=1; + else + flags->swc[i]=0; + } else { + flags->sw[i]=flags->switches[i]; + flags->swc[i]=flags->switches[i]; + } + } +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ GLATF ------------------------------ */ +/* ------------------------------------------------------------------- */ + +void glatf(double lat, double *gv, double *reff) { + double dgtr = 1.74533E-2; + double c2; + c2 = cos(2.0*dgtr*lat); + *gv = 980.616 * (1.0 - 0.0026373 * c2); + *reff = 2.0 * (*gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ CCOR ------------------------------- */ +/* ------------------------------------------------------------------- */ + +double ccor(double alt, double r, double h1, double zh) { +/* CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS + * ALT - altitude + * R - target ratio + * H1 - transition scale length + * ZH - altitude of 1/2 R + */ + double e; + double ex; + e = (alt - zh) / h1; + if (e>70) + return exp(0); + if (e<-70) + return exp(r); + ex = exp(e); + e = r / (1.0 + ex); + return exp(e); +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- SCALH ----------------------------- */ +/* ------------------------------------------------------------------- */ + +double scalh(double alt, double xm, double temp) { + double g; + double rgas=831.4; + g = gsurf / (pow((1.0 + alt/re),2.0)); + g = rgas * temp / (g * xm); + return g; +} + + + +/* ------------------------------------------------------------------- */ +/* -------------------------------- DNET ----------------------------- */ +/* ------------------------------------------------------------------- */ + +double dnet (double dd, double dm, double zhm, double xmm, double xm) { +/* TURBOPAUSE CORRECTION FOR MSIS MODELS + * Root mean density + * DD - diffusive density + * DM - full mixed density + * ZHM - transition scale length + * XMM - full mixed molecular weight + * XM - species molecular weight + * DNET - combined density + */ + double a; + double ylog; + a = zhm / (xmm-xm); + if (!((dm>0) && (dd>0))) { + printf("dnet log error %e %e %e\n",dm,dd,xm); + if ((dd==0) && (dm==0)) + dd=1; + if (dm==0) + return dd; + if (dd==0) + return dm; + } + ylog = a * log(dm/dd); + if (ylog<-10) + return dd; + if (ylog>10) + return dm; + a = dd*pow((1.0 + exp(ylog)),(1.0/a)); + return a; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- SPLINI ---------------------------- */ +/* ------------------------------------------------------------------- */ + +void splini (double *xa, double *ya, double *y2a, int n, double x, double *y) { +/* INTEGRATE CUBIC SPLINE FUNCTION FROM XA(1) TO X + * XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X + * Y2A: ARRAY OF SECOND DERIVATIVES + * N: SIZE OF ARRAYS XA,YA,Y2A + * X: ABSCISSA ENDPOINT FOR INTEGRATION + * Y: OUTPUT VALUE + */ + double yi=0; + int klo=0; + int khi=1; + double xx, h, a, b, a2, b2; + while ((x>xa[klo]) && (khi<n)) { + xx=x; + if (khi<(n-1)) { + if (x<xa[khi]) + xx=x; + else + xx=xa[khi]; + } + h = xa[khi] - xa[klo]; + a = (xa[khi] - xx)/h; + b = (xx - xa[klo])/h; + a2 = a*a; + b2 = b*b; + 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; + klo++; + khi++; + } + *y = yi; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- SPLINT ---------------------------- */ +/* ------------------------------------------------------------------- */ + +void splint (double *xa, double *ya, double *y2a, int n, double x, double *y) { +/* CALCULATE CUBIC SPLINE INTERP VALUE + * ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL. + * XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X + * Y2A: ARRAY OF SECOND DERIVATIVES + * N: SIZE OF ARRAYS XA,YA,Y2A + * X: ABSCISSA FOR INTERPOLATION + * Y: OUTPUT VALUE + */ + int klo=0; + int khi=n-1; + int k; + double h; + double a, b, yi; + while ((khi-klo)>1) { + k=(khi+klo)/2; + if (xa[k]>x) + khi=k; + else + klo=k; + } + h = xa[khi] - xa[klo]; + if (h==0.0) + printf("bad XA input to splint"); + a = (xa[khi] - x)/h; + b = (x - xa[klo])/h; + yi = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0; + *y = yi; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- SPLINE ---------------------------- */ +/* ------------------------------------------------------------------- */ + +void spline (double *x, double *y, int n, double yp1, double ypn, double *y2) { +/* CALCULATE 2ND DERIVATIVES OF CUBIC SPLINE INTERP FUNCTION + * ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL + * X,Y: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X + * N: SIZE OF ARRAYS X,Y + * YP1,YPN: SPECIFIED DERIVATIVES AT X[0] AND X[N-1]; VALUES + * >= 1E30 SIGNAL SIGNAL SECOND DERIVATIVE ZERO + * Y2: OUTPUT ARRAY OF SECOND DERIVATIVES + */ + double *u; + double sig, p, qn, un; + int i, k; + u=malloc(sizeof(double)*n); + if (u==NULL) { + printf("Out Of Memory in spline - ERROR"); + return; + } + if (yp1>0.99E30) { + y2[0]=0; + u[0]=0; + } else { + y2[0]=-0.5; + u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1); + } + for (i=1;i<(n-1);i++) { + sig = (x[i]-x[i-1])/(x[i+1] - x[i-1]); + p = sig * y2[i-1] + 2.0; + y2[i] = (sig - 1.0) / p; + 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; + } + if (ypn>0.99E30) { + qn = 0; + un = 0; + } else { + qn = 0.5; + un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])/(x[n-1] - x[n-2])); + } + y2[n-1] = (un - qn * u[n-2]) / (qn * y2[n-2] + 1.0); + for (k=n-2;k>=0;k--) + y2[k] = y2[k] * y2[k+1] + u[k]; + + free(u); +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- DENSM ----------------------------- */ +/* ------------------------------------------------------------------- */ + +__inline_double zeta(double zz, double zl) { + return ((zz-zl)*(re+zl)/(re+zz)); +} + +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) { +/* Calculate Temperature and Density Profiles for lower atmos. */ + double xs[10], ys[10], y2out[10]; + double rgas = 831.4; + double z, z1, z2, t1, t2, zg, zgdif; + double yd1, yd2; + double x, y, yi; + double expl, gamm, glb; + double densm_tmp; + int mn; + int k; + densm_tmp=d0; + if (alt>zn2[0]) { + if (xm==0.0) + return *tz; + else + return d0; + } + + /* STRATOSPHERE/MESOSPHERE TEMPERATURE */ + if (alt>zn2[mn2-1]) + z=alt; + else + z=zn2[mn2-1]; + mn=mn2; + z1=zn2[0]; + z2=zn2[mn-1]; + t1=tn2[0]; + t2=tn2[mn-1]; + zg = zeta(z, z1); + zgdif = zeta(z2, z1); + + /* set up spline nodes */ + for (k=0;k<mn;k++) { + xs[k]=zeta(zn2[k],z1)/zgdif; + ys[k]=1.0 / tn2[k]; + } + yd1=-tgn2[0] / (t1*t1) * zgdif; + yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0)); + + /* calculate spline coefficients */ + spline (xs, ys, mn, yd1, yd2, y2out); + x = zg/zgdif; + splint (xs, ys, y2out, mn, x, &y); + + /* temperature at altitude */ + *tz = 1.0 / y; + if (xm!=0.0) { + /* calaculate stratosphere / mesospehere density */ + glb = gsurf / (pow((1.0 + z1/re),2.0)); + gamm = xm * glb * zgdif / rgas; + + /* Integrate temperature profile */ + splini(xs, ys, y2out, mn, x, &yi); + expl=gamm*yi; + if (expl>50.0) + expl=50.0; + + /* Density at altitude */ + densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl); + } + + if (alt>zn3[0]) { + if (xm==0.0) + return *tz; + else + return densm_tmp; + } + + /* troposhere / stratosphere temperature */ + z = alt; + mn = mn3; + z1=zn3[0]; + z2=zn3[mn-1]; + t1=tn3[0]; + t2=tn3[mn-1]; + zg=zeta(z,z1); + zgdif=zeta(z2,z1); + + /* set up spline nodes */ + for (k=0;k<mn;k++) { + xs[k] = zeta(zn3[k],z1) / zgdif; + ys[k] = 1.0 / tn3[k]; + } + yd1=-tgn3[0] / (t1*t1) * zgdif; + yd2=-tgn3[1] / (t2*t2) * zgdif * (pow(((re+z2)/(re+z1)),2.0)); + + /* calculate spline coefficients */ + spline (xs, ys, mn, yd1, yd2, y2out); + x = zg/zgdif; + splint (xs, ys, y2out, mn, x, &y); + + /* temperature at altitude */ + *tz = 1.0 / y; + if (xm!=0.0) { + /* calaculate tropospheric / stratosphere density */ + glb = gsurf / (pow((1.0 + z1/re),2.0)); + gamm = xm * glb * zgdif / rgas; + + /* Integrate temperature profile */ + splini(xs, ys, y2out, mn, x, &yi); + expl=gamm*yi; + if (expl>50.0) + expl=50.0; + + /* Density at altitude */ + densm_tmp = densm_tmp * (t1 / *tz) * exp(-expl); + } + if (xm==0.0) + return *tz; + else + return densm_tmp; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- 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) { +/* Calculate Temperature and Density Profiles for MSIS models + * New lower thermo polynomial + */ + double yd2, yd1, x, y; + double rgas=831.4; + double densu_temp=1.0; + double za, z, zg2, tt, ta; + double dta, z1, z2, t1, t2, zg, zgdif; + int mn; + int k; + double glb; + double expl; + double yi; + double densa; + double gamma, gamm; + double xs[5], ys[5], y2out[5]; + /* joining altitudes of Bates and spline */ + za=zn1[0]; + if (alt>za) + z=alt; + else + z=za; + + /* geopotential altitude difference from ZLB */ + zg2 = zeta(z, zlb); + + /* Bates temperature */ + tt = tinf - (tinf - tlb) * exp(-s2*zg2); + ta = tt; + *tz = tt; + densu_temp = *tz; + + if (alt<za) { + /* calculate temperature below ZA + * temperature gradient at ZA from Bates profile */ + dta = (tinf - ta) * s2 * pow(((re+zlb)/(re+za)),2.0); + tgn1[0]=dta; + tn1[0]=ta; + if (alt>zn1[mn1-1]) + z=alt; + else + z=zn1[mn1-1]; + mn=mn1; + z1=zn1[0]; + z2=zn1[mn-1]; + t1=tn1[0]; + t2=tn1[mn-1]; + /* geopotental difference from z1 */ + zg = zeta (z, z1); + zgdif = zeta(z2, z1); + /* set up spline nodes */ + for (k=0;k<mn;k++) { + xs[k] = zeta(zn1[k], z1) / zgdif; + ys[k] = 1.0 / tn1[k]; + } + /* end node derivatives */ + yd1 = -tgn1[0] / (t1*t1) * zgdif; + yd2 = -tgn1[1] / (t2*t2) * zgdif * pow(((re+z2)/(re+z1)),2.0); + /* calculate spline coefficients */ + spline (xs, ys, mn, yd1, yd2, y2out); + x = zg / zgdif; + splint (xs, ys, y2out, mn, x, &y); + /* temperature at altitude */ + *tz = 1.0 / y; + densu_temp = *tz; + } + if (xm==0) + return densu_temp; + + /* calculate density above za */ + glb = gsurf / pow((1.0 + zlb/re),2.0); + gamma = xm * glb / (s2 * rgas * tinf); + expl = exp(-s2 * gamma * zg2); + if (expl>50.0) + expl=50.0; + if (tt<=0) + expl=50.0; + + /* density at altitude */ + densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl; + densu_temp=densa; + if (alt>=za) + return densu_temp; + + /* calculate density below za */ + glb = gsurf / pow((1.0 + z1/re),2.0); + gamm = xm * glb * zgdif / rgas; + + /* integrate spline temperatures */ + splini (xs, ys, y2out, mn, x, &yi); + expl = gamm * yi; + if (expl>50.0) + expl=50.0; + if (*tz<=0) + expl=50.0; + + /* density at altitude */ + densu_temp = densu_temp * pow ((t1 / *tz),(1.0 + alpha)) * exp(-expl); + return densu_temp; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GLOBE7 ---------------------------- */ +/* ------------------------------------------------------------------- */ + +/* 3hr Magnetic activity functions */ +/* Eq. A24d */ +__inline_double g0(double a, double *p) { + 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]))); +} + +/* Eq. A24c */ +__inline_double sumex(double ex) { + return (1.0 + (1.0 - pow(ex,19.0)) / (1.0 - ex) * pow(ex,0.5)); +} + +/* Eq. A24a */ +__inline_double sg0(double ex, double *p, double *ap) { + return (g0(ap[1],p) + (g0(ap[2],p)*ex + g0(ap[3],p)*ex*ex + \ + g0(ap[4],p)*pow(ex,3.0) + (g0(ap[5],p)*pow(ex,4.0) + \ + g0(ap[6],p)*pow(ex,12.0))*(1.0-pow(ex,8.0))/(1.0-ex)))/sumex(ex); +} + +double globe7(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) { +/* CALCULATE G(L) FUNCTION + * Upper Thermosphere Parameters */ + double t[15]; + int i,j; + int sw9=1; + double apd; + double xlong; + double tloc; + double c, s, c2, c4, s2; + double sr = 7.2722E-5; + double dgtr = 1.74533E-2; + double dr = 1.72142E-2; + double hr = 0.2618; + double cd32, cd18, cd14, cd39; + double p32, p18, p14, p39; + double df, dfa; + double f1, f2; + double tinf; + struct ap_array *ap; + + tloc=input->lst; + for (j=0;j<14;j++) + t[j]=0; + if (flags->sw[9]>0) + sw9=1; + else if (flags->sw[9]<0) + sw9=-1; + xlong = input->g_long; + + /* calculate legendre polynomials */ + c = sin(input->g_lat * dgtr); + s = cos(input->g_lat * dgtr); + c2 = c*c; + c4 = c2*c2; + s2 = s*s; + + plg[0][1] = c; + plg[0][2] = 0.5*(3.0*c2 -1.0); + plg[0][3] = 0.5*(5.0*c*c2-3.0*c); + plg[0][4] = (35.0*c4 - 30.0*c2 + 3.0)/8.0; + plg[0][5] = (63.0*c2*c2*c - 70.0*c2*c + 15.0*c)/8.0; + plg[0][6] = (11.0*c*plg[0][5] - 5.0*plg[0][4])/6.0; +/* plg[0][7] = (13.0*c*plg[0][6] - 6.0*plg[0][5])/7.0; */ + plg[1][1] = s; + plg[1][2] = 3.0*c*s; + plg[1][3] = 1.5*(5.0*c2-1.0)*s; + plg[1][4] = 2.5*(7.0*c2*c-3.0*c)*s; + plg[1][5] = 1.875*(21.0*c4 - 14.0*c2 +1.0)*s; + plg[1][6] = (11.0*c*plg[1][5]-6.0*plg[1][4])/5.0; +/* plg[1][7] = (13.0*c*plg[1][6]-7.0*plg[1][5])/6.0; */ +/* plg[1][8] = (15.0*c*plg[1][7]-8.0*plg[1][6])/7.0; */ + plg[2][2] = 3.0*s2; + plg[2][3] = 15.0*s2*c; + plg[2][4] = 7.5*(7.0*c2 -1.0)*s2; + plg[2][5] = 3.0*c*plg[2][4]-2.0*plg[2][3]; + plg[2][6] =(11.0*c*plg[2][5]-7.0*plg[2][4])/4.0; + plg[2][7] =(13.0*c*plg[2][6]-8.0*plg[2][5])/5.0; + plg[3][3] = 15.0*s2*s; + plg[3][4] = 105.0*s2*s*c; + plg[3][5] =(9.0*c*plg[3][4]-7.*plg[3][3])/2.0; + plg[3][6] =(11.0*c*plg[3][5]-8.*plg[3][4])/3.0; + + if (!(((flags->sw[7]==0)&&(flags->sw[8]==0))&&(flags->sw[14]==0))) { + stloc = sin(hr*tloc); + ctloc = cos(hr*tloc); + s2tloc = sin(2.0*hr*tloc); + c2tloc = cos(2.0*hr*tloc); + s3tloc = sin(3.0*hr*tloc); + c3tloc = cos(3.0*hr*tloc); + } + + cd32 = cos(dr*(input->doy-p[31])); + cd18 = cos(2.0*dr*(input->doy-p[17])); + cd14 = cos(dr*(input->doy-p[13])); + cd39 = cos(2.0*dr*(input->doy-p[38])); + p32=p[31]; + p18=p[17]; + p14=p[13]; + p39=p[38]; + + /* F10.7 EFFECT */ + df = input->f107 - input->f107A; + dfa = input->f107A - 150.0; + t[0] = p[19]*df*(1.0+p[59]*dfa) + p[20]*df*df + p[21]*dfa + p[29]*pow(dfa,2.0); + f1 = 1.0 + (p[47]*dfa +p[19]*df+p[20]*df*df)*flags->swc[1]; + f2 = 1.0 + (p[49]*dfa+p[19]*df+p[20]*df*df)*flags->swc[1]; + + /* TIME INDEPENDENT */ + t[1] = (p[1]*plg[0][2]+ p[2]*plg[0][4]+p[22]*plg[0][6]) + \ + (p[14]*plg[0][2])*dfa*flags->swc[1] +p[26]*plg[0][1]; + + /* SYMMETRICAL ANNUAL */ + t[2] = p[18]*cd32; + + /* SYMMETRICAL SEMIANNUAL */ + t[3] = (p[15]+p[16]*plg[0][2])*cd18; + + /* ASYMMETRICAL ANNUAL */ + t[4] = f1*(p[9]*plg[0][1]+p[10]*plg[0][3])*cd14; + + /* ASYMMETRICAL SEMIANNUAL */ + t[5] = p[37]*plg[0][1]*cd39; + + /* DIURNAL */ + if (flags->sw[7]) { + double t71, t72; + t71 = (p[11]*plg[1][2])*cd14*flags->swc[5]; + t72 = (p[12]*plg[1][2])*cd14*flags->swc[5]; + t[6] = f2*((p[3]*plg[1][1] + p[4]*plg[1][3] + p[27]*plg[1][5] + t71) * \ + ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + p[28]*plg[1][5] \ + + t72)*stloc); +} + + /* SEMIDIURNAL */ + if (flags->sw[8]) { + double t81, t82; + t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5]; + t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5]; + 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); + } + + /* TERDIURNAL */ + if (flags->sw[14]) { + 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); +} + + /* magnetic activity based on daily ap */ + if (flags->sw[9]==-1) { + ap = input->ap_a; + if (p[51]!=0) { + double exp1; + exp1 = exp(-10800.0*sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-sqrt(input->g_lat*input->g_lat)))); + if (exp1>0.99999) + exp1=0.99999; + if (p[24]<1.0E-4) + p[24]=1.0E-4; + apt[0]=sg0(exp1,p,ap->a); + /* apt[1]=sg2(exp1,p,ap->a); + apt[2]=sg0(exp2,p,ap->a); + apt[3]=sg2(exp2,p,ap->a); + */ + if (flags->sw[9]) { + t[8] = apt[0]*(p[50]+p[96]*plg[0][2]+p[54]*plg[0][4]+ \ + (p[125]*plg[0][1]+p[126]*plg[0][3]+p[127]*plg[0][5])*cd14*flags->swc[5]+ \ + (p[128]*plg[1][1]+p[129]*plg[1][3]+p[130]*plg[1][5])*flags->swc[7]* \ + cos(hr*(tloc-p[131]))); + } + } + } else { + double p44, p45; + apd=input->ap-4.0; + p44=p[43]; + p45=p[44]; + if (p44<0) + p44 = 1.0E-5; + apdf = apd + (p45-1.0)*(apd + (exp(-p44 * apd) - 1.0)/p44); + if (flags->sw[9]) { + t[8]=apdf*(p[32]+p[45]*plg[0][2]+p[34]*plg[0][4]+ \ + (p[100]*plg[0][1]+p[101]*plg[0][3]+p[102]*plg[0][5])*cd14*flags->swc[5]+ + (p[121]*plg[1][1]+p[122]*plg[1][3]+p[123]*plg[1][5])*flags->swc[7]* + cos(hr*(tloc-p[124]))); + } + } + + if ((flags->sw[10])&&(input->g_long>-1000.0)) { + + /* longitudinal */ + if (flags->sw[11]) { + t[10] = (1.0 + p[80]*dfa*flags->swc[1])* \ + ((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\ + +p[103]*plg[1][1]+p[104]*plg[1][3]+p[105]*plg[1][5]\ + +flags->swc[5]*(p[109]*plg[1][1]+p[110]*plg[1][3]+p[111]*plg[1][5])*cd14)* \ + cos(dgtr*input->g_long) \ + +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\ + +p[106]*plg[1][1]+p[107]*plg[1][3]+p[108]*plg[1][5]\ + +flags->swc[5]*(p[112]*plg[1][1]+p[113]*plg[1][3]+p[114]*plg[1][5])*cd14)* \ + sin(dgtr*input->g_long)); + } + + /* ut and mixed ut, longitude */ + if (flags->sw[12]){ + t[11]=(1.0+p[95]*plg[0][1])*(1.0+p[81]*dfa*flags->swc[1])*\ + (1.0+p[119]*plg[0][1]*flags->swc[5]*cd14)*\ + ((p[68]*plg[0][1]+p[69]*plg[0][3]+p[70]*plg[0][5])*\ + cos(sr*(input->sec-p[71]))); + t[11]+=flags->swc[11]*\ + (p[76]*plg[2][3]+p[77]*plg[2][5]+p[78]*plg[2][7])*\ + cos(sr*(input->sec-p[79])+2.0*dgtr*input->g_long)*(1.0+p[137]*dfa*flags->swc[1]); + } + + /* ut, longitude magnetic activity */ + if (flags->sw[13]) { + if (flags->sw[9]==-1) { + if (p[51]) { + t[12]=apt[0]*flags->swc[11]*(1.+p[132]*plg[0][1])*\ + ((p[52]*plg[1][2]+p[98]*plg[1][4]+p[67]*plg[1][6])*\ + cos(dgtr*(input->g_long-p[97])))\ + +apt[0]*flags->swc[11]*flags->swc[5]*\ + (p[133]*plg[1][1]+p[134]*plg[1][3]+p[135]*plg[1][5])*\ + cd14*cos(dgtr*(input->g_long-p[136])) \ + +apt[0]*flags->swc[12]* \ + (p[55]*plg[0][1]+p[56]*plg[0][3]+p[57]*plg[0][5])*\ + cos(sr*(input->sec-p[58])); + } + } else { + t[12] = apdf*flags->swc[11]*(1.0+p[120]*plg[0][1])*\ + ((p[60]*plg[1][2]+p[61]*plg[1][4]+p[62]*plg[1][6])*\ + cos(dgtr*(input->g_long-p[63])))\ + +apdf*flags->swc[11]*flags->swc[5]* \ + (p[115]*plg[1][1]+p[116]*plg[1][3]+p[117]*plg[1][5])* \ + cd14*cos(dgtr*(input->g_long-p[118])) \ + + apdf*flags->swc[12]* \ + (p[83]*plg[0][1]+p[84]*plg[0][3]+p[85]*plg[0][5])* \ + cos(sr*(input->sec-p[75])); + } + } + } + + /* parms not used: 82, 89, 99, 139-149 */ + tinf = p[30]; + for (i=0;i<14;i++) + tinf = tinf + abs(flags->sw[i+1])*t[i]; + return tinf; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GLOB7S ---------------------------- */ +/* ------------------------------------------------------------------- */ + +double glob7s(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) { +/* VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99 + */ + double pset=2.0; + double t[14]; + double tt; + double cd32, cd18, cd14, cd39; + double p32, p18, p14, p39; + int i,j; + double dr=1.72142E-2; + double dgtr=1.74533E-2; + /* confirm parameter set */ + if (p[99]==0) + p[99]=pset; + if (p[99]!=pset) { + printf("Wrong parameter set for glob7s\n"); + return -1; + } + for (j=0;j<14;j++) + t[j]=0.0; + cd32 = cos(dr*(input->doy-p[31])); + cd18 = cos(2.0*dr*(input->doy-p[17])); + cd14 = cos(dr*(input->doy-p[13])); + cd39 = cos(2.0*dr*(input->doy-p[38])); + p32=p[31]; + p18=p[17]; + p14=p[13]; + p39=p[38]; + + /* F10.7 */ + t[0] = p[21]*dfa; + + /* time independent */ + 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]; + + /* SYMMETRICAL ANNUAL */ + t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32; + + /* SYMMETRICAL SEMIANNUAL */ + t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18; + + /* ASYMMETRICAL ANNUAL */ + t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14; + + /* ASYMMETRICAL SEMIANNUAL */ + t[5]=(p[37]*plg[0][1])*cd39; + + /* DIURNAL */ + if (flags->sw[7]) { + double t71, t72; + t71 = p[11]*plg[1][2]*cd14*flags->swc[5]; + t72 = p[12]*plg[1][2]*cd14*flags->swc[5]; + 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) ; + } + + /* SEMIDIURNAL */ + if (flags->sw[8]) { + double t81, t82; + t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5]; + t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5]; + 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); + } + + /* TERDIURNAL */ + if (flags->sw[14]) { + t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc; + } + + /* MAGNETIC ACTIVITY */ + if (flags->sw[9]) { + if (flags->sw[9]==1) + t[8] = apdf * (p[32] + p[45] * plg[0][2] * flags->swc[2]); + if (flags->sw[9]==-1) + t[8]=(p[50]*apt[0] + p[96]*plg[0][2] * apt[0]*flags->swc[2]); + } + + /* LONGITUDINAL */ + if (!((flags->sw[10]==0) || (flags->sw[11]==0) || (input->g_long<=-1000.0))) { + t[10] = (1.0 + plg[0][1]*(p[80]*flags->swc[5]*cos(dr*(input->doy-p[81]))\ + +p[85]*flags->swc[6]*cos(2.0*dr*(input->doy-p[86])))\ + +p[83]*flags->swc[3]*cos(dr*(input->doy-p[84]))\ + +p[87]*flags->swc[4]*cos(2.0*dr*(input->doy-p[88])))\ + *((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\ + +p[74]*plg[1][1]+p[75]*plg[1][3]+p[76]*plg[1][5]\ + )*cos(dgtr*input->g_long)\ + +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\ + +p[77]*plg[1][1]+p[78]*plg[1][3]+p[79]*plg[1][5]\ + )*sin(dgtr*input->g_long)); + } + tt=0; + for (i=0;i<14;i++) + tt+=abs(flags->sw[i+1])*t[i]; + return tt; +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GTD7 ------------------------------ */ +/* ------------------------------------------------------------------- */ + +void gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) { + double xlat; + double xmm; + int mn3 = 5; + double zn3[5]={32.5,20.0,15.0,10.0,0.0}; + int mn2 = 4; + double zn2[4]={72.5,55.0,45.0,32.5}; + double altt; + double zmix=62.5; + double tmp; + double dm28m; + double tz; + double dmc; + double dmr; + double dz28; + struct nrlmsise_output soutput; + int i; + + tselec(flags); + + /* Latitude variation of gravity (none for sw[2]=0) */ + xlat=input->g_lat; + if (flags->sw[2]==0) + xlat=45.0; + glatf(xlat, &gsurf, &re); + + xmm = pdm[2][4]; + + /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */ + if (input->alt>zn2[0]) + altt=input->alt; + else + altt=zn2[0]; + + tmp=input->alt; + input->alt=altt; + gts7(input, flags, &soutput); + altt=input->alt; + input->alt=tmp; + if (flags->sw[0]) /* metric adjustment */ + dm28m=dm28*1.0E6; + else + dm28m=dm28; + output->t[0]=soutput.t[0]; + output->t[1]=soutput.t[1]; + if (input->alt>=zn2[0]) { + for (i=0;i<9;i++) + output->d[i]=soutput.d[i]; + return; + } + +/* LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0]) + * Temperature at nodes and gradients at end nodes + * Inverse temperature a linear function of spherical harmonics + */ + meso_tgn2[0]=meso_tgn1[1]; + meso_tn2[0]=meso_tn1[4]; + meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags)); + meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags)); + meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags)); + 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)); + meso_tn3[0]=meso_tn2[3]; + + if (input->alt<zn3[0]) { +/* LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0]) + * Temperature at nodes and gradients at end nodes + * Inverse temperature a linear function of spherical harmonics + */ + meso_tgn3[0]=meso_tgn2[1]; + meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags->sw[22]*glob7s(pma[3], input, flags)); + meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags)); + meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags)); + meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags)); + 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)); + } + + /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */ + + dmc=0; + if (input->alt>zmix) + dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix); + dz28=soutput.d[2]; + + /**** N2 density ****/ + dmr=soutput.d[2] / dm28m - 1.0; + output->d[2]=densm(input->alt,dm28m,xmm, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2); + output->d[2]=output->d[2] * (1.0 + dmr*dmc); + + /**** HE density ****/ + dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0; + output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc); + + /**** O density ****/ + output->d[1] = 0; + output->d[8] = 0; + + /**** O2 density ****/ + dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0; + output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc); + + /**** AR density ***/ + dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0; + output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc); + + /**** Hydrogen density ****/ + output->d[6] = 0; + + /**** Atomic nitrogen density ****/ + output->d[7] = 0; + + /**** Total mass density */ + 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]); + + if (flags->sw[0]) + output->d[5]=output->d[5]/1000; + + /**** temperature at altitude ****/ + dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2); + output->t[1]=tz; + +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GTD7D ----------------------------- */ +/* ------------------------------------------------------------------- */ + +void gtd7d(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) { + gtd7(input, flags, output); + 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]); +} + + + +/* ------------------------------------------------------------------- */ +/* -------------------------------- GHP7 ----------------------------- */ +/* ------------------------------------------------------------------- */ + +void ghp7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output, double press) { + double bm = 1.3806E-19; + double rgas = 831.4; + double test = 0.00043; + double ltest = 12; + double pl, p; + double zi; + double z; + double cl, cl2; + double ca, cd; + double xn, xm, diff; + double g, sh; + int l; + pl = log10(press); + if (pl >= -5.0) { + if (pl>2.5) + zi = 18.06 * (3.00 - pl); + else if ((pl>0.075) && (pl<=2.5)) + zi = 14.98 * (3.08 - pl); + else if ((pl>-1) && (pl<=0.075)) + zi = 17.80 * (2.72 - pl); + else if ((pl>-2) && (pl<=-1)) + zi = 14.28 * (3.64 - pl); + else if ((pl>-4) && (pl<=-2)) + zi = 12.72 * (4.32 -pl); + else if (pl<=-4) + zi = 25.3 * (0.11 - pl); + cl = input->g_lat/90.0; + cl2 = cl*cl; + if (input->doy<182) + cd = (1.0 - (double) input->doy) / 91.25; + else + cd = ((double) input->doy) / 91.25 - 3.0; + ca = 0; + if ((pl > -1.11) && (pl<=-0.23)) + ca = 1.0; + if (pl > -0.23) + ca = (2.79 - pl) / (2.79 + 0.23); + if ((pl <= -1.11) && (pl>-3)) + ca = (-2.93 - pl)/(-2.93 + 1.11); + z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl; + } else + z = 22.0 * pow((pl + 4.0),2.0) + 110.0; + + /* iteration loop */ + l = 0; + do { + l++; + input->alt = z; + gtd7(input, flags, output); + z = input->alt; + xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7]; + p = bm * xn * output->t[1]; + if (flags->sw[0]) + p = p*1.0E-6; + diff = pl - log10(p); + if (sqrt(diff*diff)<test) + return; + if (l==ltest) { + printf("ERROR: ghp7 not converging for press %e, diff %e",press,diff); + return; + } + xm = output->d[5] / xn / 1.66E-24; + if (flags->sw[0]) + xm = xm * 1.0E3; + g = gsurf / (pow((1.0 + z/re),2.0)); + sh = rgas * output->t[1] / (xm * g); + + /* new altitude estimate using scale height */ + if (l < 6) + z = z - sh * diff * 2.302; + else + z = z - sh * diff; + } while (1==1); +} + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- GTS7 ------------------------------ */ +/* ------------------------------------------------------------------- */ + +void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) { +/* Thermospheric portion of NRLMSISE-00 + * See GTD7 for more extensive comments + * alt > 72.5 km! + */ + double za; + int i, j; + double ddum, z; + double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5}; + double tinf; + int mn1 = 5; + double g0; + double tlb; + double s, z0, t0, tr12; + double db01, db04, db14, db16, db28, db32, db40, db48; + double zh28, zh04, zh16, zh32, zh40, zh01, zh14; + double zhm28, zhm04, zhm16, zhm32, zhm40, zhm01, zhm14; + double xmd; + double b28, b04, b16, b32, b40, b01, b14; + double tz; + double g28, g4, g16, g32, g40, g1, g14; + double zhf, xmm; + double zc04, zc16, zc32, zc40, zc01, zc14; + double uc04; + double hc04, hc16, hc32, hc40, hc01, hc14; + double hcc16, hcc32, hcc01, hcc14; + double zcc16, zcc32, zcc01, zcc14; + double rc16, rc32, rc01, rc14; + double rl; + double g16h, db16h, tho, zsht, zmho, zsho; + double dgtr=1.74533E-2; + double dr=1.72142E-2; + double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0}; + double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0}; + double dd; + za = pdl[1][15]; + zn1[0] = za; + for (j=0;j<9;j++) + output->d[j]=0; + + /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */ + if (input->alt>zn1[0]) + tinf = ptm[0]*pt[0] * \ + (1.0+flags->sw[16]*globe7(pt,input,flags)); + else + tinf = ptm[0]*pt[0]; + output->t[0]=tinf; + + /* GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */ + if (input->alt>zn1[4]) + g0 = ptm[3]*ps[0] * \ + (1.0+flags->sw[19]*globe7(ps,input,flags)); + else + g0 = ptm[3]*ps[0]; + tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0]; + s = g0 / (tinf - tlb); + +/* Lower thermosphere temp variations not significant for + * density above 300 km */ + if (input->alt<300.0) { + meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags)); + meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags)); + meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags)); + meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags)); + 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)); + } else { + meso_tn1[1]=ptm[6]*ptl[0][0]; + meso_tn1[2]=ptm[2]*ptl[1][0]; + meso_tn1[3]=ptm[7]*ptl[2][0]; + meso_tn1[4]=ptm[4]*ptl[3][0]; + meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(pow((ptm[4]*ptl[3][0]),2.0)); + } + + z0 = zn1[3]; + t0 = meso_tn1[3]; + tr12 = 1.0; + + /* N2 variation factor at Zlb */ + g28=flags->sw[21]*globe7(pd[2], input, flags); + + /* VARIATION OF TURBOPAUSE HEIGHT */ + zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*sin(dgtr*input->g_lat)*cos(dr*(input->doy-pt[13]))); + output->t[0]=tinf; + xmm = pdm[2][4]; + z = input->alt; + + + /**** N2 DENSITY ****/ + + /* Diffusive density at Zlb */ + db28 = pdm[2][0]*exp(g28)*pd[2][0]; + /* Diffusive density at Alt */ + output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[2]; + /* Turbopause */ + zh28=pdm[2][2]*zhf; + zhm28=pdm[2][3]*pdl[1][5]; + xmd=28.0-xmm; + /* Mixed density at Zlb */ + b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + if ((flags->sw[15])&&(z<=altl[2])) { + /* Mixed density at Alt */ + dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Net density at Alt */ + output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0); + } + + + /**** HE DENSITY ****/ + + /* Density variation factor at Zlb */ + g4 = flags->sw[21]*globe7(pd[0], input, flags); + /* Diffusive density at Zlb */ + db04 = pdm[0][0]*exp(g4)*pd[0][0]; + /* Diffusive density at Alt */ + output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[0]; + if ((flags->sw[15]) && (z<=altl[0])) { + /* Turbopause */ + zh04=pdm[0][2]; + /* Mixed density at Zlb */ + b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm04=zhm28; + /* Net density at Alt */ + output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[0][1]/b04); + zc04=pdm[0][4]*pdl[1][0]; + uc04=pdm[0][5]*pdl[1][1]; + /* Net density corrected at Alt */ + output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04); + } + + + /**** O DENSITY ****/ + + /* Density variation factor at Zlb */ + g16= flags->sw[21]*globe7(pd[1],input,flags); + /* Diffusive density at Zlb */ + db16 = pdm[1][0]*exp(g16)*pd[1][0]; + /* Diffusive density at Alt */ + output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + dd=output->d[1]; + if ((flags->sw[15]) && (z<=altl[1])) { + /* Turbopause */ + zh16=pdm[1][2]; + /* Mixed density at Zlb */ + 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); + /* Mixed density at Alt */ + dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm16=zhm28; + /* Net density at Alt */ + output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.); + rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0)); + hc16=pdm[1][5]*pdl[1][3]; + zc16=pdm[1][4]*pdl[1][2]; + output->d[1]=output->d[1]*ccor(z,rl,hc16,zc16); + /* Chemistry correction */ + hcc16=pdm[1][7]*pdl[1][13]; + zcc16=pdm[1][6]*pdl[1][12]; + rc16=pdm[1][3]*pdl[1][14]; + /* Net density corrected at Alt */ + output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16); + } + + + /**** O2 DENSITY ****/ + + /* Density variation factor at Zlb */ + g32= flags->sw[21]*globe7(pd[4], input, flags); + /* Diffusive density at Zlb */ + db32 = pdm[3][0]*exp(g32)*pd[4][0]; + /* Diffusive density at Alt */ + output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + dd=output->d[3]; + if (flags->sw[15]) { + if (z<=altl[3]) { + /* Turbopause */ + zh32=pdm[3][2]; + /* Mixed density at Zlb */ + b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm32=zhm28; + /* Net density at Alt */ + output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[3][1]/b32); + hc32=pdm[3][5]*pdl[1][7]; + zc32=pdm[3][4]*pdl[1][6]; + output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32); + } + /* Correction for general departure from diffusive equilibrium above Zlb */ + hcc32=pdm[3][7]*pdl[1][22]; + zcc32=pdm[3][6]*pdl[1][21]; + rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.)); + /* Net density corrected at Alt */ + output->d[3]=output->d[3]*ccor(z,rc32,hcc32,zcc32); + } + + + /**** AR DENSITY ****/ + + /* Density variation factor at Zlb */ + g40= flags->sw[20]*globe7(pd[5],input,flags); + /* Diffusive density at Zlb */ + db40 = pdm[4][0]*exp(g40)*pd[5][0]; + /* Diffusive density at Alt */ + output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[4]; + if ((flags->sw[15]) && (z<=altl[4])) { + /* Turbopause */ + zh40=pdm[4][2]; + /* Mixed density at Zlb */ + b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm40=zhm28; + /* Net density at Alt */ + output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[4][1]/b40); + hc40=pdm[4][5]*pdl[1][9]; + zc40=pdm[4][4]*pdl[1][8]; + /* Net density corrected at Alt */ + output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40); + } + + + /**** HYDROGEN DENSITY ****/ + + /* Density variation factor at Zlb */ + g1 = flags->sw[21]*globe7(pd[6], input, flags); + /* Diffusive density at Zlb */ + db01 = pdm[5][0]*exp(g1)*pd[6][0]; + /* Diffusive density at Alt */ + output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[6]; + if ((flags->sw[15]) && (z<=altl[6])) { + /* Turbopause */ + zh01=pdm[5][2]; + /* Mixed density at Zlb */ + b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm01=zhm28; + /* Net density at Alt */ + output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[5][1]*sqrt(pdl[1][17]*pdl[1][17])/b01); + hc01=pdm[5][5]*pdl[1][11]; + zc01=pdm[5][4]*pdl[1][10]; + output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01); + /* Chemistry correction */ + hcc01=pdm[5][7]*pdl[1][19]; + zcc01=pdm[5][6]*pdl[1][18]; + rc01=pdm[5][3]*pdl[1][20]; + /* Net density corrected at Alt */ + output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01); +} + + + /**** ATOMIC NITROGEN DENSITY ****/ + + /* Density variation factor at Zlb */ + g14 = flags->sw[21]*globe7(pd[7],input,flags); + /* Diffusive density at Zlb */ + db14 = pdm[6][0]*exp(g14)*pd[7][0]; + /* Diffusive density at Alt */ + output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + dd=output->d[7]; + if ((flags->sw[15]) && (z<=altl[7])) { + /* Turbopause */ + zh14=pdm[6][2]; + /* Mixed density at Zlb */ + b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + /* Mixed density at Alt */ + dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1); + zhm14=zhm28; + /* Net density at Alt */ + output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.); + /* Correction to specified mixing ratio at ground */ + rl=log(b28*pdm[6][1]*sqrt(pdl[0][2]*pdl[0][2])/b14); + hc14=pdm[6][5]*pdl[0][1]; + zc14=pdm[6][4]*pdl[0][0]; + output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14); + /* Chemistry correction */ + hcc14=pdm[6][7]*pdl[0][4]; + zcc14=pdm[6][6]*pdl[0][3]; + rc14=pdm[6][3]*pdl[0][5]; + /* Net density corrected at Alt */ + output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14); + } + + + /**** Anomalous OXYGEN DENSITY ****/ + + g16h = flags->sw[21]*globe7(pd[8],input,flags); + db16h = pdm[7][0]*exp(g16h)*pd[8][0]; + tho = pdm[7][9]*pdl[0][6]; + dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1); + zsht=pdm[7][5]; + zmho=pdm[7][4]; + zsho=scalh(zmho,16.0,tho); + output->d[8]=dd*exp(-zsht/zsho*(exp(-(z-zmho)/zsht)-1.)); + + + /* total mass density */ + 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]); + db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14); + + + + /* temperature */ + z = sqrt(input->alt*input->alt); + ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1); + if (flags->sw[0]) { + for(i=0;i<9;i++) + output->d[i]=output->d[i]*1.0E6; + output->d[5]=output->d[5]/1000; + } +} diff --git a/nrlmsise-00.h b/nrlmsise-00.h new file mode 100644 index 0000000..842a4b0 --- /dev/null +++ b/nrlmsise-00.h @@ -0,0 +1,222 @@ +/* -------------------------------------------------------------------- */ +/* --------- N R L M S I S E - 0 0 M O D E L 2 0 0 1 ---------- */ +/* -------------------------------------------------------------------- */ + +/* This file is part of the NRLMSISE-00 C source code package - release + * 20020305 + * + * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and + * Doug Drob. They also wrote a NRLMSISE-00 distribution package in + * FORTRAN which is available at + * http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm + * + * Dominik Brodowski implemented and maintains this C version. You can + * reach him at devel@brodo.de. See the file "DOCUMENTATION" for details, + * and check http://www.brodo.de/english/pub/nrlmsise/index.html for + * updated releases of this package. + */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------- INPUT ----------------------------- */ +/* ------------------------------------------------------------------- */ + +struct nrlmsise_flags { + int switches[24]; + double sw[24]; + double swc[24]; +}; +/* + * Switches: to turn on and off particular variations use these switches. + * 0 is off, 1 is on, and 2 is main effects off but cross terms on. + * + * Standard values are 0 for switch 0 and 1 for switches 1 to 23. The + * array "switches" needs to be set accordingly by the calling program. + * The arrays sw and swc are set internally. + * + * switches[i]: + * i - explanation + * ----------------- + * 0 - output in centimeters instead of meters + * 1 - F10.7 effect on mean + * 2 - time independent + * 3 - symmetrical annual + * 4 - symmetrical semiannual + * 5 - asymmetrical annual + * 6 - asymmetrical semiannual + * 7 - diurnal + * 8 - semidiurnal + * 9 - daily ap [when this is set to -1 (!) the pointer + * ap_a in struct nrlmsise_input must + * point to a struct ap_array] + * 10 - all UT/long effects + * 11 - longitudinal + * 12 - UT and mixed UT/long + * 13 - mixed AP/UT/LONG + * 14 - terdiurnal + * 15 - departures from diffusive equilibrium + * 16 - all TINF var + * 17 - all TLB var + * 18 - all TN1 var + * 19 - all S var + * 20 - all TN2 var + * 21 - all NLB var + * 22 - all TN3 var + * 23 - turbo scale height var + */ + +struct ap_array { + double a[7]; +}; +/* Array containing the following magnetic values: + * 0 : daily AP + * 1 : 3 hr AP index for current time + * 2 : 3 hr AP index for 3 hrs before current time + * 3 : 3 hr AP index for 6 hrs before current time + * 4 : 3 hr AP index for 9 hrs before current time + * 5 : Average of eight 3 hr AP indicies from 12 to 33 hrs + * prior to current time + * 6 : Average of eight 3 hr AP indicies from 36 to 57 hrs + * prior to current time + */ + + +struct nrlmsise_input { + int year; /* year, currently ignored */ + int doy; /* day of year */ + double sec; /* seconds in day (UT) */ + double alt; /* altitude in kilometes */ + double g_lat; /* geodetic latitude */ + double g_long; /* geodetic longitude */ + double lst; /* local apparent solar time (hours), see note below */ + double f107A; /* 81 day average of F10.7 flux (centered on doy) */ + double f107; /* daily F10.7 flux for previous day */ + double ap; /* magnetic index(daily) */ + struct ap_array *ap_a; /* see above */ +}; +/* + * NOTES ON INPUT VARIABLES: + * UT, Local Time, and Longitude are used independently in the + * model and are not of equal importance for every situation. + * For the most physically realistic calculation these three + * variables should be consistent (lst=sec/3600 + g_long/15). + * The Equation of Time departures from the above formula + * for apparent local time can be included if available but + * are of minor importance. + * + * f107 and f107A values used to generate the model correspond + * to the 10.7 cm radio flux at the actual distance of the Earth + * from the Sun rather than the radio flux at 1 AU. The following + * site provides both classes of values: + * ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/ + * + * f107, f107A, and ap effects are neither large nor well + * established below 80 km and these parameters should be set to + * 150., 150., and 4. respectively. + */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ OUTPUT ----------------------------- */ +/* ------------------------------------------------------------------- */ + +struct nrlmsise_output { + double d[9]; /* densities */ + double t[2]; /* temperatures */ +}; +/* + * OUTPUT VARIABLES: + * d[0] - HE NUMBER DENSITY(CM-3) + * d[1] - O NUMBER DENSITY(CM-3) + * d[2] - N2 NUMBER DENSITY(CM-3) + * d[3] - O2 NUMBER DENSITY(CM-3) + * d[4] - AR NUMBER DENSITY(CM-3) + * d[5] - TOTAL MASS DENSITY(GM/CM3) [includes d[8] in td7d] + * d[6] - H NUMBER DENSITY(CM-3) + * d[7] - N NUMBER DENSITY(CM-3) + * d[8] - Anomalous oxygen NUMBER DENSITY(CM-3) + * t[0] - EXOSPHERIC TEMPERATURE + * t[1] - TEMPERATURE AT ALT + * + * + * O, H, and N are set to zero below 72.5 km + * + * t[0], Exospheric temperature, is set to global average for + * altitudes below 120 km. The 120 km gradient is left at global + * average value for altitudes below 72 km. + * + * d[5], TOTAL MASS DENSITY, is NOT the same for subroutines GTD7 + * and GTD7D + * + * SUBROUTINE GTD7 -- d[5] is the sum of the mass densities of the + * species labeled by indices 0-4 and 6-7 in output variable d. + * This includes He, O, N2, O2, Ar, H, and N but does NOT include + * anomalous oxygen (species index 8). + * + * SUBROUTINE GTD7D -- d[5] is the "effective total mass density + * for drag" and is the sum of the mass densities of all species + * in this model, INCLUDING anomalous oxygen. + */ + + + +/* ------------------------------------------------------------------- */ +/* --------------------------- PROTOTYPES ---------------------------- */ +/* ------------------------------------------------------------------- */ + +/* GTD7 */ +/* Neutral Atmosphere Empircial Model from the surface to lower + * exosphere. + */ +void gtd7 (struct nrlmsise_input *input, \ + struct nrlmsise_flags *flags, \ + struct nrlmsise_output *output); + + +/* GTD7D */ +/* This subroutine provides Effective Total Mass Density for output + * d[5] which includes contributions from "anomalous oxygen" which can + * affect satellite drag above 500 km. See the section "output" for + * additional details. + */ +void gtd7d(struct nrlmsise_input *input, \ + struct nrlmsise_flags *flags, \ + struct nrlmsise_output *output); + + +/* GTS7 */ +/* Thermospheric portion of NRLMSISE-00 + */ +void gts7 (struct nrlmsise_input *input, \ + struct nrlmsise_flags *flags, \ + struct nrlmsise_output *output); + + +/* GHP7 */ +/* To specify outputs at a pressure level (press) rather than at + * an altitude. + */ +void ghp7 (struct nrlmsise_input *input, \ + struct nrlmsise_flags *flags, \ + struct nrlmsise_output *output, \ + double press); + + + +/* ------------------------------------------------------------------- */ +/* ----------------------- COMPILATION TWEAKS ------------------------ */ +/* ------------------------------------------------------------------- */ + +/* "inlining" of functions */ +/* Some compilers (e.g. gcc) allow the inlining of functions into the + * calling routine. This means a lot of overhead can be removed, and + * the execution of the program runs much faster. However, the filesize + * and thus the loading time is increased. + */ +#ifdef INLINE +#define __inline_double inline double +#else +#define __inline_double double +#endif diff --git a/nrlmsise-00_data.c b/nrlmsise-00_data.c new file mode 100644 index 0000000..610cd16 --- /dev/null +++ b/nrlmsise-00_data.c @@ -0,0 +1,740 @@ +/* -------------------------------------------------------------------- */ +/* --------- N R L M S I S E - 0 0 M O D E L 2 0 0 1 ---------- */ +/* -------------------------------------------------------------------- */ + +/* This file is part of the NRLMSISE-00 C source code package - release + * 20020305 + * + * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and + * Doug Drob. They also wrote a NRLMSISE-00 distribution package in + * FORTRAN which is available at + * http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm + * + * Dominik Brodowski implemented and maintains this C version. You can + * reach him at devel@brodo.de. See the file "DOCUMENTATION" for details, + * and check http://www.brodo.de/english/pub/nrlmsise/index.html for + * updated releases of this package. + */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------ BLOCK DATA GTD7BK ------------------------ */ +/* ------------------------------------------------------------------- */ + +/* TEMPERATURE */ +double pt[150] = { + 9.86573E-01, 1.62228E-02, 1.55270E-02,-1.04323E-01,-3.75801E-03, + -1.18538E-03,-1.24043E-01, 4.56820E-03, 8.76018E-03,-1.36235E-01, + -3.52427E-02, 8.84181E-03,-5.92127E-03,-8.61650E+00, 0.00000E+00, + 1.28492E-02, 0.00000E+00, 1.30096E+02, 1.04567E-02, 1.65686E-03, + -5.53887E-06, 2.97810E-03, 0.00000E+00, 5.13122E-03, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.27026E-06, + 0.00000E+00, 6.74494E+00, 4.93933E-03, 2.21656E-03, 2.50802E-03, + 0.00000E+00, 0.00000E+00,-2.08841E-02,-1.79873E+00, 1.45103E-03, + 2.81769E-04,-1.44703E-03,-5.16394E-05, 8.47001E-02, 1.70147E-01, + 5.72562E-03, 5.07493E-05, 4.36148E-03, 1.17863E-04, 4.74364E-03, + 6.61278E-03, 4.34292E-05, 1.44373E-03, 2.41470E-05, 2.84426E-03, + 8.56560E-04, 2.04028E-03, 0.00000E+00,-3.15994E+03,-2.46423E-03, + 1.13843E-03, 4.20512E-04, 0.00000E+00,-9.77214E+01, 6.77794E-03, + 5.27499E-03, 1.14936E-03, 0.00000E+00,-6.61311E-03,-1.84255E-02, + -1.96259E-02, 2.98618E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 6.44574E+02, 8.84668E-04, 5.05066E-04, 0.00000E+00, 4.02881E+03, + -1.89503E-03, 0.00000E+00, 0.00000E+00, 8.21407E-04, 2.06780E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -1.20410E-02,-3.63963E-03, 9.92070E-05,-1.15284E-04,-6.33059E-05, + -6.05545E-01, 8.34218E-03,-9.13036E+01, 3.71042E-04, 0.00000E+00, + 4.19000E-04, 2.70928E-03, 3.31507E-03,-4.44508E-03,-4.96334E-03, + -1.60449E-03, 3.95119E-03, 2.48924E-03, 5.09815E-04, 4.05302E-03, + 2.24076E-03, 0.00000E+00, 6.84256E-03, 4.66354E-04, 0.00000E+00, + -3.68328E-04, 0.00000E+00, 0.00000E+00,-1.46870E+02, 0.00000E+00, + 0.00000E+00, 1.09501E-03, 4.65156E-04, 5.62583E-04, 3.21596E+00, + 6.43168E-04, 3.14860E-03, 3.40738E-03, 1.78481E-03, 9.62532E-04, + 5.58171E-04, 3.43731E+00,-2.33195E-01, 5.10289E-04, 0.00000E+00, + 0.00000E+00,-9.25347E+04, 0.00000E+00,-1.99639E-03, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}; + +double pd[9][150] = { +/* HE DENSITY */ { + 1.09979E+00,-4.88060E-02,-1.97501E-01,-9.10280E-02,-6.96558E-03, + 2.42136E-02, 3.91333E-01,-7.20068E-03,-3.22718E-02, 1.41508E+00, + 1.68194E-01, 1.85282E-02, 1.09384E-01,-7.24282E+00, 0.00000E+00, + 2.96377E-01,-4.97210E-02, 1.04114E+02,-8.61108E-02,-7.29177E-04, + 1.48998E-06, 1.08629E-03, 0.00000E+00, 0.00000E+00, 8.31090E-02, + 1.12818E-01,-5.75005E-02,-1.29919E-02,-1.78849E-02,-2.86343E-06, + 0.00000E+00,-1.51187E+02,-6.65902E-03, 0.00000E+00,-2.02069E-03, + 0.00000E+00, 0.00000E+00, 4.32264E-02,-2.80444E+01,-3.26789E-03, + 2.47461E-03, 0.00000E+00, 0.00000E+00, 9.82100E-02, 1.22714E-01, + -3.96450E-02, 0.00000E+00,-2.76489E-03, 0.00000E+00, 1.87723E-03, + -8.09813E-03, 4.34428E-05,-7.70932E-03, 0.00000E+00,-2.28894E-03, + -5.69070E-03,-5.22193E-03, 6.00692E-03,-7.80434E+03,-3.48336E-03, + -6.38362E-03,-1.82190E-03, 0.00000E+00,-7.58976E+01,-2.17875E-02, + -1.72524E-02,-9.06287E-03, 0.00000E+00, 2.44725E-02, 8.66040E-02, + 1.05712E-01, 3.02543E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -6.01364E+03,-5.64668E-03,-2.54157E-03, 0.00000E+00, 3.15611E+02, + -5.69158E-03, 0.00000E+00, 0.00000E+00,-4.47216E-03,-4.49523E-03, + 4.64428E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 4.51236E-02, 2.46520E-02, 6.17794E-03, 0.00000E+00, 0.00000E+00, + -3.62944E-01,-4.80022E-02,-7.57230E+01,-1.99656E-03, 0.00000E+00, + -5.18780E-03,-1.73990E-02,-9.03485E-03, 7.48465E-03, 1.53267E-02, + 1.06296E-02, 1.18655E-02, 2.55569E-03, 1.69020E-03, 3.51936E-02, + -1.81242E-02, 0.00000E+00,-1.00529E-01,-5.10574E-03, 0.00000E+00, + 2.10228E-03, 0.00000E+00, 0.00000E+00,-1.73255E+02, 5.07833E-01, + -2.41408E-01, 8.75414E-03, 2.77527E-03,-8.90353E-05,-5.25148E+00, + -5.83899E-03,-2.09122E-02,-9.63530E-03, 9.77164E-03, 4.07051E-03, + 2.53555E-04,-5.52875E+00,-3.55993E-01,-2.49231E-03, 0.00000E+00, + 0.00000E+00, 2.86026E+01, 0.00000E+00, 3.42722E-04, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* O DENSITY */ { + 1.02315E+00,-1.59710E-01,-1.06630E-01,-1.77074E-02,-4.42726E-03, + 3.44803E-02, 4.45613E-02,-3.33751E-02,-5.73598E-02, 3.50360E-01, + 6.33053E-02, 2.16221E-02, 5.42577E-02,-5.74193E+00, 0.00000E+00, + 1.90891E-01,-1.39194E-02, 1.01102E+02, 8.16363E-02, 1.33717E-04, + 6.54403E-06, 3.10295E-03, 0.00000E+00, 0.00000E+00, 5.38205E-02, + 1.23910E-01,-1.39831E-02, 0.00000E+00, 0.00000E+00,-3.95915E-06, + 0.00000E+00,-7.14651E-01,-5.01027E-03, 0.00000E+00,-3.24756E-03, + 0.00000E+00, 0.00000E+00, 4.42173E-02,-1.31598E+01,-3.15626E-03, + 1.24574E-03,-1.47626E-03,-1.55461E-03, 6.40682E-02, 1.34898E-01, + -2.42415E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 6.13666E-04, + -5.40373E-03, 2.61635E-05,-3.33012E-03, 0.00000E+00,-3.08101E-03, + -2.42679E-03,-3.36086E-03, 0.00000E+00,-1.18979E+03,-5.04738E-02, + -2.61547E-03,-1.03132E-03, 1.91583E-04,-8.38132E+01,-1.40517E-02, + -1.14167E-02,-4.08012E-03, 1.73522E-04,-1.39644E-02,-6.64128E-02, + -6.85152E-02,-1.34414E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 6.07916E+02,-4.12220E-03,-2.20996E-03, 0.00000E+00, 1.70277E+03, + -4.63015E-03, 0.00000E+00, 0.00000E+00,-2.25360E-03,-2.96204E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 3.92786E-02, 1.31186E-02,-1.78086E-03, 0.00000E+00, 0.00000E+00, + -3.90083E-01,-2.84741E-02,-7.78400E+01,-1.02601E-03, 0.00000E+00, + -7.26485E-04,-5.42181E-03,-5.59305E-03, 1.22825E-02, 1.23868E-02, + 6.68835E-03,-1.03303E-02,-9.51903E-03, 2.70021E-04,-2.57084E-02, + -1.32430E-02, 0.00000E+00,-3.81000E-02,-3.16810E-03, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-9.05762E-04,-2.14590E-03,-1.17824E-03, 3.66732E+00, + -3.79729E-04,-6.13966E-03,-5.09082E-03,-1.96332E-03,-3.08280E-03, + -9.75222E-04, 4.03315E+00,-2.52710E-01, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* N2 DENSITY */ { + 1.16112E+00, 0.00000E+00, 0.00000E+00, 3.33725E-02, 0.00000E+00, + 3.48637E-02,-5.44368E-03, 0.00000E+00,-6.73940E-02, 1.74754E-01, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00, + 1.26733E-01, 0.00000E+00, 1.03154E+02, 5.52075E-02, 0.00000E+00, + 0.00000E+00, 8.13525E-04, 0.00000E+00, 0.00000E+00, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-2.50482E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.48894E-03, + 6.16053E-04,-5.79716E-04, 2.95482E-03, 8.47001E-02, 1.70147E-01, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* TLB */ { + 9.44846E-01, 0.00000E+00, 0.00000E+00,-3.08617E-02, 0.00000E+00, + -2.44019E-02, 6.48607E-03, 0.00000E+00, 3.08181E-02, 4.59392E-02, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00, + 2.13260E-02, 0.00000E+00,-3.56958E+02, 0.00000E+00, 1.82278E-04, + 0.00000E+00, 3.07472E-04, 0.00000E+00, 0.00000E+00, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 3.83054E-03, 0.00000E+00, 0.00000E+00, + -1.93065E-03,-1.45090E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-1.23493E-03, 1.36736E-03, 8.47001E-02, 1.70147E-01, + 3.71469E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 5.10250E-03, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 3.68756E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* O2 DENSITY */ { + 1.38720E+00, 1.44816E-01, 0.00000E+00, 6.07767E-02, 0.00000E+00, + 2.94777E-02, 7.46900E-02, 0.00000E+00,-9.23822E-02, 8.57342E-02, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.38636E+01, 0.00000E+00, + 7.71653E-02, 0.00000E+00, 8.18751E+01, 1.87736E-02, 0.00000E+00, + 0.00000E+00, 1.49667E-02, 0.00000E+00, 0.00000E+00, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-3.67874E+02, 5.48158E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01, + 1.22631E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 8.17187E-03, 3.71617E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.10826E-03, + -3.13640E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -7.35742E-02,-5.00266E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 1.94965E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* AR DENSITY */ { + 1.04761E+00, 2.00165E-01, 2.37697E-01, 3.68552E-02, 0.00000E+00, + 3.57202E-02,-2.14075E-01, 0.00000E+00,-1.08018E-01,-3.73981E-01, + 0.00000E+00, 3.10022E-02,-1.16305E-03,-2.07596E+01, 0.00000E+00, + 8.64502E-02, 0.00000E+00, 9.74908E+01, 5.16707E-02, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 3.46193E+02, 1.34297E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.48509E-03, + -1.54689E-04, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01, + 1.47753E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 1.89320E-02, 3.68181E-05, 1.32570E-02, 0.00000E+00, 0.00000E+00, + 3.59719E-03, 7.44328E-03,-1.00023E-03,-6.50528E+03, 0.00000E+00, + 1.03485E-02,-1.00983E-03,-4.06916E-03,-6.60864E+01,-1.71533E-02, + 1.10605E-02, 1.20300E-02,-5.20034E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -2.62769E+03, 7.13755E-03, 4.17999E-03, 0.00000E+00, 1.25910E+04, + 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.23595E-03, 4.60217E-03, + 5.71794E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -3.18353E-02,-2.35526E-02,-1.36189E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.03522E-02,-6.67837E+01,-1.09724E-03, 0.00000E+00, + -1.38821E-02, 1.60468E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.51574E-02, + -5.44470E-04, 0.00000E+00, 7.28224E-02, 6.59413E-02, 0.00000E+00, + -5.15692E-03, 0.00000E+00, 0.00000E+00,-3.70367E+03, 0.00000E+00, + 0.00000E+00, 1.36131E-02, 5.38153E-03, 0.00000E+00, 4.76285E+00, + -1.75677E-02, 2.26301E-02, 0.00000E+00, 1.76631E-02, 4.77162E-03, + 0.00000E+00, 5.39354E+00, 0.00000E+00,-7.51710E-03, 0.00000E+00, + 0.00000E+00,-8.82736E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* H DENSITY */ { + 1.26376E+00,-2.14304E-01,-1.49984E-01, 2.30404E-01, 2.98237E-02, + 2.68673E-02, 2.96228E-01, 2.21900E-02,-2.07655E-02, 4.52506E-01, + 1.20105E-01, 3.24420E-02, 4.24816E-02,-9.14313E+00, 0.00000E+00, + 2.47178E-02,-2.88229E-02, 8.12805E+01, 5.10380E-02,-5.80611E-03, + 2.51236E-05,-1.24083E-02, 0.00000E+00, 0.00000E+00, 8.66784E-02, + 1.58727E-01,-3.48190E-02, 0.00000E+00, 0.00000E+00, 2.89885E-05, + 0.00000E+00, 1.53595E+02,-1.68604E-02, 0.00000E+00, 1.01015E-02, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.84552E-04, + -1.22181E-03, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01, + -1.04927E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,-5.91313E-03, + -2.30501E-02, 3.14758E-05, 0.00000E+00, 0.00000E+00, 1.26956E-02, + 8.35489E-03, 3.10513E-04, 0.00000E+00, 3.42119E+03,-2.45017E-03, + -4.27154E-04, 5.45152E-04, 1.89896E-03, 2.89121E+01,-6.49973E-03, + -1.93855E-02,-1.48492E-02, 0.00000E+00,-5.10576E-02, 7.87306E-02, + 9.51981E-02,-1.49422E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 2.65503E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 6.37110E-03, 3.24789E-04, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 6.14274E-02, 1.00376E-02,-8.41083E-04, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-1.27099E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -3.94077E-03,-1.28601E-02,-7.97616E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-6.71465E-03,-1.69799E-03, 1.93772E-03, 3.81140E+00, + -7.79290E-03,-1.82589E-02,-1.25860E-02,-1.04311E-02,-3.02465E-03, + 2.43063E-03, 3.63237E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* N DENSITY */ { + 7.09557E+01,-3.26740E-01, 0.00000E+00,-5.16829E-01,-1.71664E-03, + 9.09310E-02,-6.71500E-01,-1.47771E-01,-9.27471E-02,-2.30862E-01, + -1.56410E-01, 1.34455E-02,-1.19717E-01, 2.52151E+00, 0.00000E+00, + -2.41582E-01, 5.92939E-02, 4.39756E+00, 9.15280E-02, 4.41292E-03, + 0.00000E+00, 8.66807E-03, 0.00000E+00, 0.00000E+00, 8.66784E-02, + 1.58727E-01, 9.74701E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 6.70217E+01,-1.31660E-03, 0.00000E+00,-1.65317E-02, + 0.00000E+00, 0.00000E+00, 8.50247E-02, 2.77428E+01, 4.98658E-03, + 6.15115E-03, 9.50156E-03,-2.12723E-02, 8.47001E-02, 1.70147E-01, + -2.38645E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.37380E-03, + -8.41918E-03, 2.80145E-05, 7.12383E-03, 0.00000E+00,-1.66209E-02, + 1.03533E-04,-1.68898E-02, 0.00000E+00, 3.64526E+03, 0.00000E+00, + 6.54077E-03, 3.69130E-04, 9.94419E-04, 8.42803E+01,-1.16124E-02, + -7.74414E-03,-1.68844E-03, 1.42809E-03,-1.92955E-03, 1.17225E-01, + -2.41512E-02, 1.50521E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 1.60261E+03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.54403E-04,-1.87270E-02, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 2.76439E-02, 6.43207E-03,-3.54300E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-2.80221E-02, 8.11228E+01,-6.75255E-04, 0.00000E+00, + -1.05162E-02,-3.48292E-03,-6.97321E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-1.45546E-03,-1.31970E-02,-3.57751E-03,-1.09021E+00, + -1.50181E-02,-7.12841E-03,-6.64590E-03,-3.52610E-03,-1.87773E-02, + -2.22432E-03,-3.93895E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}, /* HOT O DENSITY */ { + 6.04050E-02, 1.57034E+00, 2.99387E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.51018E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00,-8.61650E+00, 1.26454E-02, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 5.50878E-03, 0.00000E+00, 0.00000E+00, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 6.23881E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01, + -9.45934E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}}; +/* S PARAM */ +double ps[150] = { + 9.56827E-01, 6.20637E-02, 3.18433E-02, 0.00000E+00, 0.00000E+00, + 3.94900E-02, 0.00000E+00, 0.00000E+00,-9.24882E-03,-7.94023E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.74677E-03, 0.00000E+00, 1.54951E-02, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00,-6.99007E-04, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 1.24362E-02,-5.28756E-03, 8.47001E-02, 1.70147E-01, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}; + +/* TURBO */ +double pdl[2][25] = { + { 1.09930E+00, 3.90631E+00, 3.07165E+00, 9.86161E-01, 1.63536E+01, + 4.63830E+00, 1.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.93318E-02, 1.18339E-01 }, + { 1.22732E+00, 1.02669E-01, 1.17681E+00, 2.12185E+00, 1.00000E+00, + 1.00000E+00, 1.08607E+00, 1.34836E+00, 1.10016E+00, 7.34129E-01, + 1.15241E+00, 2.22784E+00, 7.95907E-01, 4.03601E+00, 4.39732E+00, + 1.23435E+02,-4.52411E-02, 1.68986E-06, 7.44294E-01, 1.03604E+00, + 1.72783E+02, 1.17681E+00, 2.12185E+00,-7.83697E-01, 9.49154E-01 } +}; +/* LOWER BOUNDARY */ +double ptm[50] = { + 1.04130E+03, 3.86000E+02, 1.95000E+02, 1.66728E+01, 2.13000E+02, + 1.20000E+02, 2.40000E+02, 1.87000E+02,-2.00000E+00, 0.00000E+00 +}; +double pdm[8][10] = { +{ 2.45600E+07, 6.71072E-06, 1.00000E+02, 0.00000E+00, 1.10000E+02, + 1.00000E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 },\ +{ 8.59400E+10, 1.00000E+00, 1.05000E+02,-8.00000E+00, 1.10000E+02, + 1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00 },\ +{ 2.81000E+11, 0.00000E+00, 1.05000E+02, 2.80000E+01, 2.89500E+01, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 }, +{ 3.30000E+10, 2.68270E-01, 1.05000E+02, 1.00000E+00, 1.10000E+02, + 1.00000E+01, 1.10000E+02,-1.00000E+01, 0.00000E+00, 0.00000E+00 }, +{ 1.33000E+09, 1.19615E-02, 1.05000E+02, 0.00000E+00, 1.10000E+02, + 1.00000E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 }, +{ 1.76100E+05, 1.00000E+00, 9.50000E+01,-8.00000E+00, 1.10000E+02, + 1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00, }, +{ 1.00000E+07, 1.00000E+00, 1.05000E+02,-8.00000E+00, 1.10000E+02, + 1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00 }, +{ 1.00000E+06, 1.00000E+00, 1.05000E+02,-8.00000E+00, 5.50000E+02, + 7.60000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 4.00000E+03 }}; + + +double ptl[4][100] = { +/* TN1(2) */ { + 1.00858E+00, 4.56011E-02,-2.22972E-02,-5.44388E-02, 5.23136E-04, + -1.88849E-02, 5.23707E-02,-9.43646E-03, 6.31707E-03,-7.80460E-02, + -4.88430E-02, 0.00000E+00, 0.00000E+00,-7.60250E+00, 0.00000E+00, + -1.44635E-02,-1.76843E-02,-1.21517E+02, 2.85647E-02, 0.00000E+00, + 0.00000E+00, 6.31792E-04, 0.00000E+00, 5.77197E-03, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-8.90272E+03, 3.30611E-03, 3.02172E-03, 0.00000E+00, + -2.13673E-03,-3.20910E-04, 0.00000E+00, 0.00000E+00, 2.76034E-03, + 2.82487E-03,-2.97592E-04,-4.21534E-03, 8.47001E-02, 1.70147E-01, + 8.96456E-03, 0.00000E+00,-1.08596E-02, 0.00000E+00, 0.00000E+00, + 5.57917E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 9.65405E-03, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN1(3) */ { + 9.39664E-01, 8.56514E-02,-6.79989E-03, 2.65929E-02,-4.74283E-03, + 1.21855E-02,-2.14905E-02, 6.49651E-03,-2.05477E-02,-4.24952E-02, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.19148E+01, 0.00000E+00, + 1.18777E-02,-7.28230E-02,-8.15965E+01, 1.73887E-02, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-1.44691E-02, 2.80259E-04, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.16584E+02, 3.18713E-03, 7.37479E-03, 0.00000E+00, + -2.55018E-03,-3.92806E-03, 0.00000E+00, 0.00000E+00,-2.89757E-03, + -1.33549E-03, 1.02661E-03, 3.53775E-04, 8.47001E-02, 1.70147E-01, + -9.17497E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 3.56082E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-1.00902E-02, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN1(4) */ { + 9.85982E-01,-4.55435E-02, 1.21106E-02, 2.04127E-02,-2.40836E-03, + 1.11383E-02,-4.51926E-02, 1.35074E-02,-6.54139E-03, 1.15275E-01, + 1.28247E-01, 0.00000E+00, 0.00000E+00,-5.30705E+00, 0.00000E+00, + -3.79332E-02,-6.24741E-02, 7.71062E-01, 2.96315E-02, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 6.81051E-03,-4.34767E-03, 8.66784E-02, + 1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 1.07003E+01,-2.76907E-03, 4.32474E-04, 0.00000E+00, + 1.31497E-03,-6.47517E-04, 0.00000E+00,-2.20621E+01,-1.10804E-03, + -8.09338E-04, 4.18184E-04, 4.29650E-03, 8.47001E-02, 1.70147E-01, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -4.04337E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-9.52550E-04, + 8.56253E-04, 4.33114E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.21223E-03, + 2.38694E-04, 9.15245E-04, 1.28385E-03, 8.67668E-04,-5.61425E-06, + 1.04445E+00, 3.41112E+01, 0.00000E+00,-8.40704E-01,-2.39639E+02, + 7.06668E-01,-2.05873E+01,-3.63696E-01, 2.39245E+01, 0.00000E+00, + -1.06657E-03,-7.67292E-04, 1.54534E-04, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN1(5) TN2(1) */ { + 1.00320E+00, 3.83501E-02,-2.38983E-03, 2.83950E-03, 4.20956E-03, + 5.86619E-04, 2.19054E-02,-1.00946E-02,-3.50259E-03, 4.17392E-02, + -8.44404E-03, 0.00000E+00, 0.00000E+00, 4.96949E+00, 0.00000E+00, + -7.06478E-03,-1.46494E-02, 3.13258E+01,-1.86493E-03, 0.00000E+00, + -1.67499E-02, 0.00000E+00, 0.00000E+00, 5.12686E-04, 8.66784E-02, + 1.58727E-01,-4.64167E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 4.37353E-03,-1.99069E+02, 0.00000E+00,-5.34884E-03, 0.00000E+00, + 1.62458E-03, 2.93016E-03, 2.67926E-03, 5.90449E+02, 0.00000E+00, + 0.00000E+00,-1.17266E-03,-3.58890E-04, 8.47001E-02, 1.70147E-01, + 0.00000E+00, 0.00000E+00, 1.38673E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.60571E-03, + 6.28078E-04, 5.05469E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.57829E-03, + -4.00855E-04, 5.04077E-05,-1.39001E-03,-2.33406E-03,-4.81197E-04, + 1.46758E+00, 6.20332E+00, 0.00000E+00, 3.66476E-01,-6.19760E+01, + 3.09198E-01,-1.98999E+01, 0.00000E+00,-3.29933E+02, 0.00000E+00, + -1.10080E-03,-9.39310E-05, 1.39638E-04, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +} }; + +double pma[10][100] = { +/* TN2(2) */ { + 9.81637E-01,-1.41317E-03, 3.87323E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.58707E-02, + -8.63658E-03, 0.00000E+00, 0.00000E+00,-2.02226E+00, 0.00000E+00, + -8.69424E-03,-1.91397E-02, 8.76779E+01, 4.52188E-03, 0.00000E+00, + 2.23760E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-7.07572E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, + -4.11210E-03, 3.50060E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-8.36657E-03, 1.61347E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-1.45130E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.24152E-03, + 6.43365E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.33255E-03, + 2.42657E-03, 1.60666E-03,-1.85728E-03,-1.46874E-03,-4.79163E-06, + 1.22464E+00, 3.53510E+01, 0.00000E+00, 4.49223E-01,-4.77466E+01, + 4.70681E-01, 8.41861E+00,-2.88198E-01, 1.67854E+02, 0.00000E+00, + 7.11493E-04, 6.05601E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN2(3) */ { + 1.00422E+00,-7.11212E-03, 5.24480E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-5.28914E-02, + -2.41301E-02, 0.00000E+00, 0.00000E+00,-2.12219E+01,-1.03830E-02, + -3.28077E-03, 1.65727E-02, 1.68564E+00,-6.68154E-03, 0.00000E+00, + 1.45155E-02, 0.00000E+00, 8.42365E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-4.34645E-03, 0.00000E+00, 0.00000E+00, 2.16780E-02, + 0.00000E+00,-1.38459E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 7.04573E-03,-4.73204E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 1.08767E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-8.08279E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.21769E-04, + -2.27387E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.26769E-03, + 3.16901E-03, 4.60316E-04,-1.01431E-04, 1.02131E-03, 9.96601E-04, + 1.25707E+00, 2.50114E+01, 0.00000E+00, 4.24472E-01,-2.77655E+01, + 3.44625E-01, 2.75412E+01, 0.00000E+00, 7.94251E+02, 0.00000E+00, + 2.45835E-03, 1.38871E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN2(4) TN3(1) */ { + 1.01890E+00,-2.46603E-02, 1.00078E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-6.70977E-02, + -4.02286E-02, 0.00000E+00, 0.00000E+00,-2.29466E+01,-7.47019E-03, + 2.26580E-03, 2.63931E-02, 3.72625E+01,-6.39041E-03, 0.00000E+00, + 9.58383E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-1.85291E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 1.39717E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 9.19771E-03,-3.69121E+02, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-1.57067E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.07265E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.92953E-03, + -2.77739E-03,-4.40092E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.47280E-03, + 2.95035E-04,-1.81246E-03, 2.81945E-03, 4.27296E-03, 9.78863E-04, + 1.40545E+00,-6.19173E+00, 0.00000E+00, 0.00000E+00,-7.93632E+01, + 4.44643E-01,-4.03085E+02, 0.00000E+00, 1.15603E+01, 0.00000E+00, + 2.25068E-03, 8.48557E-04,-2.98493E-04, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN3(2) */ { + 9.75801E-01, 3.80680E-02,-3.05198E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.85575E-02, + 5.04057E-02, 0.00000E+00, 0.00000E+00,-1.76046E+02, 1.44594E-02, + -1.48297E-03,-3.68560E-03, 3.02185E+01,-3.23338E-03, 0.00000E+00, + 1.53569E-02, 0.00000E+00,-1.15558E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 4.89620E-03, 0.00000E+00, 0.00000E+00,-1.00616E-02, + -8.21324E-03,-1.57757E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 6.63564E-03, 4.58410E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-2.51280E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 9.91215E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-8.73148E-04, + -1.29648E-03,-7.32026E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-4.68110E-03, + -4.66003E-03,-1.31567E-03,-7.39390E-04, 6.32499E-04,-4.65588E-04, + -1.29785E+00,-1.57139E+02, 0.00000E+00, 2.58350E-01,-3.69453E+01, + 4.10672E-01, 9.78196E+00,-1.52064E-01,-3.85084E+03, 0.00000E+00, + -8.52706E-04,-1.40945E-03,-7.26786E-04, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN3(3) */ { + 9.60722E-01, 7.03757E-02,-3.00266E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.22671E-02, + 4.10423E-02, 0.00000E+00, 0.00000E+00,-1.63070E+02, 1.06073E-02, + 5.40747E-04, 7.79481E-03, 1.44908E+02, 1.51484E-04, 0.00000E+00, + 1.97547E-02, 0.00000E+00,-1.41844E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 5.77884E-03, 0.00000E+00, 0.00000E+00, 9.74319E-03, + 0.00000E+00,-2.88015E+03, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-4.44902E-03,-2.92760E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 2.34419E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.36685E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-4.65325E-04, + -5.50628E-04, 3.31465E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.06179E-03, + -3.08575E-03,-7.93589E-04,-1.08629E-04, 5.95511E-04,-9.05050E-04, + 1.18997E+00, 4.15924E+01, 0.00000E+00,-4.72064E-01,-9.47150E+02, + 3.98723E-01, 1.98304E+01, 0.00000E+00, 3.73219E+03, 0.00000E+00, + -1.50040E-03,-1.14933E-03,-1.56769E-04, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN3(4) */ { + 1.03123E+00,-7.05124E-02, 8.71615E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.82621E-02, + -9.80975E-03, 0.00000E+00, 0.00000E+00, 2.89286E+01, 9.57341E-03, + 0.00000E+00, 0.00000E+00, 8.66153E+01, 7.91938E-04, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 4.68917E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 7.86638E-03, 0.00000E+00, 0.00000E+00, 9.90827E-03, + 0.00000E+00, 6.55573E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00,-4.00200E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 7.07457E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.72268E-03, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.04970E-04, + 1.21560E-03,-8.05579E-06, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.49941E-03, + -4.57256E-04,-1.59311E-04, 2.96481E-04,-1.77318E-03,-6.37918E-04, + 1.02395E+00, 1.28172E+01, 0.00000E+00, 1.49903E-01,-2.63818E+01, + 0.00000E+00, 4.70628E+01,-2.22139E-01, 4.82292E-02, 0.00000E+00, + -8.67075E-04,-5.86479E-04, 5.32462E-04, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TN3(5) SURFACE TEMP TSL */ { + 1.00828E+00,-9.10404E-02,-2.26549E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.32420E-02, + -9.08925E-03, 0.00000E+00, 0.00000E+00, 3.36105E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-1.24957E+01,-5.87939E-03, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.79765E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.01237E+03, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-1.75553E-02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.29699E-03, + 1.26659E-03, 2.68402E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.17894E-03, + 1.48746E-03, 1.06478E-04, 1.34743E-04,-2.20939E-03,-6.23523E-04, + 6.36539E-01, 1.13621E+01, 0.00000E+00,-3.93777E-01, 2.38687E+03, + 0.00000E+00, 6.61865E+02,-1.21434E-01, 9.27608E+00, 0.00000E+00, + 1.68478E-04, 1.24892E-03, 1.71345E-03, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TGN3(2) SURFACE GRAD TSLG */ { + 1.57293E+00,-6.78400E-01, 6.47500E-01, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.62974E-02, + -3.60423E-01, 0.00000E+00, 0.00000E+00, 1.28358E+02, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 4.68038E+01, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-1.67898E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.90994E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.15706E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TGN2(1) TGN1(2) */ { + 8.60028E-01, 3.77052E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.17570E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 7.77757E-03, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 1.01024E+02, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 6.54251E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.56959E-02, + 1.91001E-02, 3.15971E-02, 1.00982E-02,-6.71565E-03, 2.57693E-03, + 1.38692E+00, 2.82132E-01, 0.00000E+00, 0.00000E+00, 3.81511E+02, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +}, /* TGN3(1) TGN2(2) */ { + 1.06029E+00,-5.25231E-02, 3.73034E-01, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.31072E-02, + -3.88409E-01, 0.00000E+00, 0.00000E+00,-1.65295E+02,-2.13801E-01, + -4.38916E-02,-3.22716E-01,-8.82393E+01, 1.18458E-01, 0.00000E+00, + -4.35863E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00,-1.19782E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 2.62229E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00,-5.37443E+01, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00,-4.55788E-01, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.84009E-02, + 3.96733E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.05494E-02, + 7.39617E-02, 1.92200E-02,-8.46151E-03,-1.34244E-02, 1.96338E-02, + 1.50421E+00, 1.88368E+01, 0.00000E+00, 0.00000E+00,-5.13114E+01, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 5.11923E-02, 3.61225E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00 +} }; + +/* SEMIANNUAL MULT SAM */ +double sam[100] = { + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, + 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 +}; + + +/* MIDDLE ATMOSPHERE AVERAGES */ +double pavgm[10] = { + 2.61000E+02, 2.64000E+02, 2.29000E+02, 2.17000E+02, 2.17000E+02, + 2.23000E+02, 2.86760E+02,-2.93940E+00, 2.50000E+00, 0.00000E+00 }; + diff --git a/nrlmsise-00_test.c b/nrlmsise-00_test.c new file mode 100644 index 0000000..28f547f --- /dev/null +++ b/nrlmsise-00_test.c @@ -0,0 +1,168 @@ +/* -------------------------------------------------------------------- */ +/* --------- N R L M S I S E - 0 0 M O D E L 2 0 0 1 ---------- */ +/* -------------------------------------------------------------------- */ + +/* This file is part of the NRLMSISE-00 C source code package - release + * 20020305 + * + * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and + * Doug Drob. They also wrote a NRLMSISE-00 distribution package in + * FORTRAN which is available at + * http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm + * + * Dominik Brodowski implemented and maintains this C version. You can + * reach him at devel@brodo.de. See the file "DOCUMENTATION" for details, + * and check http://www.brodo.de/english/pub/nrlmsise/index.html for + * updated releases of this package. + */ + + + +/* ------------------------------------------------------------------- */ +/* ------------------------------ INCLUDES --------------------------- */ +/* ------------------------------------------------------------------- */ + +#include <stdio.h> +#include "nrlmsise-00.h" + + + +/* ------------------------------------------------------------------- */ +/* ----------------------------- TEST_GTD7 --------------------------- */ +/* ------------------------------------------------------------------- */ + +void test_gtd7(void) { + struct nrlmsise_output output[17]; + struct nrlmsise_input input[17]; + struct nrlmsise_flags flags; + struct ap_array aph; + int i; + int j; + /* input values */ + for (i=0;i<7;i++) + aph.a[i]=100; + flags.switches[0]=0; + for (i=1;i<24;i++) + flags.switches[i]=1; + for (i=0;i<17;i++) { + input[i].doy=172; + input[i].year=0; /* without effect */ + input[i].sec=29000; + input[i].alt=400; + input[i].g_lat=60; + input[i].g_long=-70; + input[i].lst=16; + input[i].f107A=150; + input[i].f107=150; + input[i].ap=4; + } + input[1].doy=81; + input[2].sec=75000; + input[2].alt=1000; + input[3].alt=100; + input[10].alt=0; + input[11].alt=10; + input[12].alt=30; + input[13].alt=50; + input[14].alt=70; + input[16].alt=100; + input[4].g_lat=0; + input[5].g_long=0; + input[6].lst=4; + input[7].f107A=70; + input[8].f107=180; + input[9].ap=40; + input[15].ap_a=&aph; + input[16].ap_a=&aph; + /* evaluate 0 to 14 */ + for (i=0;i<15;i++) + gtd7(&input[i], &flags, &output[i]); + /* evaluate 15 and 16 */ + flags.switches[9]=-1; + for (i=15;i<17;i++) + gtd7(&input[i], &flags, &output[i]); + /* output type 1 */ + for (i=0;i<17;i++) { + printf("\n"); + for (j=0;j<9;j++) + printf("%E ",output[i].d[j]); + printf("%E ",output[i].t[0]); + printf("%E \n",output[i].t[1]); + /* DL omitted */ + } + + /* output type 2 */ + for (i=0;i<3;i++) { + printf("\n"); + printf("\nDAY "); + for (j=0;j<5;j++) + printf(" %3i",input[i*5+j].doy); + printf("\nUT "); + for (j=0;j<5;j++) + printf(" %5.0f",input[i*5+j].sec); + printf("\nALT "); + for (j=0;j<5;j++) + printf(" %4.0f",input[i*5+j].alt); + printf("\nLAT "); + for (j=0;j<5;j++) + printf(" %3.0f",input[i*5+j].g_lat); + printf("\nLONG "); + for (j=0;j<5;j++) + printf(" %3.0f",input[i*5+j].g_long); + printf("\nLST "); + for (j=0;j<5;j++) + printf(" %5.0f",input[i*5+j].lst); + printf("\nF107A "); + for (j=0;j<5;j++) + printf(" %3.0f",input[i*5+j].f107A); + printf("\nF107 "); + for (j=0;j<5;j++) + printf(" %3.0f",input[i*5+j].f107); + printf("\n\n"); + printf("\nTINF "); + for (j=0;j<5;j++) + printf(" %7.2f",output[i*5+j].t[0]); + printf("\nTG "); + for (j=0;j<5;j++) + printf(" %7.2f",output[i*5+j].t[1]); + printf("\nHE "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[0]); + printf("\nO "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[1]); + printf("\nN2 "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[2]); + printf("\nO2 "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[3]); + printf("\nAR "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[4]); + printf("\nH "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[6]); + printf("\nN "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[7]); + printf("\nANM 0 "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[8]); + printf("\nRHO "); + for (j=0;j<5;j++) + printf(" %1.3e",output[i*5+j].d[5]); + printf("\n"); + } + printf("\n"); +} + + +/* ------------------------------------------------------------------- */ +/* -------------------------------- MAIN ----------------------------- */ +/* ------------------------------------------------------------------- */ + +int main(void) { + test_gtd7(); + return 0; +} |