Msise00Drag.cpp

Go to the documentation of this file.
00001 #pragma ident "$Id: Msise00Drag.cpp 2457 2010-08-18 14:20:12Z coandrei $"
00002 
00009 //============================================================================
00010 //
00011 //  This file is part of GPSTk, the GPS Toolkit.
00012 //
00013 //  The GPSTk is free software; you can redistribute it and/or modify
00014 //  it under the terms of the GNU Lesser General Public License as published
00015 //  by the Free Software Foundation; either version 2.1 of the License, or
00016 //  any later version.
00017 //
00018 //  The GPSTk is distributed in the hope that it will be useful,
00019 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00020 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00021 //  GNU Lesser General Public License for more details.
00022 //
00023 //  You should have received a copy of the GNU Lesser General Public
00024 //  License along with GPSTk; if not, write to the Free Software Foundation,
00025 //  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00026 //
00027 //  Wei Yan - Chinese Academy of Sciences . 2009, 2010
00028 //
00029 //============================================================================
00030 
00031 
00032 #include "Msise00Drag.hpp"
00033 #include <string>
00034 #include <cmath>
00035 #include <stdio.h>
00036 #include <stdlib.h>
00037 #include "ReferenceFrames.hpp"
00038 #include "IERS.hpp"
00039 #include "ASConstant.hpp"
00040 
00041 namespace gpstk
00042 {
00043 
00044    // test the model
00045    void Msise00Drag::test()
00046    {
00047       struct nrlmsise_output output[17];
00048       struct nrlmsise_input input[17];
00049       struct nrlmsise_flags flags;
00050       struct ap_array aph;
00051 
00052       int i;
00053       int j;
00054       
00055       /* input values */
00056       for (i=0;i<7;i++)
00057          aph.a[i]=100;
00058 
00059       flags.switches[0]=0;
00060       for (i=1;i<24;i++)
00061          flags.switches[i]=1;
00062       for (i=0;i<17;i++) {
00063          input[i].doy=172;
00064          input[i].year=0; /* without effect */
00065          input[i].sec=29000;
00066          input[i].alt=400;
00067          input[i].g_lat=60;
00068          input[i].g_long=-70;
00069          input[i].lst=16;
00070          input[i].f107A=150;
00071          input[i].f107=150;
00072          input[i].ap=4;
00073       }
00074       input[1].doy=81;
00075       input[2].sec=75000;
00076       input[2].alt=1000;
00077       input[3].alt=100;
00078       input[10].alt=0;
00079       input[11].alt=10;
00080       input[12].alt=30;
00081       input[13].alt=50;
00082       input[14].alt=70;
00083       input[16].alt=100;
00084       input[4].g_lat=0;
00085       input[5].g_long=0;
00086       input[6].lst=4;
00087       input[7].f107A=70;
00088       input[8].f107=180;
00089       input[9].ap=40;
00090       input[15].ap_a=&aph;
00091       input[16].ap_a=&aph;
00092       /* evaluate 0 to 14 */
00093       for (i=0;i<15;i++)
00094          gtd7(&input[i], &flags, &output[i]);
00095       /* evaluate 15 and 16 */
00096       flags.switches[9]=-1;
00097       for (i=15;i<17;i++)
00098          gtd7(&input[i], &flags, &output[i]);
00099       /* output type 1 */
00100       for (i=0;i<17;i++) {
00101          printf("\n");
00102          for (j=0;j<9;j++)
00103             printf("%E ",output[i].d[j]);
00104          printf("%E ",output[i].t[0]);
00105          printf("%E \n",output[i].t[1]);
00106          /* DL omitted */
00107       }
00108 
00109       /* output type 2 */
00110       for (i=0;i<3;i++) 
00111       {
00112          printf("\n");
00113          printf("\nDAY   ");
00114          for (j=0;j<5;j++)
00115             printf("         %3i",input[i*5+j].doy);
00116          printf("\nUT    ");
00117          for (j=0;j<5;j++)
00118             printf("       %5.0f",input[i*5+j].sec);
00119          printf("\nALT   ");
00120          for (j=0;j<5;j++)
00121             printf("        %4.0f",input[i*5+j].alt);
00122          printf("\nLAT   ");
00123          for (j=0;j<5;j++)
00124             printf("         %3.0f",input[i*5+j].g_lat);
00125          printf("\nLONG  ");
00126          for (j=0;j<5;j++)
00127             printf("         %3.0f",input[i*5+j].g_long);
00128          printf("\nLST   ");
00129          for (j=0;j<5;j++)
00130             printf("       %5.0f",input[i*5+j].lst);
00131          printf("\nF107A ");
00132          for (j=0;j<5;j++)
00133             printf("         %3.0f",input[i*5+j].f107A);
00134          printf("\nF107  ");
00135          for (j=0;j<5;j++)
00136             printf("         %3.0f",input[i*5+j].f107);
00137          printf("\n\n");
00138          printf("\nTINF  ");
00139          for (j=0;j<5;j++)
00140             printf("     %7.2f",output[i*5+j].t[0]);
00141          printf("\nTG    ");
00142          for (j=0;j<5;j++)
00143             printf("     %7.2f",output[i*5+j].t[1]);
00144          printf("\nHE    ");
00145          for (j=0;j<5;j++)
00146             printf("   %1.3e",output[i*5+j].d[0]);
00147          printf("\nO     ");
00148          for (j=0;j<5;j++)
00149             printf("   %1.3e",output[i*5+j].d[1]);
00150          printf("\nN2    ");
00151          for (j=0;j<5;j++)
00152             printf("   %1.3e",output[i*5+j].d[2]);
00153          printf("\nO2    ");
00154          for (j=0;j<5;j++)
00155             printf("   %1.3e",output[i*5+j].d[3]);
00156          printf("\nAR    ");
00157          for (j=0;j<5;j++)
00158             printf("   %1.3e",output[i*5+j].d[4]);
00159          printf("\nH     ");
00160          for (j=0;j<5;j++)
00161             printf("   %1.3e",output[i*5+j].d[6]);
00162          printf("\nN     ");
00163          for (j=0;j<5;j++)
00164             printf("   %1.3e",output[i*5+j].d[7]);
00165          printf("\nANM 0 ");
00166          for (j=0;j<5;j++)
00167             printf("   %1.3e",output[i*5+j].d[8]);
00168          printf("\nRHO   ");
00169          for (j=0;j<5;j++)
00170             printf("   %1.3e",output[i*5+j].d[5]);
00171          printf("\n");
00172       }
00173       printf("\n");
00174       
00175       /*
00176       nrlmsise-test should generate the following output:
00177 
00178       6.665177E+05 1.138806E+08 1.998211E+07 4.022764E+05 3.557465E+03 4.074714E-15 3.475312E+04 4.095913E+06 2.667273E+04 1.250540E+03 1.241416E+03 
00179 
00180       3.407293E+06 1.586333E+08 1.391117E+07 3.262560E+05 1.559618E+03 5.001846E-15 4.854208E+04 4.380967E+06 6.956682E+03 1.166754E+03 1.161710E+03 
00181 
00182       1.123767E+05 6.934130E+04 4.247105E+01 1.322750E-01 2.618848E-05 2.756772E-18 2.016750E+04 5.741256E+03 2.374394E+04 1.239892E+03 1.239891E+03 
00183 
00184       5.411554E+07 1.918893E+11 6.115826E+12 1.225201E+12 6.023212E+10 3.584426E-10 1.059880E+07 2.615737E+05 2.819879E-42 1.027318E+03 2.068878E+02 
00185 
00186       1.851122E+06 1.476555E+08 1.579356E+07 2.633795E+05 1.588781E+03 4.809630E-15 5.816167E+04 5.478984E+06 1.264446E+03 1.212396E+03 1.208135E+03 
00187 
00188       8.673095E+05 1.278862E+08 1.822577E+07 2.922214E+05 2.402962E+03 4.355866E-15 3.686389E+04 3.897276E+06 2.667273E+04 1.220146E+03 1.212712E+03 
00189 
00190       5.776251E+05 6.979139E+07 1.236814E+07 2.492868E+05 1.405739E+03 2.470651E-15 5.291986E+04 1.069814E+06 2.667273E+04 1.116385E+03 1.112999E+03 
00191 
00192       3.740304E+05 4.782720E+07 5.240380E+06 1.759875E+05 5.501649E+02 1.571889E-15 8.896776E+04 1.979741E+06 9.121815E+03 1.031247E+03 1.024848E+03 
00193 
00194       6.748339E+05 1.245315E+08 2.369010E+07 4.911583E+05 4.578781E+03 4.564420E-15 3.244595E+04 5.370833E+06 2.667273E+04 1.306052E+03 1.293374E+03 
00195 
00196       5.528601E+05 1.198041E+08 3.495798E+07 9.339618E+05 1.096255E+04 4.974543E-15 2.686428E+04 4.889974E+06 2.805445E+04 1.361868E+03 1.347389E+03 
00197 
00198       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 
00199 
00200       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 
00201 
00202       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 
00203 
00204       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 
00205 
00206       1.254884E+10 0.000000E+00 1.874533E+15 4.923051E+14 2.239685E+13 1.147668E-07 0.000000E+00 0.000000E+00 0.000000E+00 1.027318E+03 2.190732E+02 
00207 
00208       5.196477E+05 1.274494E+08 4.850450E+07 1.720838E+06 2.354487E+04 5.881940E-15 2.500078E+04 6.279210E+06 2.667273E+04 1.426412E+03 1.408608E+03 
00209 
00210       4.260860E+07 1.241342E+11 4.929562E+12 1.048407E+12 4.993465E+10 2.914304E-10 8.831229E+06 2.252516E+05 2.415246E-42 1.027318E+03 1.934071E+02 
00211 
00212 
00213       DAY            172          81         172         172         172
00214       UT           29000       29000       75000       29000       29000
00215       ALT            400         400        1000         100         400
00216       LAT             60          60          60          60           0
00217       LONG           -70         -70         -70         -70         -70
00218       LST             16          16          16          16          16
00219       F107A          150         150         150         150         150
00220       F107           150         150         150         150         150
00221 
00222 
00223       TINF       1250.54     1166.75     1239.89     1027.32     1212.40
00224       TG         1241.42     1161.71     1239.89      206.89     1208.14
00225       HE       6.665e+05   3.407e+06   1.124e+05   5.412e+07   1.851e+06
00226       O        1.139e+08   1.586e+08   6.934e+04   1.919e+11   1.477e+08
00227       N2       1.998e+07   1.391e+07   4.247e+01   6.116e+12   1.579e+07
00228       O2       4.023e+05   3.263e+05   1.323e-01   1.225e+12   2.634e+05
00229       AR       3.557e+03   1.560e+03   2.619e-05   6.023e+10   1.589e+03
00230       H        3.475e+04   4.854e+04   2.017e+04   1.060e+07   5.816e+04
00231       N        4.096e+06   4.381e+06   5.741e+03   2.616e+05   5.479e+06
00232       ANM 0    2.667e+04   6.957e+03   2.374e+04   2.820e-42   1.264e+03
00233       RHO      4.075e-15   5.002e-15   2.757e-18   3.584e-10   4.810e-15
00234 
00235 
00236       DAY            172         172         172         172         172
00237       UT           29000       29000       29000       29000       29000
00238       ALT            400         400         400         400         400
00239       LAT             60          60          60          60          60
00240       LONG             0         -70         -70         -70         -70
00241       LST             16           4          16          16          16
00242       F107A          150         150          70         150         150
00243       F107           150         150         150         180         150
00244 
00245 
00246       TINF       1220.15     1116.39     1031.25     1306.05     1361.87
00247       TG         1212.71     1113.00     1024.85     1293.37     1347.39
00248       HE       8.673e+05   5.776e+05   3.740e+05   6.748e+05   5.529e+05
00249       O        1.279e+08   6.979e+07   4.783e+07   1.245e+08   1.198e+08
00250       N2       1.823e+07   1.237e+07   5.240e+06   2.369e+07   3.496e+07
00251       O2       2.922e+05   2.493e+05   1.760e+05   4.912e+05   9.340e+05
00252       AR       2.403e+03   1.406e+03   5.502e+02   4.579e+03   1.096e+04
00253       H        3.686e+04   5.292e+04   8.897e+04   3.245e+04   2.686e+04
00254       N        3.897e+06   1.070e+06   1.980e+06   5.371e+06   4.890e+06
00255       ANM 0    2.667e+04   2.667e+04   9.122e+03   2.667e+04   2.805e+04
00256       RHO      4.356e-15   2.471e-15   1.572e-15   4.564e-15   4.975e-15
00257 
00258 
00259       DAY            172         172         172         172         172
00260       UT           29000       29000       29000       29000       29000
00261       ALT              0          10          30          50          70
00262       LAT             60          60          60          60          60
00263       LONG           -70         -70         -70         -70         -70
00264       LST             16          16          16          16          16
00265       F107A          150         150         150         150         150
00266       F107           150         150         150         150         150
00267 
00268 
00269       TINF       1027.32     1027.32     1027.32     1027.32     1027.32
00270       TG          281.46      227.42      237.44      279.56      219.07
00271       HE       1.375e+14   4.427e+13   2.128e+12   1.412e+11   1.255e+10
00272       O        0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
00273       N2       2.050e+19   6.598e+18   3.171e+17   2.104e+16   1.875e+15
00274       O2       5.499e+18   1.770e+18   8.506e+16   5.645e+15   4.923e+14
00275       AR       2.452e+17   7.892e+16   3.793e+15   2.517e+14   2.240e+13
00276       H        0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
00277       N        0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
00278       ANM 0    0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
00279       RHO      1.261e-03   4.059e-04   1.951e-05   1.295e-06   1.148e-07
00280 
00281 
00282       Note: These values equal those of the official FORTRAN package with
00283       one notable exception: the FORTRAN version reports for "anomalous
00284       oxygen" in test-run 4 exactly 0.000E-00, while my C compiler
00285       generates code which calculates 2.820e-42. When only 16-bit wide
00286       double variables are used, this value reduces to 0.000E-00 as well.
00287 
00288       */
00289    }
00290 
00291    
00292       /* Abstract class requires the subclass to compute the atmospheric density.
00293        * @param utc epoch in UTC
00294        * @param rb  EarthRef object.
00295        * @param r   Position vector.
00296        * @param v   Velocity vector
00297        * @return Atmospheric density in kg/m^3
00298        */
00299    double Msise00Drag::computeDensity(UTCTime utc, 
00300                                       EarthBody& rb, 
00301                                       Vector<double> r, 
00302                                       Vector<double> v)
00303    {
00304       struct nrlmsise_output output;
00305       struct nrlmsise_input input; 
00306       struct nrlmsise_flags flags;
00307       struct ap_array aph;
00308 
00309       //* Get the J2000 to TOD transformation
00310       Matrix<double> N = ReferenceFrames::J2kToTODMatrix(utc);
00311 
00312       //* Transform r from J2000 to TOD
00313       Vector<double> r_tod = N * r;
00314 
00315 
00316       Matrix<double> eci2ecef = ReferenceFrames::J2kToECEFMatrix(utc);
00317 
00318       Vector<double> r_ecef = eci2ecef * r;
00319       
00320       Position geoidPos(r_ecef(0),r_ecef(1),r_ecef(3),Position::Cartesian);
00321       double alt = geoidPos.getAltitude() / 1000.0;    //* [km]
00322 
00323       if (alt > 1000.0) 
00324       {
00325          string msg("Msise00Drag only valid from 0 to 1000 km");
00326          Exception e(msg);
00327          //GPSTK_THROW(e);
00328       }
00329       
00330       //double dist2sun = norm(ReferenceFrames::getJ2kPosition(utc.asTDB(), SolarSystem::Sun))*1000.0;
00331       //double f107_in = this.f107_opt*Math.pow(dist2sun/Constants.AU,2);
00332       //double f107_in = this.f107_opt*Math.pow(Constants.AU/dist2sun,2);
00333       double f107_in = this->f107_opt;
00334             
00335       /* input values */
00336       //for (i=0;i<7;i++)
00337       //   aph.a[i]=13.853964381;//100;
00338 
00339       flags.switches[0]=0;
00340       for (int i=1; i<24; i++)
00341       {
00342          flags.switches[i] = 1;
00343       }
00344 
00345       input.doy = utc.DOYday();
00346       input.year = 2004;         // without effect 
00347       input.sec = utc.DOYsecond();
00348       input.alt= alt;
00349       input.g_lat = geoidPos.getGeodeticLatitude();
00350       input.g_long = geoidPos.getLongitude();
00351       input.lst = input.sec/3600.0 + input.g_long/15.0;
00352       input.f107A =f107_in;
00353       input.f107 = f107_in;
00354       input.ap = this->ap_opt;    //14.924291;//13.853964381; //???
00355 
00356       if(alt > 500.0)
00357       {
00358          gtd7d(&input, &flags, &output);
00359       }else
00360       {
00361          gtd7(&input, &flags, &output);
00362       }
00363 
00364       return output.d[5]*1000.0; //[kg/m^3]  
00365 
00366    }  // End of method 'Msise00Drag::computeDensity()'
00367 
00368 
00369    // LOCAL FUNCTIONS
00370    //-----------------------------------------------------------------------
00371 
00372    void Msise00Drag::tselec(struct nrlmsise_flags *flags) 
00373    {
00374       for (int i=0; i<24; i++) 
00375       {
00376          if (i!=9) 
00377          {
00378             if (flags->switches[i]==1)
00379                flags->sw[i]=1;
00380             else
00381                flags->sw[i]=0;
00382             if (flags->switches[i]>0)
00383                flags->swc[i]=1;
00384             else
00385                flags->swc[i]=0;
00386          } 
00387          else
00388          {
00389             flags->sw[i]=flags->switches[i];
00390             flags->swc[i]=flags->switches[i];
00391          }
00392       }
00393 
00394    }  // 'tselec()'
00395 
00396 
00397    void Msise00Drag::glatf(double lat, double *gv, double *reff) 
00398    {
00399       double dgtr = 1.74533E-2;
00400       double c2;
00401       c2 = std::cos(2.0*dgtr*lat);
00402       *gv = 980.616 * (1.0 - 0.0026373 * c2);
00403       *reff = 2.0 * (*gv) / (3.085462E-6 + 2.27E-9 * c2) * 1.0E-5;
00404 
00405    }  // 'glatf()'
00406 
00407 
00408    double Msise00Drag::ccor(double alt, double r, double h1, double zh) 
00409    {
00410       double e;
00411       double ex;
00412       e = (alt - zh) / h1;
00413       if (e>70)
00414          return std::exp(0.0);
00415       if (e<-70)
00416          return std::exp(r);
00417       ex = std::exp(e);
00418       e = r / (1.0 + ex);
00419       return std::exp(e);
00420 
00421    }  // 'ccor()'
00422 
00423 
00424    double Msise00Drag::ccor2(double alt, double r, double h1, double zh, double h2)
00425    {
00426       double e1, e2;
00427       double ex1, ex2;
00428       double ccor2v;
00429       e1 = (alt - zh) / h1;
00430       e2 = (alt - zh) / h2;
00431       if ((e1 > 70) || (e2 > 70))
00432          return std::exp(0.0);
00433       if ((e1 < -70) && (e2 < -70))
00434          return std::exp(r);
00435       ex1 = std::exp(e1);
00436       ex2 = std::exp(e2);
00437       ccor2v = r / (1.0 + 0.5 * (ex1 + ex2));
00438       return std::exp(ccor2v);
00439 
00440    }  // 'ccor2()'
00441 
00442 
00443    double Msise00Drag::scalh(double alt, double xm, double temp)
00444    {
00445       double g;
00446       double rgas=831.4;
00447       g = gsurf / (std::pow((1.0 + alt/re),2.0));
00448       g = rgas * temp / (g * xm);
00449       return g;
00450 
00451    }  // 'scalh()'
00452 
00453 
00454    double Msise00Drag::dnet (double dd, double dm, double zhm, double xmm, double xm)
00455    {
00456       double a;
00457       double ylog;
00458       a  = zhm / (xmm-xm);
00459       if (!((dm>0) && (dd>0))) 
00460       {
00461          printf("dnet log error %e %e %e\n",dm,dd,xm);
00462          if ((dd==0) && (dm==0))
00463             dd=1;
00464          if (dm==0)
00465             return dd;
00466          if (dd==0)
00467             return dm;
00468       } 
00469       ylog = a * std::log(dm/dd);
00470       if (ylog<-10)
00471          return dd;
00472       if (ylog>10)
00473          return dm;
00474       a = dd*std::pow((1.0 + std::exp(ylog)),(1.0/a));
00475       return a;
00476 
00477    }  // 'dnet()'
00478 
00479 
00480    void Msise00Drag::splini (double *xa, double *ya, double *y2a, int n, double x, double *y)
00481    {
00482       double yi=0;
00483       int klo=0;
00484       int khi=1;
00485       double xx, h, a, b, a2, b2;
00486       while ((x>xa[klo]) && (khi<n))
00487       {
00488          xx=x;
00489          if (khi<(n-1)) 
00490          {
00491             if (x<xa[khi])
00492                xx=x;
00493             else 
00494                xx=xa[khi];
00495          }
00496          h = xa[khi] - xa[klo];
00497          a = (xa[khi] - xx)/h;
00498          b = (xx - xa[klo])/h;
00499          a2 = a*a;
00500          b2 = b*b;
00501          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;
00502          klo++;
00503          khi++;
00504       }
00505       *y = yi;
00506 
00507    }  // 'splini()'
00508 
00509 
00510    void Msise00Drag::splint (double *xa, double *ya, double *y2a, int n, double x, double *y)
00511    {
00512       /*      CALCULATE CUBIC SPLINE INTERP VALUE
00513       *       ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL.
00514       *       XA,YA: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
00515       *       Y2A: ARRAY OF SECOND DERIVATIVES
00516       *       N: SIZE OF ARRAYS XA,YA,Y2A
00517       *       X: ABSCISSA FOR INTERPOLATION
00518       *       Y: OUTPUT VALUE
00519       */
00520       int klo=0;
00521       int khi=n-1;
00522       int k;
00523       double h;
00524       double a, b, yi;
00525       while ((khi-klo)>1) 
00526       {
00527          k=(khi+klo)/2;
00528          if (xa[k]>x)
00529             khi=k;
00530          else
00531             klo=k;
00532       }
00533       h = xa[khi] - xa[klo];
00534       if (h==0.0)
00535          printf("bad XA input to splint");
00536       a = (xa[khi] - x)/h;
00537       b = (x - xa[klo])/h;
00538       yi = a * ya[klo] + b * ya[khi] + ((a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi]) * h * h/6.0;
00539       *y = yi;
00540 
00541    }  // 'splint()'
00542 
00543 
00544    void Msise00Drag::spline (double *x, double *y, int n, double yp1, double ypn, double *y2) 
00545    {
00546       /*       CALCULATE 2ND DERIVATIVES OF CUBIC SPLINE INTERP FUNCTION
00547       *       ADAPTED FROM NUMERICAL RECIPES BY PRESS ET AL
00548       *       X,Y: ARRAYS OF TABULATED FUNCTION IN ASCENDING ORDER BY X
00549       *       N: SIZE OF ARRAYS X,Y
00550       *       YP1,YPN: SPECIFIED DERIVATIVES AT X[0] AND X[N-1]; VALUES
00551       *                >= 1E30 SIGNAL SIGNAL SECOND DERIVATIVE ZERO
00552       *       Y2: OUTPUT ARRAY OF SECOND DERIVATIVES
00553       */
00554       double *u;
00555       double sig, p, qn, un;
00556       int i, k;
00557       u=(double*)malloc(sizeof(double)*n);
00558       if (u==NULL) 
00559       {
00560          printf("Out Of Memory in spline - ERROR");
00561          return;
00562       }
00563       if (yp1>0.99E30) 
00564       {
00565          y2[0]=0;
00566          u[0]=0;
00567       } else 
00568       {
00569          y2[0]=-0.5;
00570          u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
00571       }
00572       for (i=1;i<(n-1);i++) 
00573       {
00574          sig = (x[i]-x[i-1])/(x[i+1] - x[i-1]);
00575          p = sig * y2[i-1] + 2.0;
00576          y2[i] = (sig - 1.0) / p;
00577          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;
00578       }
00579       if (ypn>0.99E30) 
00580       {
00581          qn = 0;
00582          un = 0;
00583       } 
00584       else 
00585       {
00586          qn = 0.5;
00587          un = (3.0 / (x[n-1] - x[n-2])) * (ypn - (y[n-1] - y[n-2])/(x[n-1] - x[n-2]));
00588       }
00589       y2[n-1] = (un - qn * u[n-2]) / (qn * y2[n-2] + 1.0);
00590       for (k=n-2;k>=0;k--)
00591          y2[k] = y2[k] * y2[k+1] + u[k];
00592 
00593       free(u);
00594 
00595    }  // 'spline()'
00596 
00597 
00598    double Msise00Drag::zeta(double zz, double zl) 
00599    {
00600       return ((zz-zl)*(re+zl)/(re+zz));
00601    }
00602 
00603    double Msise00Drag::densm (double alt, double d0, double xm, double *tz, int mn3, 
00604       double *zn3, double *tn3, double *tgn3, int mn2, double *zn2, double *tn2, double *tgn2) 
00605    {
00606       /*      Calculate Temperature and Density Profiles for lower atmos.  */
00607       double xs[10], ys[10], y2out[10];
00608       double rgas = 831.4;
00609       double z, z1, z2, t1, t2, zg, zgdif;
00610       double yd1, yd2;
00611       double x, y, yi;
00612       double expl, gamm, glb;
00613       double densm_tmp;
00614       int mn;
00615       int k;
00616       densm_tmp=d0;
00617       if (alt>zn2[0]) 
00618       {
00619          if (xm==0.0)
00620             return *tz;
00621          else
00622             return d0;
00623       }
00624 
00625       /* STRATOSPHERE/MESOSPHERE TEMPERATURE */
00626       if (alt>zn2[mn2-1])
00627          z=alt;
00628       else
00629          z=zn2[mn2-1];
00630       mn=mn2;
00631       z1=zn2[0];
00632       z2=zn2[mn-1];
00633       t1=tn2[0];
00634       t2=tn2[mn-1];
00635       zg = zeta(z, z1);
00636       zgdif = zeta(z2, z1);
00637 
00638       /* set up spline nodes */
00639       for (k=0;k<mn;k++) 
00640       {
00641          xs[k]=zeta(zn2[k],z1)/zgdif;
00642          ys[k]=1.0 / tn2[k];
00643       }
00644       yd1=-tgn2[0] / (t1*t1) * zgdif;
00645       yd2=-tgn2[1] / (t2*t2) * zgdif * (std::pow(((re+z2)/(re+z1)),2.0));
00646 
00647       /* calculate spline coefficients */
00648       spline (xs, ys, mn, yd1, yd2, y2out);
00649       x = zg/zgdif;
00650       splint (xs, ys, y2out, mn, x, &y);
00651 
00652       /* temperature at altitude */
00653       *tz = 1.0 / y;
00654       if (xm!=0.0) 
00655       {
00656          /* calaculate stratosphere / mesospehere density */
00657          glb = gsurf / (std::pow((1.0 + z1/re),2.0));
00658          gamm = xm * glb * zgdif / rgas;
00659 
00660          /* Integrate temperature profile */
00661          splini(xs, ys, y2out, mn, x, &yi);
00662          expl=gamm*yi;
00663          if (expl>50.0)
00664             expl=50.0;
00665 
00666          /* Density at altitude */
00667          densm_tmp = densm_tmp * (t1 / *tz) * std::exp(-expl);
00668       }
00669 
00670       if (alt>zn3[0]) 
00671       {
00672          if (xm==0.0)
00673             return *tz;
00674          else
00675             return densm_tmp;
00676       }
00677 
00678       /* troposhere / stratosphere temperature */
00679       z = alt;
00680       mn = mn3;
00681       z1=zn3[0];
00682       z2=zn3[mn-1];
00683       t1=tn3[0];
00684       t2=tn3[mn-1];
00685       zg=zeta(z,z1);
00686       zgdif=zeta(z2,z1);
00687 
00688       /* set up spline nodes */
00689       for (k=0;k<mn;k++) 
00690       {
00691          xs[k] = zeta(zn3[k],z1) / zgdif;
00692          ys[k] = 1.0 / tn3[k];
00693       }
00694       yd1=-tgn3[0] / (t1*t1) * zgdif;
00695       yd2=-tgn3[1] / (t2*t2) * zgdif * (std::pow(((re+z2)/(re+z1)),2.0));
00696 
00697       /* calculate spline coefficients */
00698       spline (xs, ys, mn, yd1, yd2, y2out);
00699       x = zg/zgdif;
00700       splint (xs, ys, y2out, mn, x, &y);
00701 
00702       /* temperature at altitude */
00703       *tz = 1.0 / y;
00704       if (xm!=0.0) 
00705       {
00706          /* calaculate tropospheric / stratosphere density */
00707          glb = gsurf / (std::pow((1.0 + z1/re),2.0));
00708          gamm = xm * glb * zgdif / rgas;
00709 
00710          /* Integrate temperature profile */
00711          splini(xs, ys, y2out, mn, x, &yi);
00712          expl=gamm*yi;
00713          if (expl>50.0)
00714             expl=50.0;
00715 
00716          /* Density at altitude */
00717          densm_tmp = densm_tmp * (t1 / *tz) * std::exp(-expl);
00718       }
00719       if (xm==0.0)
00720          return *tz;
00721       else
00722          return densm_tmp;
00723 
00724    }  // 'densm'
00725 
00726 
00727    double Msise00Drag::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)
00728    {
00729       /*      Calculate Temperature and Density Profiles for MSIS models
00730       *      New lower thermo polynomial
00731       */
00732       double yd2, yd1, x, y;
00733       double rgas=831.4;
00734       double densu_temp=1.0;
00735       double za, z, zg2, tt, ta;
00736       double dta, z1, z2, t1, t2, zg, zgdif;
00737       int mn;
00738       int k;
00739       double glb;
00740       double expl;
00741       double yi;
00742       double densa;
00743       double gamma, gamm;
00744       double xs[5], ys[5], y2out[5];
00745       /* joining altitudes of Bates and spline */
00746       za=zn1[0];
00747       if (alt>za)
00748          z=alt;
00749       else
00750          z=za;
00751 
00752       /* geopotential altitude difference from ZLB */
00753       zg2 = zeta(z, zlb);
00754 
00755       /* Bates temperature */
00756       tt = tinf - (tinf - tlb) * std::exp(-s2*zg2);
00757       ta = tt;
00758       *tz = tt;
00759       densu_temp = *tz;
00760 
00761       if (alt<za) 
00762       {
00763          /* calculate temperature below ZA
00764          * temperature gradient at ZA from Bates profile */
00765          dta = (tinf - ta) * s2 * std::pow(((re+zlb)/(re+za)),2.0);
00766          tgn1[0]=dta;
00767          tn1[0]=ta;
00768          if (alt>zn1[mn1-1])
00769             z=alt;
00770          else
00771             z=zn1[mn1-1];
00772          mn=mn1;
00773          z1=zn1[0];
00774          z2=zn1[mn-1];
00775          t1=tn1[0];
00776          t2=tn1[mn-1];
00777          /* geopotental difference from z1 */
00778          zg = zeta (z, z1);
00779          zgdif = zeta(z2, z1);
00780          /* set up spline nodes */
00781          for (k=0;k<mn;k++) 
00782          {
00783             xs[k] = zeta(zn1[k], z1) / zgdif;
00784             ys[k] = 1.0 / tn1[k];
00785          }
00786          /* end node derivatives */
00787          yd1 = -tgn1[0] / (t1*t1) * zgdif;
00788          yd2 = -tgn1[1] / (t2*t2) * zgdif * std::pow(((re+z2)/(re+z1)),2.0);
00789          /* calculate spline coefficients */
00790          spline (xs, ys, mn, yd1, yd2, y2out);
00791          x = zg / zgdif;
00792          splint (xs, ys, y2out, mn, x, &y);
00793          /* temperature at altitude */
00794          *tz = 1.0 / y;
00795          densu_temp = *tz;
00796       }
00797       if (xm==0)
00798          return densu_temp;
00799 
00800       /* calculate density above za */
00801       glb = gsurf / std::pow((1.0 + zlb/re),2.0);
00802       gamma = xm * glb / (s2 * rgas * tinf);
00803       expl = std::exp(-s2 * gamma * zg2);
00804       if (expl>50.0)
00805          expl=50.0;
00806       if (tt<=0)
00807          expl=50.0;
00808 
00809       /* density at altitude */
00810       densa = dlb * std::pow((tlb/tt),((1.0+alpha+gamma))) * expl;
00811       densu_temp=densa;
00812       if (alt>=za)
00813          return densu_temp;
00814 
00815       /* calculate density below za */
00816       glb = gsurf / std::pow((1.0 + z1/re),2.0);
00817       gamm = xm * glb * zgdif / rgas;
00818 
00819       /* integrate spline temperatures */
00820       splini (xs, ys, y2out, mn, x, &yi);
00821       expl = gamm * yi;
00822       if (expl>50.0)
00823          expl=50.0;
00824       if (*tz<=0)
00825          expl=50.0;
00826 
00827       /* density at altitude */
00828       densu_temp = densu_temp * std::pow ((t1 / *tz),(1.0 + alpha)) * std::exp(-expl);
00829       return densu_temp;
00830 
00831 
00832    }  // 'densu()'
00833 
00834 
00835    double Msise00Drag::g0(double a, double *p) 
00836    {
00837       return (a - 4.0 + (p[25] - 1.0) * (a - 4.0 + (std::exp(-std::sqrt(p[24]*p[24]) * (a - 4.0)) - 1.0) / std::sqrt(p[24]*p[24])));
00838    }
00839 
00840    double Msise00Drag::sumex(double ex) 
00841    {
00842       return (1.0 + (1.0 - std::pow(ex,19.0)) / (1.0 - ex) * std::pow(ex,0.5));
00843    }
00844 
00845    double Msise00Drag::sg0(double ex, double *p, double *ap) 
00846    {
00847       return (g0(ap[1],p) + (g0(ap[2],p)*ex + g0(ap[3],p)*ex*ex + \
00848          g0(ap[4],p)*std::pow(ex,3.0)   + (g0(ap[5],p)*std::pow(ex,4.0) + \
00849          g0(ap[6],p)*std::pow(ex,12.0))*(1.0-std::pow(ex,8.0))/(1.0-ex)))/sumex(ex);
00850    }
00851 
00852    double Msise00Drag::globe7(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) 
00853    {
00854       /*       CALCULATE G(L) FUNCTION 
00855       *       Upper Thermosphere Parameters */
00856       double t[15];
00857       int i,j;
00858       int sw9=1;
00859       double apd;
00860       double xlong;
00861       double tloc;
00862       double c, s, c2, c4, s2;
00863       double sr = 7.2722E-5;
00864       double dgtr = 1.74533E-2;
00865       double dr = 1.72142E-2;
00866       double hr = 0.2618;
00867       double cd32, cd18, cd14, cd39;
00868       double p32, p18, p14, p39;
00869       double df;
00870       double f1, f2;
00871       double tinf;
00872       struct ap_array *ap;
00873 
00874       tloc=input->lst;
00875       for (j=0;j<14;j++)
00876          t[j]=0;
00877       if (flags->sw[9]>0)
00878          sw9=1;
00879       else if (flags->sw[9]<0)
00880          sw9=-1;
00881       xlong = input->g_long;
00882 
00883       /* calculate legendre polynomials */
00884       c = std::sin(input->g_lat * dgtr);
00885       s = std::cos(input->g_lat * dgtr);
00886       c2 = c*c;
00887       c4 = c2*c2;
00888       s2 = s*s;
00889 
00890       plg[0][1] = c;
00891       plg[0][2] = 0.5*(3.0*c2 -1.0);
00892       plg[0][3] = 0.5*(5.0*c*c2-3.0*c);
00893       plg[0][4] = (35.0*c4 - 30.0*c2 + 3.0)/8.0;
00894       plg[0][5] = (63.0*c2*c2*c - 70.0*c2*c + 15.0*c)/8.0;
00895       plg[0][6] = (11.0*c*plg[0][5] - 5.0*plg[0][4])/6.0;
00896       /*      plg[0][7] = (13.0*c*plg[0][6] - 6.0*plg[0][5])/7.0; */
00897       plg[1][1] = s;
00898       plg[1][2] = 3.0*c*s;
00899       plg[1][3] = 1.5*(5.0*c2-1.0)*s;
00900       plg[1][4] = 2.5*(7.0*c2*c-3.0*c)*s;
00901       plg[1][5] = 1.875*(21.0*c4 - 14.0*c2 +1.0)*s;
00902       plg[1][6] = (11.0*c*plg[1][5]-6.0*plg[1][4])/5.0;
00903       /*      plg[1][7] = (13.0*c*plg[1][6]-7.0*plg[1][5])/6.0; */
00904       /*      plg[1][8] = (15.0*c*plg[1][7]-8.0*plg[1][6])/7.0; */
00905       plg[2][2] = 3.0*s2;
00906       plg[2][3] = 15.0*s2*c;
00907       plg[2][4] = 7.5*(7.0*c2 -1.0)*s2;
00908       plg[2][5] = 3.0*c*plg[2][4]-2.0*plg[2][3];
00909       plg[2][6] =(11.0*c*plg[2][5]-7.0*plg[2][4])/4.0;
00910       plg[2][7] =(13.0*c*plg[2][6]-8.0*plg[2][5])/5.0;
00911       plg[3][3] = 15.0*s2*s;
00912       plg[3][4] = 105.0*s2*s*c; 
00913       plg[3][5] =(9.0*c*plg[3][4]-7.*plg[3][3])/2.0;
00914       plg[3][6] =(11.0*c*plg[3][5]-8.*plg[3][4])/3.0;
00915 
00916       if (!(((flags->sw[7]==0)&&(flags->sw[8]==0))&&(flags->sw[14]==0))) 
00917       {
00918          stloc = std::sin(hr*tloc);
00919          ctloc = std::cos(hr*tloc);
00920          s2tloc = std::sin(2.0*hr*tloc);
00921          c2tloc = std::cos(2.0*hr*tloc);
00922          s3tloc = std::sin(3.0*hr*tloc);
00923          c3tloc = std::cos(3.0*hr*tloc);
00924       }
00925 
00926       cd32 = std::cos(dr*(input->doy-p[31]));
00927       cd18 = std::cos(2.0*dr*(input->doy-p[17]));
00928       cd14 = std::cos(dr*(input->doy-p[13]));
00929       cd39 = std::cos(2.0*dr*(input->doy-p[38]));
00930       p32=p[31];
00931       p18=p[17];
00932       p14=p[13];
00933       p39=p[38];
00934 
00935       /* F10.7 EFFECT */
00936       df = input->f107 - input->f107A;
00937       dfa = input->f107A - 150.0;
00938       t[0] =  p[19]*df*(1.0+p[59]*dfa) + p[20]*df*df + p[21]*dfa + p[29]*std::pow(dfa,2.0);
00939       f1 = 1.0 + (p[47]*dfa +p[19]*df+p[20]*df*df)*flags->swc[1];
00940       f2 = 1.0 + (p[49]*dfa+p[19]*df+p[20]*df*df)*flags->swc[1];
00941 
00942       /*  TIME INDEPENDENT */
00943       t[1] = (p[1]*plg[0][2]+ p[2]*plg[0][4]+p[22]*plg[0][6]) + \
00944          (p[14]*plg[0][2])*dfa*flags->swc[1] +p[26]*plg[0][1];
00945 
00946       /*  SYMMETRICAL ANNUAL */
00947       t[2] = p[18]*cd32;
00948 
00949       /*  SYMMETRICAL SEMIANNUAL */
00950       t[3] = (p[15]+p[16]*plg[0][2])*cd18;
00951 
00952       /*  ASYMMETRICAL ANNUAL */
00953       t[4] =  f1*(p[9]*plg[0][1]+p[10]*plg[0][3])*cd14;
00954 
00955       /*  ASYMMETRICAL SEMIANNUAL */
00956       t[5] =    p[37]*plg[0][1]*cd39;
00957 
00958       /* DIURNAL */
00959       if (flags->sw[7]) 
00960       {
00961          double t71, t72;
00962          t71 = (p[11]*plg[1][2])*cd14*flags->swc[5];
00963          t72 = (p[12]*plg[1][2])*cd14*flags->swc[5];
00964          t[6] = f2*((p[3]*plg[1][1] + p[4]*plg[1][3] + p[27]*plg[1][5] + t71) * \
00965             ctloc + (p[6]*plg[1][1] + p[7]*plg[1][3] + p[28]*plg[1][5] \
00966             + t72)*stloc);
00967       }
00968 
00969       /* SEMIDIURNAL */
00970       if (flags->sw[8]) 
00971       {
00972          double t81, t82;
00973          t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
00974          t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
00975          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);
00976       }
00977 
00978       /* TERDIURNAL */
00979       if (flags->sw[14]) 
00980       {
00981          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);
00982       }
00983 
00984       /* magnetic activity based on daily ap */
00985       if (flags->sw[9]==-1) 
00986       {
00987          ap = input->ap_a;
00988          if (p[51]!=0) 
00989          {
00990             double exp1;
00991             exp1 = std::exp(-10800.0*std::sqrt(p[51]*p[51])/(1.0+p[138]*(45.0-std::sqrt(input->g_lat*input->g_lat))));
00992             if (exp1>0.99999)
00993                exp1=0.99999;
00994             if (p[24]<1.0E-4)
00995                p[24]=1.0E-4;
00996             apt[0]=sg0(exp1,p,ap->a);
00997             /* apt[1]=sg2(exp1,p,ap->a);
00998             apt[2]=sg0(exp2,p,ap->a);
00999             apt[3]=sg2(exp2,p,ap->a);
01000             */
01001             if (flags->sw[9]) 
01002             {
01003                t[8] = apt[0]*(p[50]+p[96]*plg[0][2]+p[54]*plg[0][4]+ \
01004                   (p[125]*plg[0][1]+p[126]*plg[0][3]+p[127]*plg[0][5])*cd14*flags->swc[5]+ \
01005                   (p[128]*plg[1][1]+p[129]*plg[1][3]+p[130]*plg[1][5])*flags->swc[7]* \
01006                   std::cos(hr*(tloc-p[131])));
01007             }
01008          }
01009       } 
01010       else 
01011       {
01012          double p44, p45;
01013          apd=input->ap-4.0;
01014          p44=p[43];
01015          p45=p[44];
01016          if (p44<0)
01017             p44 = 1.0E-5;
01018          apdf = apd + (p45-1.0)*(apd + (std::exp(-p44 * apd) - 1.0)/p44);
01019          if (flags->sw[9]) 
01020          {
01021             t[8]=apdf*(p[32]+p[45]*plg[0][2]+p[34]*plg[0][4]+ \
01022                (p[100]*plg[0][1]+p[101]*plg[0][3]+p[102]*plg[0][5])*cd14*flags->swc[5]+
01023                (p[121]*plg[1][1]+p[122]*plg[1][3]+p[123]*plg[1][5])*flags->swc[7]*
01024                std::cos(hr*(tloc-p[124])));
01025          }
01026       }
01027 
01028       if ((flags->sw[10])&&(input->g_long>-1000.0))
01029       {
01030 
01031          /* longitudinal */
01032          if (flags->sw[11]) 
01033          {
01034             t[10] = (1.0 + p[80]*dfa*flags->swc[1])* \
01035                ((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
01036                +p[103]*plg[1][1]+p[104]*plg[1][3]+p[105]*plg[1][5]\
01037                +flags->swc[5]*(p[109]*plg[1][1]+p[110]*plg[1][3]+p[111]*plg[1][5])*cd14)* \
01038                std::cos(dgtr*input->g_long) \
01039                +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
01040                +p[106]*plg[1][1]+p[107]*plg[1][3]+p[108]*plg[1][5]\
01041                +flags->swc[5]*(p[112]*plg[1][1]+p[113]*plg[1][3]+p[114]*plg[1][5])*cd14)* \
01042                std::sin(dgtr*input->g_long));
01043          }
01044 
01045          /* ut and mixed ut, longitude */
01046          if (flags->sw[12])
01047          {
01048             t[11]=(1.0+p[95]*plg[0][1])*(1.0+p[81]*dfa*flags->swc[1])*\
01049                (1.0+p[119]*plg[0][1]*flags->swc[5]*cd14)*\
01050                ((p[68]*plg[0][1]+p[69]*plg[0][3]+p[70]*plg[0][5])*\
01051                std::cos(sr*(input->sec-p[71])));
01052             t[11]+=flags->swc[11]*\
01053                (p[76]*plg[2][3]+p[77]*plg[2][5]+p[78]*plg[2][7])*\
01054                std::cos(sr*(input->sec-p[79])+2.0*dgtr*input->g_long)*(1.0+p[137]*dfa*flags->swc[1]);
01055          }
01056 
01057          /* ut, longitude magnetic activity */
01058          if (flags->sw[13])
01059          {
01060             if (flags->sw[9]==-1) 
01061             {
01062                if (p[51])
01063                {
01064                   t[12]=apt[0]*flags->swc[11]*(1.+p[132]*plg[0][1])*\
01065                      ((p[52]*plg[1][2]+p[98]*plg[1][4]+p[67]*plg[1][6])*\
01066                      std::cos(dgtr*(input->g_long-p[97])))\
01067                      +apt[0]*flags->swc[11]*flags->swc[5]*\
01068                      (p[133]*plg[1][1]+p[134]*plg[1][3]+p[135]*plg[1][5])*\
01069                      cd14*std::cos(dgtr*(input->g_long-p[136])) \
01070                      +apt[0]*flags->swc[12]* \
01071                      (p[55]*plg[0][1]+p[56]*plg[0][3]+p[57]*plg[0][5])*\
01072                      std::cos(sr*(input->sec-p[58]));
01073                }
01074             } 
01075             else 
01076             {
01077                t[12] = apdf*flags->swc[11]*(1.0+p[120]*plg[0][1])*\
01078                   ((p[60]*plg[1][2]+p[61]*plg[1][4]+p[62]*plg[1][6])*\
01079                   std::cos(dgtr*(input->g_long-p[63])))\
01080                   +apdf*flags->swc[11]*flags->swc[5]* \
01081                   (p[115]*plg[1][1]+p[116]*plg[1][3]+p[117]*plg[1][5])* \
01082                   cd14*std::cos(dgtr*(input->g_long-p[118])) \
01083                   + apdf*flags->swc[12]* \
01084                   (p[83]*plg[0][1]+p[84]*plg[0][3]+p[85]*plg[0][5])* \
01085                   std::cos(sr*(input->sec-p[75]));
01086             }         
01087          }
01088       }
01089 
01090       /* parms not used: 82, 89, 99, 139-149 */
01091       tinf = p[30];
01092       for (i=0;i<14;i++)
01093          tinf = tinf + std::abs(flags->sw[i+1])*t[i];
01094 
01095       return tinf;
01096 
01097    }  //'globe7()'
01098 
01099 
01100    double Msise00Drag::glob7s(double *p, struct nrlmsise_input *input, struct nrlmsise_flags *flags) 
01101    {
01102       /*    VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99 
01103       */
01104       double pset=2.0;
01105       double t[14];
01106       double tt;
01107       double cd32, cd18, cd14, cd39;
01108       double p32, p18, p14, p39;
01109       int i,j;
01110       double dr=1.72142E-2;
01111       double dgtr=1.74533E-2;
01112       /* confirm parameter set */
01113       if (p[99]==0)
01114          p[99]=pset;
01115       if (p[99]!=pset) 
01116       {
01117          printf("Wrong parameter set for glob7s\n");
01118          return -1;
01119       }
01120       for (j=0;j<14;j++)
01121          t[j]=0.0;
01122       cd32 = std::cos(dr*(input->doy-p[31]));
01123       cd18 = std::cos(2.0*dr*(input->doy-p[17]));
01124       cd14 = std::cos(dr*(input->doy-p[13]));
01125       cd39 = std::cos(2.0*dr*(input->doy-p[38]));
01126       p32=p[31];
01127       p18=p[17];
01128       p14=p[13];
01129       p39=p[38];
01130 
01131       /* F10.7 */
01132       t[0] = p[21]*dfa;
01133 
01134       /* time independent */
01135       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];
01136 
01137       /* SYMMETRICAL ANNUAL */
01138       t[2]=(p[18]+p[47]*plg[0][2]+p[29]*plg[0][4])*cd32;
01139 
01140       /* SYMMETRICAL SEMIANNUAL */
01141       t[3]=(p[15]+p[16]*plg[0][2]+p[30]*plg[0][4])*cd18;
01142 
01143       /* ASYMMETRICAL ANNUAL */
01144       t[4]=(p[9]*plg[0][1]+p[10]*plg[0][3]+p[20]*plg[0][5])*cd14;
01145 
01146       /* ASYMMETRICAL SEMIANNUAL */
01147       t[5]=(p[37]*plg[0][1])*cd39;
01148 
01149       /* DIURNAL */
01150       if (flags->sw[7])
01151       {
01152          double t71, t72;
01153          t71 = p[11]*plg[1][2]*cd14*flags->swc[5];
01154          t72 = p[12]*plg[1][2]*cd14*flags->swc[5];
01155          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) ;
01156       }
01157 
01158       /* SEMIDIURNAL */
01159       if (flags->sw[8]) 
01160       {
01161          double t81, t82;
01162          t81 = (p[23]*plg[2][3]+p[35]*plg[2][5])*cd14*flags->swc[5];
01163          t82 = (p[33]*plg[2][3]+p[36]*plg[2][5])*cd14*flags->swc[5];
01164          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);
01165       }
01166 
01167       /* TERDIURNAL */
01168       if (flags->sw[14])
01169       {
01170          t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
01171       }
01172 
01173       /* MAGNETIC ACTIVITY */
01174       if (flags->sw[9]) 
01175       {
01176          if (flags->sw[9]==1)
01177             t[8] = apdf * (p[32] + p[45] * plg[0][2] * flags->swc[2]);
01178          if (flags->sw[9]==-1)   
01179             t[8]=(p[50]*apt[0] + p[96]*plg[0][2] * apt[0]*flags->swc[2]);
01180       }
01181 
01182       /* LONGITUDINAL */
01183       if (!((flags->sw[10]==0) || (flags->sw[11]==0) || (input->g_long<=-1000.0)))
01184       {
01185          t[10] = (1.0 + plg[0][1]*(p[80]*flags->swc[5]*std::cos(dr*(input->doy-p[81]))\
01186             +p[85]*flags->swc[6]*std::cos(2.0*dr*(input->doy-p[86])))\
01187             +p[83]*flags->swc[3]*std::cos(dr*(input->doy-p[84]))\
01188             +p[87]*flags->swc[4]*std::cos(2.0*dr*(input->doy-p[88])))\
01189             *((p[64]*plg[1][2]+p[65]*plg[1][4]+p[66]*plg[1][6]\
01190             +p[74]*plg[1][1]+p[75]*plg[1][3]+p[76]*plg[1][5]\
01191             )*std::cos(dgtr*input->g_long)\
01192             +(p[90]*plg[1][2]+p[91]*plg[1][4]+p[92]*plg[1][6]\
01193             +p[77]*plg[1][1]+p[78]*plg[1][3]+p[79]*plg[1][5]\
01194             )*std::sin(dgtr*input->g_long));
01195       }
01196       tt=0;
01197       for (i=0;i<14;i++)
01198          tt+=std::abs(flags->sw[i+1])*t[i];
01199       return tt;
01200 
01201    }  // 'glob7s()'
01202 
01203 
01204    void Msise00Drag::gtd7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output)
01205    {
01206       double xlat;
01207       double xmm;
01208       int mn3 = 5;
01209       double zn3[5]={32.5,20.0,15.0,10.0,0.0};
01210       int mn2 = 4;
01211       double zn2[4]={72.5,55.0,45.0,32.5};
01212       double altt;
01213       double zmix=62.5;
01214       double tmp;
01215       double dm28m;
01216       double tz;
01217       double dmc;
01218       double dmr;
01219       double dz28;
01220       struct nrlmsise_output soutput;
01221       int i;
01222 
01223       tselec(flags);
01224 
01225       /* Latitude variation of gravity (none for sw[2]=0) */
01226       xlat=input->g_lat;
01227       if (flags->sw[2]==0)
01228          xlat=45.0;
01229       glatf(xlat, &gsurf, &re);
01230 
01231       xmm = pdm[2][4];
01232 
01233       /* THERMOSPHERE / MESOSPHERE (above zn2[0]) */
01234       if (input->alt>zn2[0])
01235          altt=input->alt;
01236       else
01237          altt=zn2[0];
01238 
01239       tmp=input->alt;
01240       input->alt=altt;
01241       gts7(input, flags, &soutput);
01242       altt=input->alt;
01243       input->alt=tmp;
01244       if (flags->sw[0])   /* metric adjustment */
01245          dm28m=dm28*1.0E6;
01246       else
01247          dm28m=dm28;
01248       output->t[0]=soutput.t[0];
01249       output->t[1]=soutput.t[1];
01250       if (input->alt>=zn2[0]) 
01251       {
01252          for (i=0;i<9;i++)
01253             output->d[i]=soutput.d[i];
01254          return;
01255       }
01256 
01257       /*       LOWER MESOSPHERE/UPPER STRATOSPHERE (between zn3[0] and zn2[0])
01258       *         Temperature at nodes and gradients at end nodes
01259       *         Inverse temperature a linear function of spherical harmonics
01260       */
01261       meso_tgn2[0]=meso_tgn1[1];
01262       meso_tn2[0]=meso_tn1[4];
01263       meso_tn2[1]=pma[0][0]*pavgm[0]/(1.0-flags->sw[20]*glob7s(pma[0], input, flags));
01264       meso_tn2[2]=pma[1][0]*pavgm[1]/(1.0-flags->sw[20]*glob7s(pma[1], input, flags));
01265       meso_tn2[3]=pma[2][0]*pavgm[2]/(1.0-flags->sw[20]*flags->sw[22]*glob7s(pma[2], input, flags));
01266       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]/(std::pow((pma[2][0]*pavgm[2]),2.0));
01267       meso_tn3[0]=meso_tn2[3];
01268 
01269       if (input->alt<zn3[0])
01270       {
01271          /*       LOWER STRATOSPHERE AND TROPOSPHERE (below zn3[0])
01272          *         Temperature at nodes and gradients at end nodes
01273          *         Inverse temperature a linear function of spherical harmonics
01274          */
01275          meso_tgn3[0]=meso_tgn2[1];
01276          meso_tn3[1]=pma[3][0]*pavgm[3]/(1.0-flags->sw[22]*glob7s(pma[3], input, flags));
01277          meso_tn3[2]=pma[4][0]*pavgm[4]/(1.0-flags->sw[22]*glob7s(pma[4], input, flags));
01278          meso_tn3[3]=pma[5][0]*pavgm[5]/(1.0-flags->sw[22]*glob7s(pma[5], input, flags));
01279          meso_tn3[4]=pma[6][0]*pavgm[6]/(1.0-flags->sw[22]*glob7s(pma[6], input, flags));
01280          meso_tgn3[1]=pma[7][0]*pavgm[7]*(1.0+flags->sw[22]*glob7s(pma[7], input, flags)) *meso_tn3[4]*meso_tn3[4]/(std::pow((pma[6][0]*pavgm[6]),2.0));
01281       }
01282 
01283       /* LINEAR TRANSITION TO FULL MIXING BELOW zn2[0] */
01284 
01285       dmc=0;
01286       if (input->alt>zmix)
01287          dmc = 1.0 - (zn2[0]-input->alt)/(zn2[0] - zmix);
01288       dz28=soutput.d[2];
01289 
01290       /**** N2 density ****/
01291       dmr=soutput.d[2] / dm28m - 1.0;
01292       output->d[2]=densm(input->alt,dm28m,xmm, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
01293       output->d[2]=output->d[2] * (1.0 + dmr*dmc);
01294 
01295       /**** HE density ****/
01296       dmr = soutput.d[0] / (dz28 * pdm[0][1]) - 1.0;
01297       output->d[0] = output->d[2] * pdm[0][1] * (1.0 + dmr*dmc);
01298 
01299       /**** O density ****/
01300       output->d[1] = 0;
01301       output->d[8] = 0;
01302 
01303       /**** O2 density ****/
01304       dmr = soutput.d[3] / (dz28 * pdm[3][1]) - 1.0;
01305       output->d[3] = output->d[2] * pdm[3][1] * (1.0 + dmr*dmc);
01306 
01307       /**** AR density ***/
01308       dmr = soutput.d[4] / (dz28 * pdm[4][1]) - 1.0;
01309       output->d[4] = output->d[2] * pdm[4][1] * (1.0 + dmr*dmc);
01310 
01311       /**** Hydrogen density ****/
01312       output->d[6] = 0;
01313 
01314       /**** Atomic nitrogen density ****/
01315       output->d[7] = 0;
01316 
01317       /**** Total mass density */
01318       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]);
01319 
01320       if (flags->sw[0])
01321          output->d[5]=output->d[5]/1000;
01322 
01323       /**** temperature at altitude ****/
01324       dd = densm(input->alt, 1.0, 0, &tz, mn3, zn3, meso_tn3, meso_tgn3, mn2, zn2, meso_tn2, meso_tgn2);
01325       output->t[1]=tz;
01326 
01327    }  // 'gtd7()'
01328 
01329 
01330    void Msise00Drag::gtd7d(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) 
01331    {
01332       gtd7(input, flags, output);
01333       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]);
01334       if (flags->sw[0])
01335          output->d[5]=output->d[5]/1000;
01336 
01337    }  // 'gtd7d()'
01338 
01339 
01340 
01341    void Msise00Drag::ghp7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output, double press) 
01342    {
01343       double bm = 1.3806E-19;
01344       double rgas = 831.4;
01345       double test = 0.00043;
01346       double ltest = 12;
01347       double pl, p;
01348       double zi;
01349       double z;
01350       double cl, cl2;
01351       double ca, cd;
01352       double xn, xm, diff;
01353       double g, sh;
01354       int l;
01355       pl = std::log10(press);
01356       if (pl >= -5.0) 
01357       {
01358          if (pl>2.5)
01359             zi = 18.06 * (3.00 - pl);
01360          else if ((pl>0.075) && (pl<=2.5))
01361             zi = 14.98 * (3.08 - pl);
01362          else if ((pl>-1) && (pl<=0.075))
01363             zi = 17.80 * (2.72 - pl);
01364          else if ((pl>-2) && (pl<=-1))
01365             zi = 14.28 * (3.64 - pl);
01366          else if ((pl>-4) && (pl<=-2))
01367             zi = 12.72 * (4.32 -pl);
01368          else if (pl<=-4)
01369             zi = 25.3 * (0.11 - pl);
01370          cl = input->g_lat/90.0;
01371          cl2 = cl*cl;
01372          if (input->doy<182)
01373             cd = (1.0 - (double) input->doy) / 91.25;
01374          else 
01375             cd = ((double) input->doy) / 91.25 - 3.0;
01376          ca = 0;
01377          if ((pl > -1.11) && (pl<=-0.23))
01378             ca = 1.0;
01379          if (pl > -0.23)
01380             ca = (2.79 - pl) / (2.79 + 0.23);
01381          if ((pl <= -1.11) && (pl>-3))
01382             ca = (-2.93 - pl)/(-2.93 + 1.11);
01383          z = zi - 4.87 * cl * cd * ca - 1.64 * cl2 * ca + 0.31 * ca * cl;
01384       } else
01385          z = 22.0 * std::pow((pl + 4.0),2.0) + 110.0;
01386 
01387       /* iteration  loop */
01388       l = 0;
01389       do 
01390       {
01391          l++;
01392          input->alt = z;
01393          gtd7(input, flags, output);
01394          z = input->alt;
01395          xn = output->d[0] + output->d[1] + output->d[2] + output->d[3] + output->d[4] + output->d[6] + output->d[7];
01396          p = bm * xn * output->t[1];
01397          if (flags->sw[0])
01398             p = p*1.0E-6;
01399          diff = pl - std::log10(p);
01400          if (std::sqrt(diff*diff)<test)
01401             return;
01402          if (l==ltest) {
01403             printf("ERROR: ghp7 not converging for press %e, diff %e",press,diff);
01404             return;
01405          }
01406          xm = output->d[5] / xn / 1.66E-24;
01407          if (flags->sw[0])
01408             xm = xm * 1.0E3;
01409          g = gsurf / (std::pow((1.0 + z/re),2.0));
01410          sh = rgas * output->t[1] / (xm * g);
01411 
01412          /* new altitude estimate using scale height */
01413          if (l <  6)
01414             z = z - sh * diff * 2.302;
01415          else
01416             z = z - sh * diff;
01417       } while (1==1);
01418 
01419    }  // 'ghp7()'
01420 
01421 
01422    void Msise00Drag::gts7(struct nrlmsise_input *input, struct nrlmsise_flags *flags, struct nrlmsise_output *output) 
01423    {
01424       /*     Thermospheric portion of NRLMSISE-00
01425       *     See GTD7 for more extensive comments
01426       *     alt > 72.5 km! 
01427       */
01428       double za;
01429       int i, j;
01430       double ddum, z;
01431       double zn1[5] = {120.0, 110.0, 100.0, 90.0, 72.5};
01432       double tinf;
01433       int mn1 = 5;
01434       double g0;
01435       double tlb;
01436       double s, z0, t0, tr12;
01437       double db01, db04, db14, db16, db28, db32, db40, db48;
01438       double zh28, zh04, zh16, zh32, zh40, zh01, zh14;
01439       double zhm28, zhm04, zhm16, zhm32, zhm40, zhm01, zhm14;
01440       double xmd;
01441       double b28, b04, b16, b32, b40, b01, b14;
01442       double tz;
01443       double g28, g4, g16, g32, g40, g1, g14;
01444       double zhf, xmm;
01445       double zc04, zc16, zc32, zc40, zc01, zc14;
01446       double hc04, hc16, hc32, hc40, hc01, hc14;
01447       double hcc16, hcc32, hcc01, hcc14;
01448       double zcc16, zcc32, zcc01, zcc14;
01449       double rc16, rc32, rc01, rc14;
01450       double rl;
01451       double g16h, db16h, tho, zsht, zmho, zsho;
01452       double dgtr=1.74533E-2;
01453       double dr=1.72142E-2;
01454       double alpha[9]={-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
01455       double altl[8]={200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
01456       double dd;
01457       double hc216, hcc232;
01458       za = pdl[1][15];
01459       zn1[0] = za;
01460       for (j=0;j<9;j++) 
01461          output->d[j]=0;
01462 
01463       /* TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1) */
01464       if (input->alt>zn1[0])
01465          tinf = ptm[0]*pt[0] * \
01466          (1.0+flags->sw[16]*globe7(pt,input,flags));
01467       else
01468          tinf = ptm[0]*pt[0];
01469       output->t[0]=tinf;
01470 
01471       /*  GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5) */
01472       if (input->alt>zn1[4])
01473          g0 = ptm[3]*ps[0] * \
01474          (1.0+flags->sw[19]*globe7(ps,input,flags));
01475       else
01476          g0 = ptm[3]*ps[0];
01477       tlb = ptm[1] * (1.0 + flags->sw[17]*globe7(pd[3],input,flags))*pd[3][0];
01478       s = g0 / (tinf - tlb);
01479 
01480       /*      Lower thermosphere temp variations not significant for
01481       *       density above 300 km */
01482       if (input->alt<300.0) 
01483       {
01484          meso_tn1[1]=ptm[6]*ptl[0][0]/(1.0-flags->sw[18]*glob7s(ptl[0], input, flags));
01485          meso_tn1[2]=ptm[2]*ptl[1][0]/(1.0-flags->sw[18]*glob7s(ptl[1], input, flags));
01486          meso_tn1[3]=ptm[7]*ptl[2][0]/(1.0-flags->sw[18]*glob7s(ptl[2], input, flags));
01487          meso_tn1[4]=ptm[4]*ptl[3][0]/(1.0-flags->sw[18]*flags->sw[20]*glob7s(ptl[3], input, flags));
01488          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]/(std::pow((ptm[4]*ptl[3][0]),2.0));
01489       } else
01490       {
01491          meso_tn1[1]=ptm[6]*ptl[0][0];
01492          meso_tn1[2]=ptm[2]*ptl[1][0];
01493          meso_tn1[3]=ptm[7]*ptl[2][0];
01494          meso_tn1[4]=ptm[4]*ptl[3][0];
01495          meso_tgn1[1]=ptm[8]*pma[8][0]*meso_tn1[4]*meso_tn1[4]/(std::pow((ptm[4]*ptl[3][0]),2.0));
01496       }
01497 
01498       z0 = zn1[3];
01499       t0 = meso_tn1[3];
01500       tr12 = 1.0;
01501 
01502       /* N2 variation factor at Zlb */
01503       g28=flags->sw[21]*globe7(pd[2], input, flags);
01504 
01505       /* VARIATION OF TURBOPAUSE HEIGHT */
01506       zhf=pdl[1][24]*(1.0+flags->sw[5]*pdl[0][24]*std::sin(dgtr*input->g_lat)*std::cos(dr*(input->doy-pt[13])));
01507       output->t[0]=tinf;
01508       xmm = pdm[2][4];
01509       z = input->alt;
01510 
01511 
01512       /**** N2 DENSITY ****/
01513 
01514       /* Diffusive density at Zlb */
01515       db28 = pdm[2][0]*std::exp(g28)*pd[2][0];
01516       /* Diffusive density at Alt */
01517       output->d[2]=densu(z,db28,tinf,tlb,28.0,alpha[2],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01518       dd=output->d[2];
01519       /* Turbopause */
01520       zh28=pdm[2][2]*zhf;
01521       zhm28=pdm[2][3]*pdl[1][5]; 
01522       xmd=28.0-xmm;
01523       /* Mixed density at Zlb */
01524       b28=densu(zh28,db28,tinf,tlb,xmd,(alpha[2]-1.0),&tz,ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
01525       if ((flags->sw[15])&&(z<=altl[2])) 
01526       {
01527          /*  Mixed density at Alt */
01528          dm28=densu(z,b28,tinf,tlb,xmm,alpha[2],&tz,ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01529          /*  Net density at Alt */
01530          output->d[2]=dnet(output->d[2],dm28,zhm28,xmm,28.0);
01531       }
01532 
01533 
01534       /**** HE DENSITY ****/
01535 
01536       /*   Density variation factor at Zlb */
01537       g4 = flags->sw[21]*globe7(pd[0], input, flags);
01538       /*  Diffusive density at Zlb */
01539       db04 = pdm[0][0]*std::exp(g4)*pd[0][0];
01540       /*  Diffusive density at Alt */
01541       output->d[0]=densu(z,db04,tinf,tlb, 4.,alpha[0],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01542       dd=output->d[0];
01543       if ((flags->sw[15]) && (z<altl[0]))
01544       {
01545          /*  Turbopause */
01546          zh04=pdm[0][2];
01547          /*  Mixed density at Zlb */
01548          b04=densu(zh04,db04,tinf,tlb,4.-xmm,alpha[0]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01549          /*  Mixed density at Alt */
01550          dm04=densu(z,b04,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01551          zhm04=zhm28;
01552          /*  Net density at Alt */
01553          output->d[0]=dnet(output->d[0],dm04,zhm04,xmm,4.);
01554          /*  Correction to specified mixing ratio at ground */
01555          rl=std::log(b28*pdm[0][1]/b04);
01556          zc04=pdm[0][4]*pdl[1][0];
01557          hc04=pdm[0][5]*pdl[1][1];
01558          /*  Net density corrected at Alt */
01559          output->d[0]=output->d[0]*ccor(z,rl,hc04,zc04);
01560       }
01561 
01562 
01563       /**** O DENSITY ****/
01564 
01565       /*  Density variation factor at Zlb */
01566       g16= flags->sw[21]*globe7(pd[1],input,flags);
01567       /*  Diffusive density at Zlb */
01568       db16 =  pdm[1][0]*std::exp(g16)*pd[1][0];
01569       /*   Diffusive density at Alt */
01570       output->d[1]=densu(z,db16,tinf,tlb, 16.,alpha[1],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
01571       dd=output->d[1];
01572       if ((flags->sw[15]) && (z<=altl[1])) 
01573       {
01574          /*   Turbopause */
01575          zh16=pdm[1][2];
01576          /*  Mixed density at Zlb */
01577          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);
01578          /*  Mixed density at Alt */
01579          dm16=densu(z,b16,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01580          zhm16=zhm28;
01581          /*  Net density at Alt */
01582          output->d[1]=dnet(output->d[1],dm16,zhm16,xmm,16.);
01583          rl=pdm[1][1]*pdl[1][16]*(1.0+flags->sw[1]*pdl[0][23]*(input->f107A-150.0));
01584          hc16=pdm[1][5]*pdl[1][3];
01585          zc16=pdm[1][4]*pdl[1][2];
01586          hc216=pdm[1][5]*pdl[1][4];
01587          output->d[1]=output->d[1]*ccor2(z,rl,hc16,zc16,hc216);
01588          /*   Chemistry correction */
01589          hcc16=pdm[1][7]*pdl[1][13];
01590          zcc16=pdm[1][6]*pdl[1][12];
01591          rc16=pdm[1][3]*pdl[1][14];
01592          /*  Net density corrected at Alt */
01593          output->d[1]=output->d[1]*ccor(z,rc16,hcc16,zcc16);
01594       }
01595 
01596 
01597       /**** O2 DENSITY ****/
01598 
01599       /*   Density variation factor at Zlb */
01600       g32= flags->sw[21]*globe7(pd[4], input, flags);
01601       /*  Diffusive density at Zlb */
01602       db32 = pdm[3][0]*std::exp(g32)*pd[4][0];
01603       /*   Diffusive density at Alt */
01604       output->d[3]=densu(z,db32,tinf,tlb, 32.,alpha[3],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
01605       dd=output->d[3];
01606       if (flags->sw[15]) 
01607       {
01608          if (z<=altl[3]) 
01609          {
01610             /*   Turbopause */
01611             zh32=pdm[3][2];
01612             /*  Mixed density at Zlb */
01613             b32=densu(zh32,db32,tinf,tlb,32.-xmm,alpha[3]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01614             /*  Mixed density at Alt */
01615             dm32=densu(z,b32,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01616             zhm32=zhm28;
01617             /*  Net density at Alt */
01618             output->d[3]=dnet(output->d[3],dm32,zhm32,xmm,32.);
01619             /*   Correction to specified mixing ratio at ground */
01620             rl=std::log(b28*pdm[3][1]/b32);
01621             hc32=pdm[3][5]*pdl[1][7];
01622             zc32=pdm[3][4]*pdl[1][6];
01623             output->d[3]=output->d[3]*ccor(z,rl,hc32,zc32);
01624          }
01625          /*  Correction for general departure from diffusive equilibrium above Zlb */
01626          hcc32=pdm[3][7]*pdl[1][22];
01627          hcc232=pdm[3][7]*pdl[0][22];
01628          zcc32=pdm[3][6]*pdl[1][21];
01629          rc32=pdm[3][3]*pdl[1][23]*(1.+flags->sw[1]*pdl[0][23]*(input->f107A-150.));
01630          /*  Net density corrected at Alt */
01631          output->d[3]=output->d[3]*ccor2(z,rc32,hcc32,zcc32,hcc232);
01632       }
01633 
01634 
01635       /**** AR DENSITY ****/
01636 
01637       /*   Density variation factor at Zlb */
01638       g40= flags->sw[20]*globe7(pd[5],input,flags);
01639       /*  Diffusive density at Zlb */
01640       db40 = pdm[4][0]*std::exp(g40)*pd[5][0];
01641       /*   Diffusive density at Alt */
01642       output->d[4]=densu(z,db40,tinf,tlb, 40.,alpha[4],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01643       dd=output->d[4];
01644       if ((flags->sw[15]) && (z<=altl[4]))
01645       {
01646          /*   Turbopause */
01647          zh40=pdm[4][2];
01648          /*  Mixed density at Zlb */
01649          b40=densu(zh40,db40,tinf,tlb,40.-xmm,alpha[4]-1.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01650          /*  Mixed density at Alt */
01651          dm40=densu(z,b40,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01652          zhm40=zhm28;
01653          /*  Net density at Alt */
01654          output->d[4]=dnet(output->d[4],dm40,zhm40,xmm,40.);
01655          /*   Correction to specified mixing ratio at ground */
01656          rl=std::log(b28*pdm[4][1]/b40);
01657          hc40=pdm[4][5]*pdl[1][9];
01658          zc40=pdm[4][4]*pdl[1][8];
01659          /*  Net density corrected at Alt */
01660          output->d[4]=output->d[4]*ccor(z,rl,hc40,zc40);
01661       }
01662 
01663 
01664       /**** HYDROGEN DENSITY ****/
01665 
01666       /*   Density variation factor at Zlb */
01667       g1 = flags->sw[21]*globe7(pd[6], input, flags);
01668       /*  Diffusive density at Zlb */
01669       db01 = pdm[5][0]*std::exp(g1)*pd[6][0];
01670       /*   Diffusive density at Alt */
01671       output->d[6]=densu(z,db01,tinf,tlb,1.,alpha[6],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01672       dd=output->d[6];
01673       if ((flags->sw[15]) && (z<=altl[6])) 
01674       {
01675          /*   Turbopause */
01676          zh01=pdm[5][2];
01677          /*  Mixed density at Zlb */
01678          b01=densu(zh01,db01,tinf,tlb,1.-xmm,alpha[6]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01679          /*  Mixed density at Alt */
01680          dm01=densu(z,b01,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01681          zhm01=zhm28;
01682          /*  Net density at Alt */
01683          output->d[6]=dnet(output->d[6],dm01,zhm01,xmm,1.);
01684          /*   Correction to specified mixing ratio at ground */
01685          rl=std::log(b28*pdm[5][1]*std::sqrt(pdl[1][17]*pdl[1][17])/b01);
01686          hc01=pdm[5][5]*pdl[1][11];
01687          zc01=pdm[5][4]*pdl[1][10];
01688          output->d[6]=output->d[6]*ccor(z,rl,hc01,zc01);
01689          /*   Chemistry correction */
01690          hcc01=pdm[5][7]*pdl[1][19];
01691          zcc01=pdm[5][6]*pdl[1][18];
01692          rc01=pdm[5][3]*pdl[1][20];
01693          /*  Net density corrected at Alt */
01694          output->d[6]=output->d[6]*ccor(z,rc01,hcc01,zcc01);
01695       }
01696 
01697 
01698       /**** ATOMIC NITROGEN DENSITY ****/
01699 
01700       /*   Density variation factor at Zlb */
01701       g14 = flags->sw[21]*globe7(pd[7],input,flags);
01702       /*  Diffusive density at Zlb */
01703       db14 = pdm[6][0]*std::exp(g14)*pd[7][0];
01704       /*   Diffusive density at Alt */
01705       output->d[7]=densu(z,db14,tinf,tlb,14.,alpha[7],&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01706       dd=output->d[7];
01707       if ((flags->sw[15]) && (z<=altl[7])) 
01708       {
01709          /*   Turbopause */
01710          zh14=pdm[6][2];
01711          /*  Mixed density at Zlb */
01712          b14=densu(zh14,db14,tinf,tlb,14.-xmm,alpha[7]-1., &output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01713          /*  Mixed density at Alt */
01714          dm14=densu(z,b14,tinf,tlb,xmm,0.,&output->t[1],ptm[5],s,mn1,zn1,meso_tn1,meso_tgn1);
01715          zhm14=zhm28;
01716          /*  Net density at Alt */
01717          output->d[7]=dnet(output->d[7],dm14,zhm14,xmm,14.);
01718          /*   Correction to specified mixing ratio at ground */
01719          rl=std::log(b28*pdm[6][1]*std::sqrt(pdl[0][2]*pdl[0][2])/b14);
01720          hc14=pdm[6][5]*pdl[0][1];
01721          zc14=pdm[6][4]*pdl[0][0];
01722          output->d[7]=output->d[7]*ccor(z,rl,hc14,zc14);
01723          /*   Chemistry correction */
01724          hcc14=pdm[6][7]*pdl[0][4];
01725          zcc14=pdm[6][6]*pdl[0][3];
01726          rc14=pdm[6][3]*pdl[0][5];
01727          /*  Net density corrected at Alt */
01728          output->d[7]=output->d[7]*ccor(z,rc14,hcc14,zcc14);
01729       }
01730 
01731 
01732       /**** Anomalous OXYGEN DENSITY ****/
01733 
01734       g16h = flags->sw[21]*globe7(pd[8],input,flags);
01735       db16h = pdm[7][0]*std::exp(g16h)*pd[8][0];
01736       tho = pdm[7][9]*pdl[0][6];
01737       dd=densu(z,db16h,tho,tho,16.,alpha[8],&output->t[1],ptm[5],s,mn1, zn1,meso_tn1,meso_tgn1);
01738       zsht=pdm[7][5];
01739       zmho=pdm[7][4];
01740       zsho=scalh(zmho,16.0,tho);
01741       output->d[8]=dd*std::exp(-zsht/zsho*(std::exp(-(z-zmho)/zsht)-1.));
01742 
01743 
01744       /* total mass density */
01745       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]);
01746       db48=1.66E-24*(4.0*db04+16.0*db16+28.0*db28+32.0*db32+40.0*db40+db01+14.0*db14);
01747 
01748 
01749 
01750       /* temperature */
01751       z = std::sqrt(input->alt*input->alt);
01752       ddum = densu(z,1.0, tinf, tlb, 0.0, 0.0, &output->t[1], ptm[5], s, mn1, zn1, meso_tn1, meso_tgn1);
01753       if (flags->sw[0]) 
01754       {
01755          for(i=0;i<9;i++)
01756             output->d[i]=output->d[i]*1.0E6;
01757          output->d[5]=output->d[5]/1000;
01758       }
01759 
01760    }  //'gts7'
01761 
01762 
01763 
01764 
01766 
01767    /* -------------------------------------------------------------------- */
01768    /* ---------  N R L M S I S E - 0 0    M O D E L    2 0 0 1  ---------- */
01769    /* -------------------------------------------------------------------- */
01770 
01771    /* This file is part of the NRLMSISE-00  C source code package - release
01772     * 20041227
01773     *
01774     * The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and
01775     * Doug Drob. They also wrote a NRLMSISE-00 distribution package in 
01776     * FORTRAN which is available at
01777     * http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm
01778     *
01779     * Dominik Brodowski implemented and maintains this C version. You can
01780     * reach him at mail@brodo.de. See the file "DOCUMENTATION" for details,
01781     * and check http://www.brodo.de/english/pub/nrlmsise/index.html for
01782     * updated releases of this package.
01783     */
01784 
01785    /* TEMPERATURE */
01786    double Msise00Drag::pt[150] = {
01787        9.86573E-01, 1.62228E-02, 1.55270E-02,-1.04323E-01,-3.75801E-03,
01788       -1.18538E-03,-1.24043E-01, 4.56820E-03, 8.76018E-03,-1.36235E-01,
01789       -3.52427E-02, 8.84181E-03,-5.92127E-03,-8.61650E+00, 0.00000E+00,
01790        1.28492E-02, 0.00000E+00, 1.30096E+02, 1.04567E-02, 1.65686E-03,
01791       -5.53887E-06, 2.97810E-03, 0.00000E+00, 5.13122E-03, 8.66784E-02,
01792        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.27026E-06,
01793        0.00000E+00, 6.74494E+00, 4.93933E-03, 2.21656E-03, 2.50802E-03,
01794        0.00000E+00, 0.00000E+00,-2.08841E-02,-1.79873E+00, 1.45103E-03,
01795        2.81769E-04,-1.44703E-03,-5.16394E-05, 8.47001E-02, 1.70147E-01,
01796        5.72562E-03, 5.07493E-05, 4.36148E-03, 1.17863E-04, 4.74364E-03,
01797        6.61278E-03, 4.34292E-05, 1.44373E-03, 2.41470E-05, 2.84426E-03,
01798        8.56560E-04, 2.04028E-03, 0.00000E+00,-3.15994E+03,-2.46423E-03,
01799        1.13843E-03, 4.20512E-04, 0.00000E+00,-9.77214E+01, 6.77794E-03,
01800        5.27499E-03, 1.14936E-03, 0.00000E+00,-6.61311E-03,-1.84255E-02,
01801       -1.96259E-02, 2.98618E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01802        6.44574E+02, 8.84668E-04, 5.05066E-04, 0.00000E+00, 4.02881E+03,
01803       -1.89503E-03, 0.00000E+00, 0.00000E+00, 8.21407E-04, 2.06780E-03,
01804        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01805       -1.20410E-02,-3.63963E-03, 9.92070E-05,-1.15284E-04,-6.33059E-05,
01806       -6.05545E-01, 8.34218E-03,-9.13036E+01, 3.71042E-04, 0.00000E+00,
01807        4.19000E-04, 2.70928E-03, 3.31507E-03,-4.44508E-03,-4.96334E-03,
01808       -1.60449E-03, 3.95119E-03, 2.48924E-03, 5.09815E-04, 4.05302E-03,
01809        2.24076E-03, 0.00000E+00, 6.84256E-03, 4.66354E-04, 0.00000E+00,
01810       -3.68328E-04, 0.00000E+00, 0.00000E+00,-1.46870E+02, 0.00000E+00,
01811        0.00000E+00, 1.09501E-03, 4.65156E-04, 5.62583E-04, 3.21596E+00,
01812        6.43168E-04, 3.14860E-03, 3.40738E-03, 1.78481E-03, 9.62532E-04,
01813        5.58171E-04, 3.43731E+00,-2.33195E-01, 5.10289E-04, 0.00000E+00,
01814        0.00000E+00,-9.25347E+04, 0.00000E+00,-1.99639E-03, 0.00000E+00,
01815        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01816        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
01817    };
01818 
01819    double Msise00Drag::pd[9][150] = {
01820    /*  HE DENSITY */ {
01821        1.09979E+00,-4.88060E-02,-1.97501E-01,-9.10280E-02,-6.96558E-03,
01822        2.42136E-02, 3.91333E-01,-7.20068E-03,-3.22718E-02, 1.41508E+00,
01823        1.68194E-01, 1.85282E-02, 1.09384E-01,-7.24282E+00, 0.00000E+00,
01824        2.96377E-01,-4.97210E-02, 1.04114E+02,-8.61108E-02,-7.29177E-04,
01825        1.48998E-06, 1.08629E-03, 0.00000E+00, 0.00000E+00, 8.31090E-02,
01826        1.12818E-01,-5.75005E-02,-1.29919E-02,-1.78849E-02,-2.86343E-06,
01827        0.00000E+00,-1.51187E+02,-6.65902E-03, 0.00000E+00,-2.02069E-03,
01828        0.00000E+00, 0.00000E+00, 4.32264E-02,-2.80444E+01,-3.26789E-03,
01829        2.47461E-03, 0.00000E+00, 0.00000E+00, 9.82100E-02, 1.22714E-01,
01830       -3.96450E-02, 0.00000E+00,-2.76489E-03, 0.00000E+00, 1.87723E-03,
01831       -8.09813E-03, 4.34428E-05,-7.70932E-03, 0.00000E+00,-2.28894E-03,
01832       -5.69070E-03,-5.22193E-03, 6.00692E-03,-7.80434E+03,-3.48336E-03,
01833       -6.38362E-03,-1.82190E-03, 0.00000E+00,-7.58976E+01,-2.17875E-02,
01834       -1.72524E-02,-9.06287E-03, 0.00000E+00, 2.44725E-02, 8.66040E-02,
01835        1.05712E-01, 3.02543E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01836       -6.01364E+03,-5.64668E-03,-2.54157E-03, 0.00000E+00, 3.15611E+02,
01837       -5.69158E-03, 0.00000E+00, 0.00000E+00,-4.47216E-03,-4.49523E-03,
01838        4.64428E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01839        4.51236E-02, 2.46520E-02, 6.17794E-03, 0.00000E+00, 0.00000E+00,
01840       -3.62944E-01,-4.80022E-02,-7.57230E+01,-1.99656E-03, 0.00000E+00,
01841       -5.18780E-03,-1.73990E-02,-9.03485E-03, 7.48465E-03, 1.53267E-02,
01842        1.06296E-02, 1.18655E-02, 2.55569E-03, 1.69020E-03, 3.51936E-02,
01843       -1.81242E-02, 0.00000E+00,-1.00529E-01,-5.10574E-03, 0.00000E+00,
01844        2.10228E-03, 0.00000E+00, 0.00000E+00,-1.73255E+02, 5.07833E-01,
01845       -2.41408E-01, 8.75414E-03, 2.77527E-03,-8.90353E-05,-5.25148E+00,
01846       -5.83899E-03,-2.09122E-02,-9.63530E-03, 9.77164E-03, 4.07051E-03,
01847        2.53555E-04,-5.52875E+00,-3.55993E-01,-2.49231E-03, 0.00000E+00,
01848        0.00000E+00, 2.86026E+01, 0.00000E+00, 3.42722E-04, 0.00000E+00,
01849        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01850        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
01851    }, /* O DENSITY */ {
01852        1.02315E+00,-1.59710E-01,-1.06630E-01,-1.77074E-02,-4.42726E-03,
01853        3.44803E-02, 4.45613E-02,-3.33751E-02,-5.73598E-02, 3.50360E-01,
01854        6.33053E-02, 2.16221E-02, 5.42577E-02,-5.74193E+00, 0.00000E+00,
01855        1.90891E-01,-1.39194E-02, 1.01102E+02, 8.16363E-02, 1.33717E-04,
01856        6.54403E-06, 3.10295E-03, 0.00000E+00, 0.00000E+00, 5.38205E-02,
01857        1.23910E-01,-1.39831E-02, 0.00000E+00, 0.00000E+00,-3.95915E-06,
01858        0.00000E+00,-7.14651E-01,-5.01027E-03, 0.00000E+00,-3.24756E-03,
01859        0.00000E+00, 0.00000E+00, 4.42173E-02,-1.31598E+01,-3.15626E-03,
01860        1.24574E-03,-1.47626E-03,-1.55461E-03, 6.40682E-02, 1.34898E-01,
01861       -2.42415E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 6.13666E-04,
01862       -5.40373E-03, 2.61635E-05,-3.33012E-03, 0.00000E+00,-3.08101E-03,
01863       -2.42679E-03,-3.36086E-03, 0.00000E+00,-1.18979E+03,-5.04738E-02,
01864       -2.61547E-03,-1.03132E-03, 1.91583E-04,-8.38132E+01,-1.40517E-02,
01865       -1.14167E-02,-4.08012E-03, 1.73522E-04,-1.39644E-02,-6.64128E-02,
01866       -6.85152E-02,-1.34414E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01867        6.07916E+02,-4.12220E-03,-2.20996E-03, 0.00000E+00, 1.70277E+03,
01868       -4.63015E-03, 0.00000E+00, 0.00000E+00,-2.25360E-03,-2.96204E-03,
01869        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01870        3.92786E-02, 1.31186E-02,-1.78086E-03, 0.00000E+00, 0.00000E+00,
01871       -3.90083E-01,-2.84741E-02,-7.78400E+01,-1.02601E-03, 0.00000E+00,
01872       -7.26485E-04,-5.42181E-03,-5.59305E-03, 1.22825E-02, 1.23868E-02,
01873        6.68835E-03,-1.03303E-02,-9.51903E-03, 2.70021E-04,-2.57084E-02,
01874       -1.32430E-02, 0.00000E+00,-3.81000E-02,-3.16810E-03, 0.00000E+00,
01875        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01876        0.00000E+00,-9.05762E-04,-2.14590E-03,-1.17824E-03, 3.66732E+00,
01877       -3.79729E-04,-6.13966E-03,-5.09082E-03,-1.96332E-03,-3.08280E-03,
01878       -9.75222E-04, 4.03315E+00,-2.52710E-01, 0.00000E+00, 0.00000E+00,
01879        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01880        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01881        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
01882    }, /* N2 DENSITY */ {
01883        1.16112E+00, 0.00000E+00, 0.00000E+00, 3.33725E-02, 0.00000E+00,
01884        3.48637E-02,-5.44368E-03, 0.00000E+00,-6.73940E-02, 1.74754E-01,
01885        0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00,
01886        1.26733E-01, 0.00000E+00, 1.03154E+02, 5.52075E-02, 0.00000E+00,
01887        0.00000E+00, 8.13525E-04, 0.00000E+00, 0.00000E+00, 8.66784E-02,
01888        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01889        0.00000E+00,-2.50482E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01890        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.48894E-03,
01891        6.16053E-04,-5.79716E-04, 2.95482E-03, 8.47001E-02, 1.70147E-01,
01892        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01893        0.00000E+00, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01894        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01895        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01896        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01897        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01898        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01899        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01900        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01901        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01902        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01903        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01904        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01905        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01906        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01907        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01908        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01909        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01910        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01911        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01912        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
01913    }, /* TLB */ {
01914        9.44846E-01, 0.00000E+00, 0.00000E+00,-3.08617E-02, 0.00000E+00,
01915       -2.44019E-02, 6.48607E-03, 0.00000E+00, 3.08181E-02, 4.59392E-02,
01916        0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00,
01917        2.13260E-02, 0.00000E+00,-3.56958E+02, 0.00000E+00, 1.82278E-04,
01918        0.00000E+00, 3.07472E-04, 0.00000E+00, 0.00000E+00, 8.66784E-02,
01919        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01920        0.00000E+00, 0.00000E+00, 3.83054E-03, 0.00000E+00, 0.00000E+00,
01921       -1.93065E-03,-1.45090E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01922        0.00000E+00,-1.23493E-03, 1.36736E-03, 8.47001E-02, 1.70147E-01,
01923        3.71469E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01924        5.10250E-03, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01925        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01926        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01927        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01928        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01929        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01930        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01931        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01932        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01933        0.00000E+00, 3.68756E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01934        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01935        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01936        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01937        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01938        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01939        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01940        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01941        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01942        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01943        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
01944    }, /* O2 DENSITY */ {
01945        1.35580E+00, 1.44816E-01, 0.00000E+00, 6.07767E-02, 0.00000E+00,
01946        2.94777E-02, 7.46900E-02, 0.00000E+00,-9.23822E-02, 8.57342E-02,
01947        0.00000E+00, 0.00000E+00, 0.00000E+00, 2.38636E+01, 0.00000E+00,
01948        7.71653E-02, 0.00000E+00, 8.18751E+01, 1.87736E-02, 0.00000E+00,
01949        0.00000E+00, 1.49667E-02, 0.00000E+00, 0.00000E+00, 8.66784E-02,
01950        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01951        0.00000E+00,-3.67874E+02, 5.48158E-03, 0.00000E+00, 0.00000E+00,
01952        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01953        0.00000E+00, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
01954        1.22631E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01955        8.17187E-03, 3.71617E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01956        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01957        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.10826E-03,
01958       -3.13640E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01959        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01960        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01961        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01962        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01963       -7.35742E-02,-5.00266E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01964        0.00000E+00, 1.94965E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01965        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01966        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01967        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01968        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01969        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01970        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01971        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01972        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01973        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01974        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
01975    }, /* AR DENSITY */ {
01976        1.04761E+00, 2.00165E-01, 2.37697E-01, 3.68552E-02, 0.00000E+00,
01977        3.57202E-02,-2.14075E-01, 0.00000E+00,-1.08018E-01,-3.73981E-01,
01978        0.00000E+00, 3.10022E-02,-1.16305E-03,-2.07596E+01, 0.00000E+00,
01979        8.64502E-02, 0.00000E+00, 9.74908E+01, 5.16707E-02, 0.00000E+00,
01980        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 8.66784E-02,
01981        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01982        0.00000E+00, 3.46193E+02, 1.34297E-02, 0.00000E+00, 0.00000E+00,
01983        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.48509E-03,
01984       -1.54689E-04, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
01985        1.47753E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01986        1.89320E-02, 3.68181E-05, 1.32570E-02, 0.00000E+00, 0.00000E+00,
01987        3.59719E-03, 7.44328E-03,-1.00023E-03,-6.50528E+03, 0.00000E+00,
01988        1.03485E-02,-1.00983E-03,-4.06916E-03,-6.60864E+01,-1.71533E-02,
01989        1.10605E-02, 1.20300E-02,-5.20034E-03, 0.00000E+00, 0.00000E+00,
01990        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01991       -2.62769E+03, 7.13755E-03, 4.17999E-03, 0.00000E+00, 1.25910E+04,
01992        0.00000E+00, 0.00000E+00, 0.00000E+00,-2.23595E-03, 4.60217E-03,
01993        5.71794E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01994       -3.18353E-02,-2.35526E-02,-1.36189E-02, 0.00000E+00, 0.00000E+00,
01995        0.00000E+00, 2.03522E-02,-6.67837E+01,-1.09724E-03, 0.00000E+00,
01996       -1.38821E-02, 1.60468E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
01997        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.51574E-02,
01998       -5.44470E-04, 0.00000E+00, 7.28224E-02, 6.59413E-02, 0.00000E+00,
01999       -5.15692E-03, 0.00000E+00, 0.00000E+00,-3.70367E+03, 0.00000E+00,
02000        0.00000E+00, 1.36131E-02, 5.38153E-03, 0.00000E+00, 4.76285E+00,
02001       -1.75677E-02, 2.26301E-02, 0.00000E+00, 1.76631E-02, 4.77162E-03,
02002        0.00000E+00, 5.39354E+00, 0.00000E+00,-7.51710E-03, 0.00000E+00,
02003        0.00000E+00,-8.82736E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02004        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02005        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
02006    }, /* H DENSITY */ {
02007        1.26376E+00,-2.14304E-01,-1.49984E-01, 2.30404E-01, 2.98237E-02,
02008        2.68673E-02, 2.96228E-01, 2.21900E-02,-2.07655E-02, 4.52506E-01,
02009        1.20105E-01, 3.24420E-02, 4.24816E-02,-9.14313E+00, 0.00000E+00,
02010        2.47178E-02,-2.88229E-02, 8.12805E+01, 5.10380E-02,-5.80611E-03,
02011        2.51236E-05,-1.24083E-02, 0.00000E+00, 0.00000E+00, 8.66784E-02,
02012        1.58727E-01,-3.48190E-02, 0.00000E+00, 0.00000E+00, 2.89885E-05,
02013        0.00000E+00, 1.53595E+02,-1.68604E-02, 0.00000E+00, 1.01015E-02,
02014        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.84552E-04,
02015       -1.22181E-03, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
02016       -1.04927E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,-5.91313E-03,
02017       -2.30501E-02, 3.14758E-05, 0.00000E+00, 0.00000E+00, 1.26956E-02,
02018        8.35489E-03, 3.10513E-04, 0.00000E+00, 3.42119E+03,-2.45017E-03,
02019       -4.27154E-04, 5.45152E-04, 1.89896E-03, 2.89121E+01,-6.49973E-03,
02020       -1.93855E-02,-1.48492E-02, 0.00000E+00,-5.10576E-02, 7.87306E-02,
02021        9.51981E-02,-1.49422E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02022        2.65503E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02023        0.00000E+00, 0.00000E+00, 0.00000E+00, 6.37110E-03, 3.24789E-04,
02024        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02025        6.14274E-02, 1.00376E-02,-8.41083E-04, 0.00000E+00, 0.00000E+00,
02026        0.00000E+00,-1.27099E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02027       -3.94077E-03,-1.28601E-02,-7.97616E-03, 0.00000E+00, 0.00000E+00,
02028        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02029        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02030        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02031        0.00000E+00,-6.71465E-03,-1.69799E-03, 1.93772E-03, 3.81140E+00,
02032       -7.79290E-03,-1.82589E-02,-1.25860E-02,-1.04311E-02,-3.02465E-03,
02033        2.43063E-03, 3.63237E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02034        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02035        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02036        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
02037    }, /* N DENSITY */ {
02038        7.09557E+01,-3.26740E-01, 0.00000E+00,-5.16829E-01,-1.71664E-03,
02039        9.09310E-02,-6.71500E-01,-1.47771E-01,-9.27471E-02,-2.30862E-01,
02040       -1.56410E-01, 1.34455E-02,-1.19717E-01, 2.52151E+00, 0.00000E+00,
02041       -2.41582E-01, 5.92939E-02, 4.39756E+00, 9.15280E-02, 4.41292E-03,
02042        0.00000E+00, 8.66807E-03, 0.00000E+00, 0.00000E+00, 8.66784E-02,
02043        1.58727E-01, 9.74701E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02044        0.00000E+00, 6.70217E+01,-1.31660E-03, 0.00000E+00,-1.65317E-02,
02045        0.00000E+00, 0.00000E+00, 8.50247E-02, 2.77428E+01, 4.98658E-03,
02046        6.15115E-03, 9.50156E-03,-2.12723E-02, 8.47001E-02, 1.70147E-01,
02047       -2.38645E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.37380E-03,
02048       -8.41918E-03, 2.80145E-05, 7.12383E-03, 0.00000E+00,-1.66209E-02,
02049        1.03533E-04,-1.68898E-02, 0.00000E+00, 3.64526E+03, 0.00000E+00,
02050        6.54077E-03, 3.69130E-04, 9.94419E-04, 8.42803E+01,-1.16124E-02,
02051       -7.74414E-03,-1.68844E-03, 1.42809E-03,-1.92955E-03, 1.17225E-01,
02052       -2.41512E-02, 1.50521E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02053        1.60261E+03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02054        0.00000E+00, 0.00000E+00, 0.00000E+00,-3.54403E-04,-1.87270E-02,
02055        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02056        2.76439E-02, 6.43207E-03,-3.54300E-02, 0.00000E+00, 0.00000E+00,
02057        0.00000E+00,-2.80221E-02, 8.11228E+01,-6.75255E-04, 0.00000E+00,
02058       -1.05162E-02,-3.48292E-03,-6.97321E-03, 0.00000E+00, 0.00000E+00,
02059        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02060        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02061        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02062        0.00000E+00,-1.45546E-03,-1.31970E-02,-3.57751E-03,-1.09021E+00,
02063       -1.50181E-02,-7.12841E-03,-6.64590E-03,-3.52610E-03,-1.87773E-02,
02064       -2.22432E-03,-3.93895E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02065        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02066        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02067        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 
02068    }, /* HOT O DENSITY */ {
02069        6.04050E-02, 1.57034E+00, 2.99387E-02, 0.00000E+00, 0.00000E+00,
02070        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.51018E+00,
02071        0.00000E+00, 0.00000E+00, 0.00000E+00,-8.61650E+00, 1.26454E-02,
02072        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02073        0.00000E+00, 5.50878E-03, 0.00000E+00, 0.00000E+00, 8.66784E-02,
02074        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02075        0.00000E+00, 0.00000E+00, 6.23881E-02, 0.00000E+00, 0.00000E+00,
02076        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02077        0.00000E+00, 0.00000E+00, 0.00000E+00, 8.47001E-02, 1.70147E-01,
02078       -9.45934E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02079        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02080        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02081        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 
02082        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02083        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02084        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02085        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02086        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02087        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02088        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02089        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02090        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02091        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02092        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02093        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02094        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02095        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02096        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02097        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02098        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
02099    }};
02100    /* S PARAM  */
02101    double Msise00Drag::ps[150] = {
02102        9.56827E-01, 6.20637E-02, 3.18433E-02, 0.00000E+00, 0.00000E+00,
02103        3.94900E-02, 0.00000E+00, 0.00000E+00,-9.24882E-03,-7.94023E-03,
02104        0.00000E+00, 0.00000E+00, 0.00000E+00, 1.74712E+02, 0.00000E+00,
02105        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02106        0.00000E+00, 2.74677E-03, 0.00000E+00, 1.54951E-02, 8.66784E-02,
02107        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02108        0.00000E+00, 0.00000E+00, 0.00000E+00,-6.99007E-04, 0.00000E+00,
02109        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02110        0.00000E+00, 1.24362E-02,-5.28756E-03, 8.47001E-02, 1.70147E-01,
02111        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02112        0.00000E+00, 2.47425E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02113        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02114        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02115        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02116        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02117        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02118        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02119        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02120        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02121        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02122        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02123        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02124        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02125        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02126        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02127        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02128        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02129        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02130        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02131        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
02132    };
02133 
02134    /* TURBO */
02135    double Msise00Drag::pdl[2][25] = {
02136       { 1.09930E+00, 3.90631E+00, 3.07165E+00, 9.86161E-01, 1.63536E+01,
02137        4.63830E+00, 1.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02138        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02139        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02140        0.00000E+00, 0.00000E+00, 1.28840E+00, 3.10302E-02, 1.18339E-01 }, 
02141       { 1.00000E+00, 7.00000E-01, 1.15020E+00, 3.44689E+00, 1.28840E+00,
02142        1.00000E+00, 1.08738E+00, 1.22947E+00, 1.10016E+00, 7.34129E-01,
02143        1.15241E+00, 2.22784E+00, 7.95046E-01, 4.01612E+00, 4.47749E+00,
02144        1.23435E+02,-7.60535E-02, 1.68986E-06, 7.44294E-01, 1.03604E+00,
02145        1.72783E+02, 1.15020E+00, 3.44689E+00,-7.46230E-01, 9.49154E-01 }
02146    };
02147    /* LOWER BOUNDARY */
02148    double Msise00Drag::ptm[10]/* ptm[50]*/ = {
02149        1.04130E+03, 3.86000E+02, 1.95000E+02, 1.66728E+01, 2.13000E+02,
02150        1.20000E+02, 2.40000E+02, 1.87000E+02,-2.00000E+00, 0.00000E+00
02151    };
02152    double Msise00Drag::pdm[8][10] = {
02153    {    2.45600E+07, 6.71072E-06, 1.00000E+02, 0.00000E+00, 1.10000E+02,
02154        1.00000E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 },\
02155    {    8.59400E+10, 1.00000E+00, 1.05000E+02,-8.00000E+00, 1.10000E+02,
02156        1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00 },\
02157    {    2.81000E+11, 0.00000E+00, 1.05000E+02, 2.80000E+01, 2.89500E+01,
02158        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 },
02159    {    3.30000E+10, 2.68270E-01, 1.05000E+02, 1.00000E+00, 1.10000E+02,
02160        1.00000E+01, 1.10000E+02,-1.00000E+01, 0.00000E+00, 0.00000E+00 },
02161    {    1.33000E+09, 1.19615E-02, 1.05000E+02, 0.00000E+00, 1.10000E+02,
02162        1.00000E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00 },
02163    {    1.76100E+05, 1.00000E+00, 9.50000E+01,-8.00000E+00, 1.10000E+02,
02164        1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00, },
02165    {    1.00000E+07, 1.00000E+00, 1.05000E+02,-8.00000E+00, 1.10000E+02,
02166        1.00000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 0.00000E+00 },
02167    {    1.00000E+06, 1.00000E+00, 1.05000E+02,-8.00000E+00, 5.50000E+02,
02168        7.60000E+01, 9.00000E+01, 2.00000E+00, 0.00000E+00, 4.00000E+03 }};
02169 
02170 
02171    double Msise00Drag::ptl[4][100] = {
02172    /* TN1(2) */ {
02173        1.00858E+00, 4.56011E-02,-2.22972E-02,-5.44388E-02, 5.23136E-04,
02174       -1.88849E-02, 5.23707E-02,-9.43646E-03, 6.31707E-03,-7.80460E-02,
02175       -4.88430E-02, 0.00000E+00, 0.00000E+00,-7.60250E+00, 0.00000E+00,
02176       -1.44635E-02,-1.76843E-02,-1.21517E+02, 2.85647E-02, 0.00000E+00,
02177        0.00000E+00, 6.31792E-04, 0.00000E+00, 5.77197E-03, 8.66784E-02,
02178        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02179        0.00000E+00,-8.90272E+03, 3.30611E-03, 3.02172E-03, 0.00000E+00,
02180       -2.13673E-03,-3.20910E-04, 0.00000E+00, 0.00000E+00, 2.76034E-03,
02181        2.82487E-03,-2.97592E-04,-4.21534E-03, 8.47001E-02, 1.70147E-01,
02182        8.96456E-03, 0.00000E+00,-1.08596E-02, 0.00000E+00, 0.00000E+00,
02183        5.57917E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02184        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02185        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02186        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02187        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02188        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02189        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02190        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02191        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02192        0.00000E+00, 9.65405E-03, 0.00000E+00, 0.00000E+00, 2.00000E+00
02193    }, /*  TN1(3) */ {
02194        9.39664E-01, 8.56514E-02,-6.79989E-03, 2.65929E-02,-4.74283E-03,
02195        1.21855E-02,-2.14905E-02, 6.49651E-03,-2.05477E-02,-4.24952E-02,
02196        0.00000E+00, 0.00000E+00, 0.00000E+00, 1.19148E+01, 0.00000E+00,
02197        1.18777E-02,-7.28230E-02,-8.15965E+01, 1.73887E-02, 0.00000E+00,
02198        0.00000E+00, 0.00000E+00,-1.44691E-02, 2.80259E-04, 8.66784E-02,
02199        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02200        0.00000E+00, 2.16584E+02, 3.18713E-03, 7.37479E-03, 0.00000E+00,
02201       -2.55018E-03,-3.92806E-03, 0.00000E+00, 0.00000E+00,-2.89757E-03,
02202       -1.33549E-03, 1.02661E-03, 3.53775E-04, 8.47001E-02, 1.70147E-01,
02203       -9.17497E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02204        3.56082E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02205        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02206        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02207        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02208        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02209        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02210        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02211        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02212        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02213        0.00000E+00,-1.00902E-02, 0.00000E+00, 0.00000E+00, 2.00000E+00
02214    }, /* TN1(4) */ {
02215        9.85982E-01,-4.55435E-02, 1.21106E-02, 2.04127E-02,-2.40836E-03,
02216        1.11383E-02,-4.51926E-02, 1.35074E-02,-6.54139E-03, 1.15275E-01,
02217        1.28247E-01, 0.00000E+00, 0.00000E+00,-5.30705E+00, 0.00000E+00,
02218       -3.79332E-02,-6.24741E-02, 7.71062E-01, 2.96315E-02, 0.00000E+00,
02219        0.00000E+00, 0.00000E+00, 6.81051E-03,-4.34767E-03, 8.66784E-02,
02220        1.58727E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02221        0.00000E+00, 1.07003E+01,-2.76907E-03, 4.32474E-04, 0.00000E+00,
02222        1.31497E-03,-6.47517E-04, 0.00000E+00,-2.20621E+01,-1.10804E-03,
02223       -8.09338E-04, 4.18184E-04, 4.29650E-03, 8.47001E-02, 1.70147E-01,
02224        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02225       -4.04337E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02226        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02227        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-9.52550E-04,
02228        8.56253E-04, 4.33114E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02229        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.21223E-03,
02230        2.38694E-04, 9.15245E-04, 1.28385E-03, 8.67668E-04,-5.61425E-06,
02231        1.04445E+00, 3.41112E+01, 0.00000E+00,-8.40704E-01,-2.39639E+02,
02232        7.06668E-01,-2.05873E+01,-3.63696E-01, 2.39245E+01, 0.00000E+00,
02233       -1.06657E-03,-7.67292E-04, 1.54534E-04, 0.00000E+00, 0.00000E+00,
02234        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02235    }, /* TN1(5) TN2(1) */ {
02236        1.00320E+00, 3.83501E-02,-2.38983E-03, 2.83950E-03, 4.20956E-03,
02237        5.86619E-04, 2.19054E-02,-1.00946E-02,-3.50259E-03, 4.17392E-02,
02238       -8.44404E-03, 0.00000E+00, 0.00000E+00, 4.96949E+00, 0.00000E+00,
02239       -7.06478E-03,-1.46494E-02, 3.13258E+01,-1.86493E-03, 0.00000E+00,
02240       -1.67499E-02, 0.00000E+00, 0.00000E+00, 5.12686E-04, 8.66784E-02,
02241        1.58727E-01,-4.64167E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02242        4.37353E-03,-1.99069E+02, 0.00000E+00,-5.34884E-03, 0.00000E+00,
02243        1.62458E-03, 2.93016E-03, 2.67926E-03, 5.90449E+02, 0.00000E+00,
02244        0.00000E+00,-1.17266E-03,-3.58890E-04, 8.47001E-02, 1.70147E-01,
02245        0.00000E+00, 0.00000E+00, 1.38673E-02, 0.00000E+00, 0.00000E+00,
02246        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02247        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02248        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.60571E-03,
02249        6.28078E-04, 5.05469E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02250        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.57829E-03,
02251       -4.00855E-04, 5.04077E-05,-1.39001E-03,-2.33406E-03,-4.81197E-04,
02252        1.46758E+00, 6.20332E+00, 0.00000E+00, 3.66476E-01,-6.19760E+01,
02253        3.09198E-01,-1.98999E+01, 0.00000E+00,-3.29933E+02, 0.00000E+00,
02254       -1.10080E-03,-9.39310E-05, 1.39638E-04, 0.00000E+00, 0.00000E+00,
02255        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02256    } };
02257 
02258    double Msise00Drag::pma[10][100] = {
02259    /* TN2(2) */ {
02260        9.81637E-01,-1.41317E-03, 3.87323E-02, 0.00000E+00, 0.00000E+00,
02261        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.58707E-02,
02262       -8.63658E-03, 0.00000E+00, 0.00000E+00,-2.02226E+00, 0.00000E+00,
02263       -8.69424E-03,-1.91397E-02, 8.76779E+01, 4.52188E-03, 0.00000E+00,
02264        2.23760E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02265        0.00000E+00,-7.07572E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02266       -4.11210E-03, 3.50060E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02267        0.00000E+00, 0.00000E+00,-8.36657E-03, 1.61347E+01, 0.00000E+00,
02268        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02269        0.00000E+00, 0.00000E+00,-1.45130E-02, 0.00000E+00, 0.00000E+00,
02270        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02271        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02272        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.24152E-03,
02273        6.43365E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02274        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.33255E-03,
02275        2.42657E-03, 1.60666E-03,-1.85728E-03,-1.46874E-03,-4.79163E-06,
02276        1.22464E+00, 3.53510E+01, 0.00000E+00, 4.49223E-01,-4.77466E+01,
02277        4.70681E-01, 8.41861E+00,-2.88198E-01, 1.67854E+02, 0.00000E+00,
02278        7.11493E-04, 6.05601E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02279        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02280    }, /* TN2(3) */ {
02281        1.00422E+00,-7.11212E-03, 5.24480E-03, 0.00000E+00, 0.00000E+00,
02282        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-5.28914E-02,
02283       -2.41301E-02, 0.00000E+00, 0.00000E+00,-2.12219E+01,-1.03830E-02,
02284       -3.28077E-03, 1.65727E-02, 1.68564E+00,-6.68154E-03, 0.00000E+00,
02285        1.45155E-02, 0.00000E+00, 8.42365E-03, 0.00000E+00, 0.00000E+00,
02286        0.00000E+00,-4.34645E-03, 0.00000E+00, 0.00000E+00, 2.16780E-02,
02287        0.00000E+00,-1.38459E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02288        0.00000E+00, 0.00000E+00, 7.04573E-03,-4.73204E+01, 0.00000E+00,
02289        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02290        0.00000E+00, 0.00000E+00, 1.08767E-02, 0.00000E+00, 0.00000E+00,
02291        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02292        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-8.08279E-03,
02293        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.21769E-04,
02294       -2.27387E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02295        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.26769E-03,
02296        3.16901E-03, 4.60316E-04,-1.01431E-04, 1.02131E-03, 9.96601E-04,
02297        1.25707E+00, 2.50114E+01, 0.00000E+00, 4.24472E-01,-2.77655E+01,
02298        3.44625E-01, 2.75412E+01, 0.00000E+00, 7.94251E+02, 0.00000E+00,
02299        2.45835E-03, 1.38871E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02300        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02301    }, /* TN2(4) TN3(1) */ {
02302        1.01890E+00,-2.46603E-02, 1.00078E-02, 0.00000E+00, 0.00000E+00,
02303        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-6.70977E-02,
02304       -4.02286E-02, 0.00000E+00, 0.00000E+00,-2.29466E+01,-7.47019E-03,
02305        2.26580E-03, 2.63931E-02, 3.72625E+01,-6.39041E-03, 0.00000E+00,
02306        9.58383E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02307        0.00000E+00,-1.85291E-03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02308        0.00000E+00, 1.39717E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02309        0.00000E+00, 0.00000E+00, 9.19771E-03,-3.69121E+02, 0.00000E+00,
02310        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02311        0.00000E+00, 0.00000E+00,-1.57067E-02, 0.00000E+00, 0.00000E+00,
02312        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02313        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.07265E-03,
02314        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.92953E-03,
02315       -2.77739E-03,-4.40092E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02316        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.47280E-03,
02317        2.95035E-04,-1.81246E-03, 2.81945E-03, 4.27296E-03, 9.78863E-04,
02318        1.40545E+00,-6.19173E+00, 0.00000E+00, 0.00000E+00,-7.93632E+01,
02319        4.44643E-01,-4.03085E+02, 0.00000E+00, 1.15603E+01, 0.00000E+00,
02320        2.25068E-03, 8.48557E-04,-2.98493E-04, 0.00000E+00, 0.00000E+00,
02321        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02322    }, /* TN3(2) */ {
02323        9.75801E-01, 3.80680E-02,-3.05198E-02, 0.00000E+00, 0.00000E+00,
02324        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.85575E-02,
02325        5.04057E-02, 0.00000E+00, 0.00000E+00,-1.76046E+02, 1.44594E-02,
02326       -1.48297E-03,-3.68560E-03, 3.02185E+01,-3.23338E-03, 0.00000E+00,
02327        1.53569E-02, 0.00000E+00,-1.15558E-02, 0.00000E+00, 0.00000E+00,
02328        0.00000E+00, 4.89620E-03, 0.00000E+00, 0.00000E+00,-1.00616E-02,
02329       -8.21324E-03,-1.57757E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02330        0.00000E+00, 0.00000E+00, 6.63564E-03, 4.58410E+01, 0.00000E+00,
02331        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02332        0.00000E+00, 0.00000E+00,-2.51280E-02, 0.00000E+00, 0.00000E+00,
02333        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02334        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 9.91215E-03,
02335        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-8.73148E-04,
02336       -1.29648E-03,-7.32026E-05, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02337        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-4.68110E-03,
02338       -4.66003E-03,-1.31567E-03,-7.39390E-04, 6.32499E-04,-4.65588E-04,
02339       -1.29785E+00,-1.57139E+02, 0.00000E+00, 2.58350E-01,-3.69453E+01,
02340        4.10672E-01, 9.78196E+00,-1.52064E-01,-3.85084E+03, 0.00000E+00,
02341       -8.52706E-04,-1.40945E-03,-7.26786E-04, 0.00000E+00, 0.00000E+00,
02342        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02343    }, /* TN3(3) */ {
02344        9.60722E-01, 7.03757E-02,-3.00266E-02, 0.00000E+00, 0.00000E+00,
02345        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.22671E-02,
02346        4.10423E-02, 0.00000E+00, 0.00000E+00,-1.63070E+02, 1.06073E-02,
02347        5.40747E-04, 7.79481E-03, 1.44908E+02, 1.51484E-04, 0.00000E+00,
02348        1.97547E-02, 0.00000E+00,-1.41844E-02, 0.00000E+00, 0.00000E+00,
02349        0.00000E+00, 5.77884E-03, 0.00000E+00, 0.00000E+00, 9.74319E-03,
02350        0.00000E+00,-2.88015E+03, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02351        0.00000E+00, 0.00000E+00,-4.44902E-03,-2.92760E+01, 0.00000E+00,
02352        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02353        0.00000E+00, 0.00000E+00, 2.34419E-02, 0.00000E+00, 0.00000E+00,
02354        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02355        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.36685E-03,
02356        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-4.65325E-04,
02357       -5.50628E-04, 3.31465E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02358        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.06179E-03,
02359       -3.08575E-03,-7.93589E-04,-1.08629E-04, 5.95511E-04,-9.05050E-04,
02360        1.18997E+00, 4.15924E+01, 0.00000E+00,-4.72064E-01,-9.47150E+02,
02361        3.98723E-01, 1.98304E+01, 0.00000E+00, 3.73219E+03, 0.00000E+00,
02362       -1.50040E-03,-1.14933E-03,-1.56769E-04, 0.00000E+00, 0.00000E+00,
02363        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02364    }, /* TN3(4) */ {
02365        1.03123E+00,-7.05124E-02, 8.71615E-03, 0.00000E+00, 0.00000E+00,
02366        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-3.82621E-02,
02367       -9.80975E-03, 0.00000E+00, 0.00000E+00, 2.89286E+01, 9.57341E-03,
02368        0.00000E+00, 0.00000E+00, 8.66153E+01, 7.91938E-04, 0.00000E+00,
02369        0.00000E+00, 0.00000E+00, 4.68917E-03, 0.00000E+00, 0.00000E+00,
02370        0.00000E+00, 7.86638E-03, 0.00000E+00, 0.00000E+00, 9.90827E-03,
02371        0.00000E+00, 6.55573E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02372        0.00000E+00, 0.00000E+00, 0.00000E+00,-4.00200E+01, 0.00000E+00,
02373        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02374        0.00000E+00, 0.00000E+00, 7.07457E-03, 0.00000E+00, 0.00000E+00,
02375        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02376        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.72268E-03,
02377        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.04970E-04,
02378        1.21560E-03,-8.05579E-06, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02379        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.49941E-03,
02380       -4.57256E-04,-1.59311E-04, 2.96481E-04,-1.77318E-03,-6.37918E-04,
02381        1.02395E+00, 1.28172E+01, 0.00000E+00, 1.49903E-01,-2.63818E+01,
02382        0.00000E+00, 4.70628E+01,-2.22139E-01, 4.82292E-02, 0.00000E+00,
02383       -8.67075E-04,-5.86479E-04, 5.32462E-04, 0.00000E+00, 0.00000E+00,
02384        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02385    }, /* TN3(5) SURFACE TEMP TSL */ {
02386        1.00828E+00,-9.10404E-02,-2.26549E-02, 0.00000E+00, 0.00000E+00,
02387        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-2.32420E-02,
02388       -9.08925E-03, 0.00000E+00, 0.00000E+00, 3.36105E+01, 0.00000E+00,
02389        0.00000E+00, 0.00000E+00,-1.24957E+01,-5.87939E-03, 0.00000E+00,
02390        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02391        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02392        0.00000E+00, 2.79765E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02393        0.00000E+00, 0.00000E+00, 0.00000E+00, 2.01237E+03, 0.00000E+00,
02394        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02395        0.00000E+00, 0.00000E+00,-1.75553E-02, 0.00000E+00, 0.00000E+00,
02396        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02397        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02398        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.29699E-03,
02399        1.26659E-03, 2.68402E-04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02400        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 1.17894E-03,
02401        1.48746E-03, 1.06478E-04, 1.34743E-04,-2.20939E-03,-6.23523E-04,
02402        6.36539E-01, 1.13621E+01, 0.00000E+00,-3.93777E-01, 2.38687E+03,
02403        0.00000E+00, 6.61865E+02,-1.21434E-01, 9.27608E+00, 0.00000E+00,
02404        1.68478E-04, 1.24892E-03, 1.71345E-03, 0.00000E+00, 0.00000E+00,
02405        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02406    }, /* TGN3(2) SURFACE GRAD TSLG */ {
02407        1.57293E+00,-6.78400E-01, 6.47500E-01, 0.00000E+00, 0.00000E+00,
02408        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-7.62974E-02,
02409       -3.60423E-01, 0.00000E+00, 0.00000E+00, 1.28358E+02, 0.00000E+00,
02410        0.00000E+00, 0.00000E+00, 4.68038E+01, 0.00000E+00, 0.00000E+00,
02411        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02412        0.00000E+00,-1.67898E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02413        0.00000E+00, 2.90994E+04, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02414        0.00000E+00, 0.00000E+00, 0.00000E+00, 3.15706E+01, 0.00000E+00,
02415        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02416        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02417        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02418        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02419        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02420        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02421        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02422        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02423        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02424        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02425        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02426        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02427    }, /* TGN2(1) TGN1(2) */ {
02428        8.60028E-01, 3.77052E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02429        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.17570E+00,
02430        0.00000E+00, 0.00000E+00, 0.00000E+00, 7.77757E-03, 0.00000E+00,
02431        0.00000E+00, 0.00000E+00, 1.01024E+02, 0.00000E+00, 0.00000E+00,
02432        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02433        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02434        0.00000E+00, 6.54251E+02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02435        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02436        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02437        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02438        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02439        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02440        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02441        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02442        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,-1.56959E-02,
02443        1.91001E-02, 3.15971E-02, 1.00982E-02,-6.71565E-03, 2.57693E-03,
02444        1.38692E+00, 2.82132E-01, 0.00000E+00, 0.00000E+00, 3.81511E+02,
02445        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02446        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02447        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02448    }, /* TGN3(1) TGN2(2) */ {
02449        1.06029E+00,-5.25231E-02, 3.73034E-01, 0.00000E+00, 0.00000E+00,
02450        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.31072E-02,
02451       -3.88409E-01, 0.00000E+00, 0.00000E+00,-1.65295E+02,-2.13801E-01,
02452       -4.38916E-02,-3.22716E-01,-8.82393E+01, 1.18458E-01, 0.00000E+00,
02453       -4.35863E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02454        0.00000E+00,-1.19782E-01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02455        0.00000E+00, 2.62229E+01, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02456        0.00000E+00, 0.00000E+00, 0.00000E+00,-5.37443E+01, 0.00000E+00,
02457        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02458        0.00000E+00, 0.00000E+00,-4.55788E-01, 0.00000E+00, 0.00000E+00,
02459        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02460        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02461        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 3.84009E-02,
02462        3.96733E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02463        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 5.05494E-02,
02464        7.39617E-02, 1.92200E-02,-8.46151E-03,-1.34244E-02, 1.96338E-02,
02465        1.50421E+00, 1.88368E+01, 0.00000E+00, 0.00000E+00,-5.13114E+01,
02466        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02467        5.11923E-02, 3.61225E-02, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02468        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 2.00000E+00
02469    } };
02470     
02471    /* SEMIANNUAL MULT SAM */
02472    double Msise00Drag::sam[100] = {
02473        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02474        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02475        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02476        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02477        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02478        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02479        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02480        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02481        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02482        1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
02483        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02484        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02485        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02486        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02487        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02488        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02489        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02490        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02491        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00,
02492        0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00, 0.00000E+00
02493    };
02494     
02495     
02496    /* MIDDLE ATMOSPHERE AVERAGES */
02497    double Msise00Drag::pavgm[10] = {
02498        2.61000E+02, 2.64000E+02, 2.29000E+02, 2.17000E+02, 2.17000E+02,
02499        2.23000E+02, 2.86760E+02,-2.93940E+00, 2.50000E+00, 0.00000E+00 };
02500 
02501 
02503 
02504    /* PARMB */
02505    double Msise00Drag::gsurf = 0.0;
02506    double Msise00Drag::re = 0.0;
02507 
02508    /* GTS3C */
02509    double Msise00Drag::dd=0.0;
02510 
02511     /* DMIX */
02512    double Msise00Drag::dm04=0.0;
02513    double Msise00Drag::dm16=0.0;
02514    double Msise00Drag::dm28=0.0;
02515    double Msise00Drag::dm32=0.0;
02516    double Msise00Drag::dm40=0.0;
02517    double Msise00Drag::dm01=0.0;
02518    double Msise00Drag::dm14=0.0;
02519 
02520     /* MESO7 */
02521    double Msise00Drag::meso_tn1[5]={0.0};
02522    double Msise00Drag::meso_tn2[4]={0.0};
02523    double Msise00Drag::meso_tn3[5]={0.0};
02524    double Msise00Drag::meso_tgn1[2]={0.0};
02525    double Msise00Drag::meso_tgn2[2]={0.0};
02526    double Msise00Drag::meso_tgn3[2]={0.0};
02527 
02528     /* LPOLY */
02529    double Msise00Drag::dfa=0.0;
02530    double Msise00Drag::plg[4][9]={0.0};
02531    double Msise00Drag::ctloc=0.0;
02532    double Msise00Drag::stloc=0.0;
02533    double Msise00Drag::c2tloc=0.0;
02534    double Msise00Drag::s2tloc=0.0;
02535    double Msise00Drag::s3tloc=0.0;
02536    double Msise00Drag::c3tloc=0.0;
02537    double Msise00Drag::apdf=0.0;
02538    double Msise00Drag::apt[4]={0.0};
02539 
02540 
02541 }  // End of namespace 'gpstk'
02542 
02543 
02544 
02545 

Generated on Tue May 22 03:31:00 2012 for GPS ToolKit Software Library by  doxygen 1.3.9.1