Update C code to include changes up to Official Beta Release 2.0

Update C code to include changes up to Official Beta Release 2.0
(NRLMSISE-00.DIST17.TXT).

Signed-off-by: Dominik Brodowski <linux@dominikbrodowski.net>

1 files changed, 37 insertions, 6 deletions

diff --git a/nrlmsise-00.c b/nrlmsise-00.c
index 73a8062..17413d0 100644
--- a/nrlmsise-00.c
+++ b/nrlmsise-00.c
@@ -3,7 +3,7 @@
 /* -------------------------------------------------------------------- */
 
 /* This file is part of the NRLMSISE-00  C source code package - release
- * 20020305
+ * 20020503
  *
  * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and
  * Doug Drob. They also wrote a NRLMSISE-00 distribution package in 
@@ -140,6 +140,35 @@ double ccor(double alt, double r, double h1, double zh) {
 
 
 /* ------------------------------------------------------------------- */
+/* ------------------------------ CCOR ------------------------------- */
+/* ------------------------------------------------------------------- */
+
+double ccor2(double alt, double r, double h1, double zh, double h2) {
+/*        CHEMISTRY/DISSOCIATION CORRECTION FOR MSIS MODELS
+ *         ALT - altitude
+ *         R - target ratio
+ *         H1 - transition scale length
+ *         ZH - altitude of 1/2 R
+ *         H2 - transition scale length #2 ?
+ */
+	double e1, e2;
+	double ex1, ex2;
+	double ccor2v;
+	e1 = (alt - zh) / h1;
+	e2 = (alt - zh) / h2;
+	if ((e1 > 70) || (e2 > 70))
+		return exp(0);
+	if ((e1 < -70) && (e2 < -70))
+		return exp(r);
+	ex1 = exp(e1);
+	ex2 = exp(e2);
+	ccor2v = r / (1.0 + 0.5 * (ex1 + ex2));
+	return exp(ccor2v);
+}
+
+
+
+/* ------------------------------------------------------------------- */
 /* ------------------------------- SCALH ----------------------------- */
 /* ------------------------------------------------------------------- */
 
@@ -1145,7 +1174,6 @@ void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrl
 	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;
@@ -1157,6 +1185,7 @@ void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrl
 	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;
+	double hc216, hcc232;
 	za = pdl[1][15];
 	zn1[0] = za;
 	for (j=0;j<9;j++) 
@@ -1239,7 +1268,7 @@ void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrl
         /*  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])) {
+	if ((flags->sw[15]) && (z<altl[0])) {
 		/*  Turbopause */
 		zh04=pdm[0][2];
 		/*  Mixed density at Zlb */
@@ -1252,7 +1281,7 @@ void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrl
 		/*  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];
+		hc04=pdm[0][5]*pdl[1][1];
 		/*  Net density corrected at Alt */
 		output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
 	}
@@ -1280,7 +1309,8 @@ void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrl
 		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);
+		hc216=pdm[1][5]*pdl[1][4];
+		output->d[1]=output->d[1]*ccor2(z,rl,hc16,zc16,hc216);
 		/*   Chemistry correction */
 		hcc16=pdm[1][7]*pdl[1][13];
 		zcc16=pdm[1][6]*pdl[1][12];
@@ -1318,10 +1348,11 @@ void gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrl
 		}
 		/*  Correction for general departure from diffusive equilibrium above Zlb */
 		hcc32=pdm[3][7]*pdl[1][22];
+		hcc232=pdm[3][7]*pdl[0][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);
+		output->d[3]=output->d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
 	}