src/FDM/JSBSim/models/atmosphere/FGMSIS.cpp

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