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