NRLMSISE00.java
/* Copyright 2002-2017 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.forces.drag.atmosphere;
import java.util.Arrays;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinatesProvider;
/** This class implements the mathematical representation of the 2001
* Naval Research Laboratory Mass Spectrometer and Incoherent Scatter
* Radar Exosphere (NRLMSISE-00) of the MSIS® class model.
* <p>
* NRLMSISE-00 calculates the neutral atmosphere empirical model from the surface
* to lower exosphere (0 to 1000 km) and provides:
* <ul>
* <li>Exospheric Temperature above Input Position (K)</li>
* <li>Local Temperature at Input Position (K)</li>
* <li>Total Mass-Density at Input Position (kg/m³)</li>
* <li>Partial Densities at Input Position (1/m³) for:
* <ul>
* <li>He,</li>
* <li>H,</li>
* <li>N,</li>
* <li>O,</li>
* <li>Ar,</li>
* <li>N2,</li>
* <li>O2,</li>
* <li>anomalous oxygen.</li>
* </ul>
* </li>
* </ul>
* </p>
* <p>
* The model needs geographical and time information to compute general values,
* but also needs space weather data:
* <ul>
* <li>mean and daily solar flux,</li>
* <li>geomagnetic indices.</li>
* </ul>
* </p>
* <p>
* Switches can be used to turn on and off particular variations:<br>
* 0 is off, 1 is on, and 2 is main effects off but cross terms on.<br>
* The standard value is 1 for all the 23 available switches.<br>
* Function of each switch according to its number:
* <ul>
* <li>#1 - F10.7 effect on mean</li>
* <li>#2 - Independent of time</li>
* <li>#3 - Symmetrical annual</li>
* <li>#4 - Symmetrical semiannual</li>
* <li>#5 - Asymmetrical annual</li>
* <li>#6 - Asymmetrical semiannual</li>
* <li>#7 - Diurnal</li>
* <li>#8 - Semidiurnal</li>
* <li>#9 - Daily Ap [**]</li>
* <li>#10 - All UT, longitudinal effects</li>
* <li>#11 - Longitudinal</li>
* <li>#12 - UT and mixed UT, longitudinal</li>
* <li>#13 - Mixed AP, UT, longitudinal</li>
* <li>#14 - Terdiurnal</li>
* <li>#15 - Departures from diffusive equilibrium</li>
* <li>#16 - All exospheric temperature variations</li>
* <li>#17 - All variations from 120 km temperature (TLB)</li>
* <li>#18 - All lower thermosphere (TN1) temperature variations</li>
* <li>#19 - All 120 km gradient (S) variations</li>
* <li>#20 - All upper stratosphere (TN2) temperature variations</li>
* <li>#21 - All variations from 120 km values (ZLB)</li>
* <li>#22 - All lower mesosphere temperature (TN3) variations</li>
* <li>#23 - Turbopause scale height variations</li>
* </ul>
* [**] Switch #9 is a bit specific:
* <ul>
* <li>set to 1, the daily Ap only is used (first element of ap array),</li>
* <li>set to -1, the entire array of ap is used, including 3 hr ap indices.</li>
* </ul>
* </p>
* <p>
* The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and Doug Drob.<br>
* They also wrote a NRLMSISE-00 distribution package in FORTRAN available at:<br>
* ftp://hanna.ccmc.gsfc.nasa.gov/pub/modelweb/atmospheric/msis/nrlmsise00/<br>
* <br>
* Dominik Brodowski implemented a C version of the NRLMSISE-00 model available at:<br>
* http://www.brodo.de/space/nrlmsise/index.html
* </p>
* <p>
* Instances of this class are immutable.
* </p>
*
* @author Mike Picone & al (Naval Research Laboratory), 2001: FORTRAN routine
* @author Dominik Brodowski, 2004: C routine
* @author Pascal Parraud, 2016: Java translation
* @since 8.1
*/
public class NRLMSISE00 implements Atmosphere {
/** Serializable UID. */
private static final long serialVersionUID = -7923498628122574334L;
// Constants
/** Identifier for helium density. */
private static final int HELIUM = 0;
/** Identifier for atomic oxygen density. */
private static final int ATOMIC_OXYGEN = 1;
/** Identifier for molecular nitrogen density. */
private static final int MOLECULAR_NITROGEN = 2;
/** Identifier for molecular oxygen density. */
private static final int MOLECULAR_OXYGEN = 3;
/** Identifier for argon density. */
private static final int ARGON = 4;
/** Identifier for atomic nitrogen density. */
private static final int TOTAL_MASS = 5;
/** Identifier for hydrogen density. */
private static final int HYDROGEN = 6;
/** Identifier for atomic nitrogen density. */
private static final int ATOMIC_NITROGEN = 7;
/** Identifier for anomalous oxygen density. */
private static final int ANOMALOUS_OXYGEN = 8;
/** Identifier for exospheric temperature. */
private static final int EXOSPHERIC = 0;
/** Identifier for temperature at altitude. */
private static final int ALTITUDE = 1;
// CONVERSION CONSTANTS
/** Conversion from degree to radian. */
private static final double DEG_TO_RAD = 1.74533e-2;
/** Conversion from day to radian. */
private static final double DAY_TO_RAD = 1.72142e-2;
/** Conversion from hour to radian. */
private static final double HOUR_TO_RAD = 0.2618;
/** Conversion from second to radian. */
private static final double SEC_TO_RAD = 7.2722e-5;
// EARTH GEOPHYSICAL CONSTANTS
/** Reference latitude (°). */
private static final double LAT_REF = 45.;
/** Reference gravity on Earth surface at reference latitude (cm/s2). */
private static final double G_REF = 980.616;
// CHEMICAL CONSTANTS
/** Unified atomic mass unit (kg). */
private static final double AMU = 1.66e-27;
/** Gas constant (inverse of). */
private static final double R_GAS = 831.4;
/** Hydrogen atomic mass. */
private static final double H_MASS = 1.;
/** Helium atomic mass. */
private static final double HE_MASS = 4.;
/** Nitrogen atomic mass. */
private static final double N_MASS = 14.;
/** N2 molecular mass. */
private static final double N2_MASS = 2. * N_MASS;
/** Oxygen atomic mass. */
private static final double O_MASS = 16.;
/** O2 molecular mass. */
private static final double O2_MASS = 2. * O_MASS;
/** Argon atomic mass. */
private static final double AR_MASS = 40.;
// NRL MSISE 2000 SPECIFIC CONSTANTS
/** Reference average flux. */
private static final double FLUX_REF = 150.;
/** Array of altitudes #1. */
private static final double[] ZN1 = {123.435, 110.0, 100.0, 90.0, 72.5};
/** Array of altitudes #2. */
private static final double[] ZN2 = {72.5, 55.0, 45.0, 32.5};
/** Array of altitudes #3. */
private static final double[] ZN3 = {32.5, 20.0, 15.0, 10.0, 0.0};
/** Mix altitude (km). */
private static final double ZMIX = 62.5;
/** NRLMSISE-00 data: temperature pt[150]. */
private static final double[] PT = {
9.86573e-01, 1.62228e-02, 1.55270e-02, -1.04323e-01, -3.75801e-03,
-1.18538e-03, -1.24043e-01, 4.56820e-03, 8.76018e-03, -1.36235e-01,
-3.52427e-02, 8.84181e-03, -5.92127e-03, -8.61650e+00, 0.00000e+00,
1.28492e-02, 0.00000e+00, 1.30096e+02, 1.04567e-02, 1.65686e-03,
-5.53887e-06, 2.97810e-03, 0.00000e+00, 5.13122e-03, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.27026e-06,
0.00000e+00, 6.74494e+00, 4.93933e-03, 2.21656e-03, 2.50802e-03,
0.00000e+00, 0.00000e+00, -2.08841e-02, -1.79873e+00, 1.45103e-03,
2.81769e-04, -1.44703e-03, -5.16394e-05, 8.47001e-02, 1.70147e-01,
5.72562e-03, 5.07493e-05, 4.36148e-03, 1.17863e-04, 4.74364e-03,
6.61278e-03, 4.34292e-05, 1.44373e-03, 2.41470e-05, 2.84426e-03,
8.56560e-04, 2.04028e-03, 0.00000e+00, -3.15994e+03, -2.46423e-03,
1.13843e-03, 4.20512e-04, 0.00000e+00, -9.77214e+01, 6.77794e-03,
5.27499e-03, 1.14936e-03, 0.00000e+00, -6.61311e-03, -1.84255e-02,
-1.96259e-02, 2.98618e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
6.44574e+02, 8.84668e-04, 5.05066e-04, 0.00000e+00, 4.02881e+03,
-1.89503e-03, 0.00000e+00, 0.00000e+00, 8.21407e-04, 2.06780e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-1.20410e-02, -3.63963e-03, 9.92070e-05, -1.15284e-04, -6.33059e-05,
-6.05545e-01, 8.34218e-03, -9.13036e+01, 3.71042e-04, 0.00000e+00,
4.19000e-04, 2.70928e-03, 3.31507e-03, -4.44508e-03, -4.96334e-03,
-1.60449e-03, 3.95119e-03, 2.48924e-03, 5.09815e-04, 4.05302e-03,
2.24076e-03, 0.00000e+00, 6.84256e-03, 4.66354e-04, 0.00000e+00,
-3.68328e-04, 0.00000e+00, 0.00000e+00, -1.46870e+02, 0.00000e+00,
0.00000e+00, 1.09501e-03, 4.65156e-04, 5.62583e-04, 3.21596e+00,
6.43168e-04, 3.14860e-03, 3.40738e-03, 1.78481e-03, 9.62532e-04,
5.58171e-04, 3.43731e+00, -2.33195e-01, 5.10289e-04, 0.00000e+00,
0.00000e+00, -9.25347e+04, 0.00000e+00, -1.99639e-03, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
};
/** NRLMSISE-00 data: density pd[9][150]. */
private static final double[][] PD = {
// HE DENSITY
{
1.09979e+00, -4.88060e-02, -1.97501e-01, -9.10280e-02, -6.96558e-03,
2.42136e-02, 3.91333e-01, -7.20068e-03, -3.22718e-02, 1.41508e+00,
1.68194e-01, 1.85282e-02, 1.09384e-01, -7.24282e+00, 0.00000e+00,
2.96377e-01, -4.97210e-02, 1.04114e+02, -8.61108e-02, -7.29177e-04,
1.48998e-06, 1.08629e-03, 0.00000e+00, 0.00000e+00, 8.31090e-02,
1.12818e-01, -5.75005e-02, -1.29919e-02, -1.78849e-02, -2.86343e-06,
0.00000e+00, -1.51187e+02, -6.65902e-03, 0.00000e+00, -2.02069e-03,
0.00000e+00, 0.00000e+00, 4.32264e-02, -2.80444e+01, -3.26789e-03,
2.47461e-03, 0.00000e+00, 0.00000e+00, 9.82100e-02, 1.22714e-01,
-3.96450e-02, 0.00000e+00, -2.76489e-03, 0.00000e+00, 1.87723e-03,
-8.09813e-03, 4.34428e-05, -7.70932e-03, 0.00000e+00, -2.28894e-03,
-5.69070e-03, -5.22193e-03, 6.00692e-03, -7.80434e+03, -3.48336e-03,
-6.38362e-03, -1.82190e-03, 0.00000e+00, -7.58976e+01, -2.17875e-02,
-1.72524e-02, -9.06287e-03, 0.00000e+00, 2.44725e-02, 8.66040e-02,
1.05712e-01, 3.02543e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-6.01364e+03, -5.64668e-03, -2.54157e-03, 0.00000e+00, 3.15611e+02,
-5.69158e-03, 0.00000e+00, 0.00000e+00, -4.47216e-03, -4.49523e-03,
4.64428e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
4.51236e-02, 2.46520e-02, 6.17794e-03, 0.00000e+00, 0.00000e+00,
-3.62944e-01, -4.80022e-02, -7.57230e+01, -1.99656e-03, 0.00000e+00,
-5.18780e-03, -1.73990e-02, -9.03485e-03, 7.48465e-03, 1.53267e-02,
1.06296e-02, 1.18655e-02, 2.55569e-03, 1.69020e-03, 3.51936e-02,
-1.81242e-02, 0.00000e+00, -1.00529e-01, -5.10574e-03, 0.00000e+00,
2.10228e-03, 0.00000e+00, 0.00000e+00, -1.73255e+02, 5.07833e-01,
-2.41408e-01, 8.75414e-03, 2.77527e-03, -8.90353e-05, -5.25148e+00,
-5.83899e-03, -2.09122e-02, -9.63530e-03, 9.77164e-03, 4.07051e-03,
2.53555e-04, -5.52875e+00, -3.55993e-01, -2.49231e-03, 0.00000e+00,
0.00000e+00, 2.86026e+01, 0.00000e+00, 3.42722e-04, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// O DENSITY
{
1.02315e+00, -1.59710e-01, -1.06630e-01, -1.77074e-02, -4.42726e-03,
3.44803e-02, 4.45613e-02, -3.33751e-02, -5.73598e-02, 3.50360e-01,
6.33053e-02, 2.16221e-02, 5.42577e-02, -5.74193e+00, 0.00000e+00,
1.90891e-01, -1.39194e-02, 1.01102e+02, 8.16363e-02, 1.33717e-04,
6.54403e-06, 3.10295e-03, 0.00000e+00, 0.00000e+00, 5.38205e-02,
1.23910e-01, -1.39831e-02, 0.00000e+00, 0.00000e+00, -3.95915e-06,
0.00000e+00, -7.14651e-01, -5.01027e-03, 0.00000e+00, -3.24756e-03,
0.00000e+00, 0.00000e+00, 4.42173e-02, -1.31598e+01, -3.15626e-03,
1.24574e-03, -1.47626e-03, -1.55461e-03, 6.40682e-02, 1.34898e-01,
-2.42415e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 6.13666e-04,
-5.40373e-03, 2.61635e-05, -3.33012e-03, 0.00000e+00, -3.08101e-03,
-2.42679e-03, -3.36086e-03, 0.00000e+00, -1.18979e+03, -5.04738e-02,
-2.61547e-03, -1.03132e-03, 1.91583e-04, -8.38132e+01, -1.40517e-02,
-1.14167e-02, -4.08012e-03, 1.73522e-04, -1.39644e-02, -6.64128e-02,
-6.85152e-02, -1.34414e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
6.07916e+02, -4.12220e-03, -2.20996e-03, 0.00000e+00, 1.70277e+03,
-4.63015e-03, 0.00000e+00, 0.00000e+00, -2.25360e-03, -2.96204e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
3.92786e-02, 1.31186e-02, -1.78086e-03, 0.00000e+00, 0.00000e+00,
-3.90083e-01, -2.84741e-02, -7.78400e+01, -1.02601e-03, 0.00000e+00,
-7.26485e-04, -5.42181e-03, -5.59305e-03, 1.22825e-02, 1.23868e-02,
6.68835e-03, -1.03303e-02, -9.51903e-03, 2.70021e-04, -2.57084e-02,
-1.32430e-02, 0.00000e+00, -3.81000e-02, -3.16810e-03, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -9.05762e-04, -2.14590e-03, -1.17824e-03, 3.66732e+00,
-3.79729e-04, -6.13966e-03, -5.09082e-03, -1.96332e-03, -3.08280e-03,
-9.75222e-04, 4.03315e+00, -2.52710e-01, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// N2 DENSITY
{
1.16112e+00, 0.00000e+00, 0.00000e+00, 3.33725e-02, 0.00000e+00,
3.48637e-02, -5.44368e-03, 0.00000e+00, -6.73940e-02, 1.74754e-01,
0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
1.26733e-01, 0.00000e+00, 1.03154e+02, 5.52075e-02, 0.00000e+00,
0.00000e+00, 8.13525e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -2.50482e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.48894e-03,
6.16053e-04, -5.79716e-04, 2.95482e-03, 8.47001e-02, 1.70147e-01,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// TOTAL MASS
{
9.44846e-01, 0.00000e+00, 0.00000e+00, -3.08617e-02, 0.00000e+00,
-2.44019e-02, 6.48607e-03, 0.00000e+00, 3.08181e-02, 4.59392e-02,
0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
2.13260e-02, 0.00000e+00, -3.56958e+02, 0.00000e+00, 1.82278e-04,
0.00000e+00, 3.07472e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 3.83054e-03, 0.00000e+00, 0.00000e+00,
-1.93065e-03, -1.45090e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -1.23493e-03, 1.36736e-03, 8.47001e-02, 1.70147e-01,
3.71469e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
5.10250e-03, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 3.68756e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// O2 DENSITY
{
1.35580e+00, 1.44816e-01, 0.00000e+00, 6.07767e-02, 0.00000e+00,
2.94777e-02, 7.46900e-02, 0.00000e+00, -9.23822e-02, 8.57342e-02,
0.00000e+00, 0.00000e+00, 0.00000e+00, 2.38636e+01, 0.00000e+00,
7.71653e-02, 0.00000e+00, 8.18751e+01, 1.87736e-02, 0.00000e+00,
0.00000e+00, 1.49667e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -3.67874e+02, 5.48158e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
1.22631e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
8.17187e-03, 3.71617e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.10826e-03,
-3.13640e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-7.35742e-02, -5.00266e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 1.94965e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// AR DENSITY
{
1.04761e+00, 2.00165e-01, 2.37697e-01, 3.68552e-02, 0.00000e+00,
3.57202e-02, -2.14075e-01, 0.00000e+00, -1.08018e-01, -3.73981e-01,
0.00000e+00, 3.10022e-02, -1.16305e-03, -2.07596e+01, 0.00000e+00,
8.64502e-02, 0.00000e+00, 9.74908e+01, 5.16707e-02, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 3.46193e+02, 1.34297e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.48509e-03,
-1.54689e-04, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
1.47753e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
1.89320e-02, 3.68181e-05, 1.32570e-02, 0.00000e+00, 0.00000e+00,
3.59719e-03, 7.44328e-03, -1.00023e-03, -6.50528e+03, 0.00000e+00,
1.03485e-02, -1.00983e-03, -4.06916e-03, -6.60864e+01, -1.71533e-02,
1.10605e-02, 1.20300e-02, -5.20034e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-2.62769e+03, 7.13755e-03, 4.17999e-03, 0.00000e+00, 1.25910e+04,
0.00000e+00, 0.00000e+00, 0.00000e+00, -2.23595e-03, 4.60217e-03,
5.71794e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-3.18353e-02, -2.35526e-02, -1.36189e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.03522e-02, -6.67837e+01, -1.09724e-03, 0.00000e+00,
-1.38821e-02, 1.60468e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.51574e-02,
-5.44470e-04, 0.00000e+00, 7.28224e-02, 6.59413e-02, 0.00000e+00,
-5.15692e-03, 0.00000e+00, 0.00000e+00, -3.70367e+03, 0.00000e+00,
0.00000e+00, 1.36131e-02, 5.38153e-03, 0.00000e+00, 4.76285e+00,
-1.75677e-02, 2.26301e-02, 0.00000e+00, 1.76631e-02, 4.77162e-03,
0.00000e+00, 5.39354e+00, 0.00000e+00, -7.51710e-03, 0.00000e+00,
0.00000e+00, -8.82736e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// H DENSITY
{
1.26376e+00, -2.14304e-01, -1.49984e-01, 2.30404e-01, 2.98237e-02,
2.68673e-02, 2.96228e-01, 2.21900e-02, -2.07655e-02, 4.52506e-01,
1.20105e-01, 3.24420e-02, 4.24816e-02, -9.14313e+00, 0.00000e+00,
2.47178e-02, -2.88229e-02, 8.12805e+01, 5.10380e-02, -5.80611e-03,
2.51236e-05, -1.24083e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,
1.58727e-01, -3.48190e-02, 0.00000e+00, 0.00000e+00, 2.89885e-05,
0.00000e+00, 1.53595e+02, -1.68604e-02, 0.00000e+00, 1.01015e-02,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.84552e-04,
-1.22181e-03, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
-1.04927e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, -5.91313e-03,
-2.30501e-02, 3.14758e-05, 0.00000e+00, 0.00000e+00, 1.26956e-02,
8.35489e-03, 3.10513e-04, 0.00000e+00, 3.42119e+03, -2.45017e-03,
-4.27154e-04, 5.45152e-04, 1.89896e-03, 2.89121e+01, -6.49973e-03,
-1.93855e-02, -1.48492e-02, 0.00000e+00, -5.10576e-02, 7.87306e-02,
9.51981e-02, -1.49422e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
2.65503e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 6.37110e-03, 3.24789e-04,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
6.14274e-02, 1.00376e-02, -8.41083e-04, 0.00000e+00, 0.00000e+00,
0.00000e+00, -1.27099e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-3.94077e-03, -1.28601e-02, -7.97616e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -6.71465e-03, -1.69799e-03, 1.93772e-03, 3.81140e+00,
-7.79290e-03, -1.82589e-02, -1.25860e-02, -1.04311e-02, -3.02465e-03,
2.43063e-03, 3.63237e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// N DENSITY
{
7.09557e+01, -3.26740e-01, 0.00000e+00, -5.16829e-01, -1.71664e-03,
9.09310e-02, -6.71500e-01, -1.47771e-01, -9.27471e-02, -2.30862e-01,
-1.56410e-01, 1.34455e-02, -1.19717e-01, 2.52151e+00, 0.00000e+00,
-2.41582e-01, 5.92939e-02, 4.39756e+00, 9.15280e-02, 4.41292e-03,
0.00000e+00, 8.66807e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,
1.58727e-01, 9.74701e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 6.70217e+01, -1.31660e-03, 0.00000e+00, -1.65317e-02,
0.00000e+00, 0.00000e+00, 8.50247e-02, 2.77428e+01, 4.98658e-03,
6.15115e-03, 9.50156e-03, -2.12723e-02, 8.47001e-02, 1.70147e-01,
-2.38645e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.37380e-03,
-8.41918e-03, 2.80145e-05, 7.12383e-03, 0.00000e+00, -1.66209e-02,
1.03533e-04, -1.68898e-02, 0.00000e+00, 3.64526e+03, 0.00000e+00,
6.54077e-03, 3.69130e-04, 9.94419e-04, 8.42803e+01, -1.16124e-02,
-7.74414e-03, -1.68844e-03, 1.42809e-03, -1.92955e-03, 1.17225e-01,
-2.41512e-02, 1.50521e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
1.60261e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, -3.54403e-04, -1.87270e-02,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
2.76439e-02, 6.43207e-03, -3.54300e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, -2.80221e-02, 8.11228e+01, -6.75255e-04, 0.00000e+00,
-1.05162e-02, -3.48292e-03, -6.97321e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -1.45546e-03, -1.31970e-02, -3.57751e-03, -1.09021e+00,
-1.50181e-02, -7.12841e-03, -6.64590e-03, -3.52610e-03, -1.87773e-02,
-2.22432e-03, -3.93895e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
// HOT O DENSITY
{
6.04050e-02, 1.57034e+00, 2.99387e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.51018e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, -8.61650e+00, 1.26454e-02,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 5.50878e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 6.23881e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
-9.45934e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
}
};
/** NRLMSISE-00 data: ps[150]. */
private static final double[] PS = {
9.56827e-01, 6.20637e-02, 3.18433e-02, 0.00000e+00, 0.00000e+00,
3.94900e-02, 0.00000e+00, 0.00000e+00, -9.24882e-03, -7.94023e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.74677e-03, 0.00000e+00, 1.54951e-02, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, -6.99007e-04, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 1.24362e-02, -5.28756e-03, 8.47001e-02, 1.70147e-01,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
};
/** NRLMSISE-00 data: TURBO pdl[2][25]. */
private static final double[][] PDL = {
{
1.09930e+00, 3.90631e+00, 3.07165e+00, 9.86161e-01, 1.63536e+01,
4.63830e+00, 1.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 1.28840e+00, 3.10302e-02, 1.18339e-01
},
{
1.00000e+00, 7.00000e-01, 1.15020e+00, 3.44689e+00, 1.28840e+00,
1.00000e+00, 1.08738e+00, 1.22947e+00, 1.10016e+00, 7.34129e-01,
1.15241e+00, 2.22784e+00, 7.95046e-01, 4.01612e+00, 4.47749e+00,
1.23435e+02, -7.60535e-02, 1.68986e-06, 7.44294e-01, 1.03604e+00,
1.72783e+02, 1.15020e+00, 3.44689e+00, -7.46230e-01, 9.49154e-01
}
};
/** NRLMSISE-00 data: LOWER BOUNDARY ptm[10]. */
private static final double[] PTM = {
1.04130e+03, 3.86000e+02, 1.95000e+02, 1.66728e+01, 2.13000e+02,
1.20000e+02, 2.40000e+02, 1.87000e+02, -2.00000e+00, 0.00000e+00
};
/** NRLMSISE-00 data: pdm[8][10]. */
private static final double[][] PDM = {
{
2.45600e+07, 6.71072e-06, 1.00000e+02, 0.00000e+00, 1.10000e+02,
1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
{
8.59400E+10, 1.00000e+00, 1.05000e+02, -8.00000e+00, 1.10000e+02,
1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00
},
{
2.81000E+11, 0.00000e+00, 1.05000e+02, 2.80000e+01, 2.89500e+01,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
{
3.30000E+10, 2.68270e-01, 1.05000e+02, 1.00000e+00, 1.10000e+02,
1.00000e+01, 1.10000e+02, -1.00000e+01, 0.00000e+00, 0.00000e+00
},
{
1.33000e+09, 1.19615e-02, 1.05000e+02, 0.00000e+00, 1.10000e+02,
1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
},
{
1.76100e+05, 1.00000e+00, 9.50000e+01, -8.00000e+00, 1.10000e+02,
1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00,
},
{
1.00000e+07, 1.00000e+00, 1.05000e+02, -8.00000e+00, 1.10000e+02,
1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00
},
{
1.00000e+06, 1.00000e+00, 1.05000e+02, -8.00000e+00, 5.50000e+02,
7.60000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 4.00000e+03
}
};
/** NRLMSISE-00 data: ptl[4][100]. */
private static final double[][] PTL = {
// TN1(2)
{
1.00858e+00, 4.56011e-02, -2.22972e-02, -5.44388e-02, 5.23136e-04,
-1.88849e-02, 5.23707e-02, -9.43646e-03, 6.31707e-03, -7.80460e-02,
-4.88430e-02, 0.00000e+00, 0.00000e+00, -7.60250e+00, 0.00000e+00,
-1.44635e-02, -1.76843e-02, -1.21517e+02, 2.85647e-02, 0.00000e+00,
0.00000e+00, 6.31792e-04, 0.00000e+00, 5.77197e-03, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -8.90272e+03, 3.30611e-03, 3.02172e-03, 0.00000e+00,
-2.13673e-03, -3.20910e-04, 0.00000e+00, 0.00000e+00, 2.76034e-03,
2.82487e-03, -2.97592e-04, -4.21534e-03, 8.47001e-02, 1.70147e-01,
8.96456e-03, 0.00000e+00, -1.08596e-02, 0.00000e+00, 0.00000e+00,
5.57917e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 9.65405e-03, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN1(3)
{
9.39664e-01, 8.56514e-02, -6.79989e-03, 2.65929e-02, -4.74283e-03,
1.21855e-02, -2.14905e-02, 6.49651e-03, -2.05477e-02, -4.24952e-02,
0.00000e+00, 0.00000e+00, 0.00000e+00, 1.19148e+01, 0.00000e+00,
1.18777e-02, -7.28230e-02, -8.15965e+01, 1.73887e-02, 0.00000e+00,
0.00000e+00, 0.00000e+00, -1.44691e-02, 2.80259e-04, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.16584e+02, 3.18713e-03, 7.37479e-03, 0.00000e+00,
-2.55018e-03, -3.92806e-03, 0.00000e+00, 0.00000e+00, -2.89757e-03,
-1.33549e-03, 1.02661e-03, 3.53775e-04, 8.47001e-02, 1.70147e-01,
-9.17497e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
3.56082e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -1.00902e-02, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN1(4)
{
9.85982e-01, -4.55435e-02, 1.21106e-02, 2.04127e-02, -2.40836e-03,
1.11383e-02, -4.51926e-02, 1.35074e-02, -6.54139e-03, 1.15275e-01,
1.28247e-01, 0.00000e+00, 0.00000e+00, -5.30705e+00, 0.00000e+00,
-3.79332e-02, -6.24741e-02, 7.71062e-01, 2.96315e-02, 0.00000e+00,
0.00000e+00, 0.00000e+00, 6.81051e-03, -4.34767e-03, 8.66784e-02,
1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 1.07003e+01, -2.76907e-03, 4.32474e-04, 0.00000e+00,
1.31497e-03, -6.47517e-04, 0.00000e+00, -2.20621e+01, -1.10804e-03,
-8.09338e-04, 4.18184e-04, 4.29650e-03, 8.47001e-02, 1.70147e-01,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-4.04337e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -9.52550e-04,
8.56253e-04, 4.33114e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.21223e-03,
2.38694e-04, 9.15245e-04, 1.28385e-03, 8.67668e-04, -5.61425e-06,
1.04445e+00, 3.41112e+01, 0.00000e+00, -8.40704e-01, -2.39639e+02,
7.06668e-01, -2.05873e+01, -3.63696e-01, 2.39245e+01, 0.00000e+00,
-1.06657e-03, -7.67292e-04, 1.54534e-04, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN1(5) TN2(1)
{
1.00320e+00, 3.83501e-02, -2.38983e-03, 2.83950e-03, 4.20956e-03,
5.86619e-04, 2.19054e-02, -1.00946e-02, -3.50259e-03, 4.17392e-02,
-8.44404e-03, 0.00000e+00, 0.00000e+00, 4.96949e+00, 0.00000e+00,
-7.06478e-03, -1.46494e-02, 3.13258e+01, -1.86493e-03, 0.00000e+00,
-1.67499e-02, 0.00000e+00, 0.00000e+00, 5.12686e-04, 8.66784e-02,
1.58727e-01, -4.64167e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
4.37353e-03, -1.99069e+02, 0.00000e+00, -5.34884e-03, 0.00000e+00,
1.62458e-03, 2.93016e-03, 2.67926e-03, 5.90449e+02, 0.00000e+00,
0.00000e+00, -1.17266e-03, -3.58890e-04, 8.47001e-02, 1.70147e-01,
0.00000e+00, 0.00000e+00, 1.38673e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.60571e-03,
6.28078e-04, 5.05469e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.57829e-03,
-4.00855e-04, 5.04077e-05, -1.39001e-03, -2.33406e-03, -4.81197e-04,
1.46758e+00, 6.20332e+00, 0.00000e+00, 3.66476e-01, -6.19760e+01,
3.09198e-01, -1.98999e+01, 0.00000e+00, -3.29933e+02, 0.00000e+00,
-1.10080e-03, -9.39310e-05, 1.39638e-04, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
}
};
/** NRLMSISE-00 data: pma[10][100]. */
private static final double[][] PMA = {
// TN2(2)
{
9.81637e-01, -1.41317e-03, 3.87323e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.58707e-02,
-8.63658e-03, 0.00000e+00, 0.00000e+00, -2.02226e+00, 0.00000e+00,
-8.69424e-03, -1.91397e-02, 8.76779e+01, 4.52188e-03, 0.00000e+00,
2.23760e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -7.07572e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
-4.11210e-03, 3.50060e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, -8.36657e-03, 1.61347e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, -1.45130e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.24152e-03,
6.43365e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.33255e-03,
2.42657e-03, 1.60666e-03, -1.85728e-03, -1.46874e-03, -4.79163e-06,
1.22464e+00, 3.53510e+01, 0.00000e+00, 4.49223e-01, -4.77466e+01,
4.70681e-01, 8.41861e+00, -2.88198e-01, 1.67854e+02, 0.00000e+00,
7.11493e-04, 6.05601e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN2(3)
{
1.00422e+00, -7.11212e-03, 5.24480e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -5.28914e-02,
-2.41301e-02, 0.00000e+00, 0.00000e+00, -2.12219e+01, -1.03830e-02,
-3.28077e-03, 1.65727e-02, 1.68564e+00, -6.68154e-03, 0.00000e+00,
1.45155e-02, 0.00000e+00, 8.42365e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, -4.34645e-03, 0.00000e+00, 0.00000e+00, 2.16780e-02,
0.00000e+00, -1.38459e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 7.04573e-03, -4.73204e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 1.08767e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -8.08279e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.21769e-04,
-2.27387e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.26769e-03,
3.16901e-03, 4.60316e-04, -1.01431e-04, 1.02131e-03, 9.96601e-04,
1.25707e+00, 2.50114e+01, 0.00000e+00, 4.24472e-01, -2.77655e+01,
3.44625e-01, 2.75412e+01, 0.00000e+00, 7.94251e+02, 0.00000e+00,
2.45835e-03, 1.38871e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN2(4) TN3(1)
{
1.01890e+00, -2.46603e-02, 1.00078e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -6.70977e-02,
-4.02286e-02, 0.00000e+00, 0.00000e+00, -2.29466e+01, -7.47019e-03,
2.26580e-03, 2.63931e-02, 3.72625e+01, -6.39041e-03, 0.00000e+00,
9.58383e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -1.85291e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 1.39717e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 9.19771e-03, -3.69121e+02, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, -1.57067e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.07265e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.92953e-03,
-2.77739e-03, -4.40092e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.47280e-03,
2.95035e-04, -1.81246e-03, 2.81945e-03, 4.27296e-03, 9.78863e-04,
1.40545e+00, -6.19173e+00, 0.00000e+00, 0.00000e+00, -7.93632e+01,
4.44643e-01, -4.03085e+02, 0.00000e+00, 1.15603e+01, 0.00000e+00,
2.25068e-03, 8.48557e-04, -2.98493e-04, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN3(2)
{
9.75801e-01, 3.80680e-02, -3.05198e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.85575e-02,
5.04057e-02, 0.00000e+00, 0.00000e+00, -1.76046e+02, 1.44594e-02,
-1.48297e-03, -3.68560e-03, 3.02185e+01, -3.23338e-03, 0.00000e+00,
1.53569e-02, 0.00000e+00, -1.15558e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 4.89620e-03, 0.00000e+00, 0.00000e+00, -1.00616e-02,
-8.21324e-03, -1.57757e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 6.63564e-03, 4.58410e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, -2.51280e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 9.91215e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -8.73148e-04,
-1.29648e-03, -7.32026e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -4.68110e-03,
-4.66003e-03, -1.31567e-03, -7.39390e-04, 6.32499e-04, -4.65588e-04,
-1.29785e+00, -1.57139e+02, 0.00000e+00, 2.58350e-01, -3.69453e+01,
4.10672e-01, 9.78196e+00, -1.52064e-01, -3.85084e+03, 0.00000e+00,
-8.52706e-04, -1.40945e-03, -7.26786e-04, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN3(3)
{
9.60722e-01, 7.03757e-02, -3.00266e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.22671e-02,
4.10423e-02, 0.00000e+00, 0.00000e+00, -1.63070e+02, 1.06073e-02,
5.40747e-04, 7.79481e-03, 1.44908e+02, 1.51484e-04, 0.00000e+00,
1.97547e-02, 0.00000e+00, -1.41844e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 5.77884e-03, 0.00000e+00, 0.00000e+00, 9.74319e-03,
0.00000e+00, -2.88015e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, -4.44902e-03, -2.92760e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 2.34419e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.36685e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -4.65325e-04,
-5.50628e-04, 3.31465e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.06179e-03,
-3.08575e-03, -7.93589e-04, -1.08629e-04, 5.95511e-04, -9.05050e-04,
1.18997e+00, 4.15924e+01, 0.00000e+00, -4.72064e-01, -9.47150e+02,
3.98723e-01, 1.98304e+01, 0.00000e+00, 3.73219e+03, 0.00000e+00,
-1.50040e-03, -1.14933e-03, -1.56769e-04, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN3(4)
{
1.03123e+00, -7.05124e-02, 8.71615e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.82621e-02,
-9.80975e-03, 0.00000e+00, 0.00000e+00, 2.89286e+01, 9.57341e-03,
0.00000e+00, 0.00000e+00, 8.66153e+01, 7.91938e-04, 0.00000e+00,
0.00000e+00, 0.00000e+00, 4.68917e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 7.86638e-03, 0.00000e+00, 0.00000e+00, 9.90827e-03,
0.00000e+00, 6.55573e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, -4.00200e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 7.07457e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.72268e-03,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.04970e-04,
1.21560e-03, -8.05579e-06, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.49941e-03,
-4.57256e-04, -1.59311e-04, 2.96481e-04, -1.77318e-03, -6.37918e-04,
1.02395e+00, 1.28172e+01, 0.00000e+00, 1.49903e-01, -2.63818e+01,
0.00000e+00, 4.70628e+01, -2.22139e-01, 4.82292e-02, 0.00000e+00,
-8.67075e-04, -5.86479e-04, 5.32462e-04, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TN3(5) SURFACE TEMP TSL
{
1.00828e+00, -9.10404e-02, -2.26549e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.32420e-02,
-9.08925e-03, 0.00000e+00, 0.00000e+00, 3.36105e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, -1.24957e+01, -5.87939e-03, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.79765e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 2.01237e+03, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, -1.75553e-02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.29699e-03,
1.26659e-03, 2.68402e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.17894e-03,
1.48746e-03, 1.06478e-04, 1.34743e-04, -2.20939e-03, -6.23523e-04,
6.36539e-01, 1.13621e+01, 0.00000e+00, -3.93777e-01, 2.38687e+03,
0.00000e+00, 6.61865e+02, -1.21434e-01, 9.27608e+00, 0.00000e+00,
1.68478e-04, 1.24892e-03, 1.71345e-03, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TGN3(2) SURFACE GRAD TSLG
{
1.57293e+00, -6.78400e-01, 6.47500e-01, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.62974e-02,
-3.60423e-01, 0.00000e+00, 0.00000e+00, 1.28358e+02, 0.00000e+00,
0.00000e+00, 0.00000e+00, 4.68038e+01, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -1.67898e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.90994e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 3.15706e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TGN2(1) TGN1(2)
{
8.60028e-01, 3.77052e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.17570e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 7.77757e-03, 0.00000e+00,
0.00000e+00, 0.00000e+00, 1.01024e+02, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 6.54251e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.56959e-02,
1.91001e-02, 3.15971e-02, 1.00982e-02, -6.71565e-03, 2.57693e-03,
1.38692e+00, 2.82132e-01, 0.00000e+00, 0.00000e+00, 3.81511e+02,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
},
// TGN3(1) TGN2(2)
{
1.06029e+00, -5.25231e-02, 3.73034e-01, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.31072e-02,
-3.88409e-01, 0.00000e+00, 0.00000e+00, -1.65295e+02, -2.13801e-01,
-4.38916e-02, -3.22716e-01, -8.82393e+01, 1.18458e-01, 0.00000e+00,
-4.35863e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, -1.19782e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 2.62229e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, -5.37443e+01, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, -4.55788e-01, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.84009e-02,
3.96733e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.05494e-02,
7.39617e-02, 1.92200e-02, -8.46151e-03, -1.34244e-02, 1.96338e-02,
1.50421e+00, 1.88368e+01, 0.00000e+00, 0.00000e+00, -5.13114e+01,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
5.11923e-02, 3.61225e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
}
};
/** NRLMSISE-00 data: MIDDLE ATMOSPHERE AVERAGES pavgm[10]. */
private static final double[] PAVGM = {
2.61000e+02, 2.64000e+02, 2.29000e+02, 2.17000e+02, 2.17000e+02,
2.23000e+02, 2.86760e+02, -2.93940e+00, 2.50000e+00, 0.00000e+00
};
// Fields
/** External data container. */
private final NRLMSISE00InputParameters inputParams;
/** Sun position. */
private PVCoordinatesProvider sun;
/** Earth body shape. */
private final BodyShape earth;
/** Switches for main effects. */
private final int[] sw;
/** Switches for cross effects. */
private final int[] swc;
/** Constructor.
* <p>
* The model is constructed with all switches set to 1.
* </p>
* <p>
* Parameters are mandatory only for the
* {@link #getDensity(AbsoluteDate, Vector3D, Frame) getDensity()} and
* {@link #getVelocity(AbsoluteDate, Vector3D, Frame) getVelocity()} methods.
* </p>
* @param parameters the solar and magnetic activity data
* @param sun the Sun position
* @param earth the Earth body shape
*/
public NRLMSISE00(final NRLMSISE00InputParameters parameters,
final PVCoordinatesProvider sun,
final BodyShape earth) {
this(parameters, sun, earth, allOnes(), allOnes());
}
/** Constructor.
* <p>
* The model is constructed with all switches set to 1.
* </p>
* <p>
* Parameters are mandatory only for the
* {@link #getDensity(AbsoluteDate, Vector3D, Frame) getDensity()} and
* {@link #getVelocity(AbsoluteDate, Vector3D, Frame) getVelocity()} methods.
* </p>
* @param parameters the solar and magnetic activity data
* @param sun the Sun position
* @param earth the Earth body shape
* @param sw switches for main effects
* @param swc switches for cross effects
*/
private NRLMSISE00(final NRLMSISE00InputParameters parameters,
final PVCoordinatesProvider sun,
final BodyShape earth,
final int[] sw,
final int[] swc) {
this.inputParams = parameters;
this.sun = sun;
this.earth = earth;
this.sw = sw;
this.swc = swc;
}
/** Change a switch.
* <p>
* This method creates a new instance, the current instance is
* not changed at all!
* </p>
* @param number switch number between 1 and 23
* @param value switch value
* @return a <em>new</em> instance, with switch changed
* @exception OrekitException if switch number is not between 1 and 23
*/
public NRLMSISE00 withSwitch(final int number, final int value)
throws OrekitException {
if (number < 1 || number > 23) {
throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, number, 1, 23);
}
final int[] newSw = sw.clone();
final int[] newSwc = swc.clone();
if (number != 9) {
newSw[number] = (value == 1) ? 1 : 0;
newSwc[number] = (value > 0) ? 1 : 0;
} else {
if (value == -1 || value == 1) {
newSw[number] = value;
} else {
newSw[number] = 0;
}
newSwc[number] = newSw[number];
}
return new NRLMSISE00(inputParams, sun, earth, newSwc, newSwc);
}
/** Create an array of switches set to 1.
* @return array of switches
*/
private static int[] allOnes() {
final int[] array = new int[24];
Arrays.fill(array, 1);
return array;
}
/** {@inheritDoc} */
@Override
public Frame getFrame() {
return earth.getBodyFrame();
}
/** {@inheritDoc} */
@Override
public double getDensity(final AbsoluteDate date,
final Vector3D position,
final Frame frame)
throws OrekitException {
// check if data are available :
if ((date.compareTo(inputParams.getMaxDate()) > 0) ||
(date.compareTo(inputParams.getMinDate()) < 0)) {
throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
date, inputParams.getMinDate(), inputParams.getMaxDate());
}
// compute day number in current year and the seconds within the day
final DateTimeComponents dtc = date.getComponents(TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
final int doy = dtc.getDate().getDayOfYear();
final double sec = dtc.getTime().getSecondsInLocalDay();
// compute geodetic position (km and °)
final GeodeticPoint inBody = earth.transform(position, frame, date);
final double alt = inBody.getAltitude() / 1000.;
final double lon = FastMath.toDegrees(inBody.getLongitude());
final double lat = FastMath.toDegrees(inBody.getLatitude());
// compute local solar time
final double lst = localSolarTime(date, position, frame);
// get solar activity data and compute
final Output out = new Output(doy, sec, lat, lon, lst, inputParams.getAverageFlux(date),
inputParams.getDailyFlux(date), inputParams.getAp(date));
out.gtd7d(alt);
// return the local density
return out.getDensity(TOTAL_MASS);
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
final FieldVector3D<T> position,
final Frame frame)
throws OrekitException {
// check if data are available :
final AbsoluteDate dateD = date.toAbsoluteDate();
if ((dateD.compareTo(inputParams.getMaxDate()) > 0) ||
(dateD.compareTo(inputParams.getMinDate()) < 0)) {
throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
dateD, inputParams.getMinDate(), inputParams.getMaxDate());
}
// compute day number in current year and the seconds within the day
final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
final DateTimeComponents dtc = dateD.getComponents(ut1);
final int doy = dtc.getDate().getDayOfYear();
final T sec = date.durationFrom(new AbsoluteDate(dtc.getDate(), TimeComponents.H00, ut1));
// compute geodetic position (km and °)
final FieldGeodeticPoint<T> inBody = earth.transform(position, frame, date);
final T alt = inBody.getAltitude().divide(1000.);
final T lon = inBody.getLongitude().multiply(180.0 / FastMath.PI);
final T lat = inBody.getLatitude().multiply(180.0 / FastMath.PI);
// compute local solar time
final T lst = localSolarTime(dateD, position, frame);
// get solar activity data and compute
final FieldOutput<T> out = new FieldOutput<>(doy, sec, lat, lon, lst,
inputParams.getAverageFlux(dateD),
inputParams.getDailyFlux(dateD), inputParams.getAp(dateD));
out.gtd7d(alt);
// return the local density
return out.getDensity(TOTAL_MASS);
}
/** Get local solar time.
* @param date current date
* @param position current position in frame
* @param frame the frame in which is defined the position
* @return the local solar time (hour in [0, 24[)
* @throws OrekitException if position cannot be computed in given frame
*/
private double localSolarTime(final AbsoluteDate date,
final Vector3D position,
final Frame frame) throws OrekitException {
final Vector3D sunPos = sun.getPVCoordinates(date, frame).getPosition();
final double lst = FastMath.PI + FastMath.atan2(
sunPos.getX() * position.getY() - sunPos.getY() * position.getX(),
sunPos.getX() * position.getX() + sunPos.getY() * position.getY());
return lst * 12. / FastMath.PI;
}
/** Get local solar time.
* @param date current date
* @param position current position in frame
* @param frame the frame in which is defined the position
* @param <T> type of the filed elements
* @return the local solar time (hour in [0, 24[)
* @throws OrekitException if position cannot be computed in given frame
*/
private <T extends RealFieldElement<T>> T localSolarTime(final AbsoluteDate date,
final FieldVector3D<T> position,
final Frame frame)
throws OrekitException {
final Vector3D sunPos = sun.getPVCoordinates(date, frame).getPosition();
final T y = position.getY().multiply(sunPos.getX()).subtract(position.getX().multiply(sunPos.getY()));
final T x = position.getX().multiply(sunPos.getX()).add(position.getY().multiply(sunPos.getY()));
final T hl = y.atan2(x).add(FastMath.PI);
return hl.multiply(12. / FastMath.PI);
}
/**
* This class is a placeholder for the computed densities and temperatures.
* <p>
* Densities are provided as an array d such as:
* <ul>
* <li>d[0] = He number density (1/m³)</li>
* <li>d[1] = O number density (1/m³)</li>
* <li>d[2] = N2 number density (1/m³)</li>
* <li>d[3] = O2 number density (1/m³)</li>
* <li>d[4] = Ar number density (1/m³)</li>
* <li>d[5] = total mass density (kg/m³) (*)</li>
* <li>d[6] = H number density (1/m³)</li>
* <li>d[7] = N number density (1/m³)</li>
* <li>d[8] = anomalous oxygen number density (1/m³)
* </ul>
* Total mass density, d[5], is NOT the same for methods gtd7 and gtd7d:
* <ul>
* <li>For gtd7: d[5] is the sum of the mass densities of the species
* He, O, N2, O2, Ar, H and N but does NOT include anomalous oxygen.</li>
* <li>For gtd7d: d[5] is the "effective total mass density for drag" and is the sum
* of the mass densities of all species in this model, INCLUDING anomalous oxygen.</li>
* </ul>
* O, H, and N are set to zero below 72.5 km.
* </p>
* <p>
* Temperatures are provided as an array t such as:
* <ul>
* <li>t[0] = exospheric temperature (K)</li>
* <li>t[1] = temperature at altitude (K)</li>
* </ul>
* t[0] is set to global average for altitudes below 120 km.<br>
* The 120 km gradient is left at global average value for altitudes below 72 km.
* </p>
*/
private class Output {
/** Day of year (from 1 to 365 or 366). */
private final int doy;
/** Seconds in day (UT scale). */
private final double sec;
/** Geodetic latitude (°). */
private final double lat;
/** Geodetic longitude (°). */
private final double lon;
/** Local apparent solar time (hours). */
private final double hl;
/** 81 day average of F10.7 flux (centered on day). */
private final double f107a;
/** Daily F10.7 flux for previous day. */
private final double f107;
/** Array containing:
* <ul>
* <li>0: daily Ap</li>
* <li>1: 3 hr ap index for current time</li>
* <li>2: 3 hr ap index for 3 hrs before current time</li>
* <li>3: 3 hr ap index for 6 hrs before current time</li>
* <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
* <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
* <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
* </ul>. */
private final double[] ap;
/** Gravity at latitude (cm/s2). */
private final double glat;
/** Effective Earth radius at latitude (km). */
private final double rlat;
/** N2 mixed density at alt. */
private double dm28;
/** Legendre polynomials. */
private final double[][] plg;
/** Cosinus of local solar time. */
private final double ctloc;
/** Sinus of local solar time. */
private final double stloc;
/** Square of ctloc. */
private final double c2tloc;
/** Square of stloc. */
private final double s2tloc;
/** Cube of ctloc. */
private final double c3tloc;
/** Cube of stloc. */
private final double s3tloc;
/** Magnetic activity based on daily ap. */
private double apdf;
/** Magnetic activity based on daily ap. */
private double apt;
/** Temperature at nodes for ZN1 scale. */
private final double[] meso_tn1;
/** Temperature at nodes for ZN2 scale. */
private final double[] meso_tn2;
/** Temperature at nodes for ZN3 scale. */
private final double[] meso_tn3;
/** Temperature gradients at end nodes for ZN1 scale. */
private final double[] meso_tgn1;
/** Temperature gradients at end nodes for ZN2 scale. */
private final double[] meso_tgn2;
/** Temperature gradients at end nodes for ZN3 scale. */
private final double[] meso_tgn3;
/** Densities. */
private final double[] densities;
/** Temperatures. */
private final double[] temperatures;
/** Simple constructor.
* @param doy day of year (from 1 to 365 or 366)
* @param sec seconds in day (UT scale)
* @param lat geodetic latitude (°)
* @param lon geodetic longitude (°)
* @param hl local apparent solar time (hours)
* @param f107a 81 day average of F10.7 flux (centered on day)
* @param f107 daily F10.7 flux for previous day
* @param ap array containing:
* <ul>
* <li>0: daily Ap</li>
* <li>1: 3 hr ap index for current time</li>
* <li>2: 3 hr ap index for 3 hrs before current time</li>
* <li>3: 3 hr ap index for 6 hrs before current time</li>
* <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
* <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
* <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
* </ul>
*/
Output(final int doy, final double sec,
final double lat, final double lon, final double hl,
final double f107a, final double f107, final double[] ap) {
this.doy = doy;
this.sec = sec;
this.lat = lat;
this.lon = lon;
this.hl = hl;
this.f107a = f107a;
this.f107 = f107;
this.ap = ap.clone();
this.plg = new double[4][8];
this.meso_tn1 = new double[ZN1.length];
this.meso_tn2 = new double[ZN2.length];
this.meso_tn3 = new double[ZN3.length];
this.meso_tgn1 = new double[2];
this.meso_tgn2 = new double[2];
this.meso_tgn3 = new double[2];
densities = new double[9];
temperatures = new double[2];
// Calculates latitude variable gravity and effective radius
final double xlat = (sw[2] == 0) ? LAT_REF : lat;
final double c2 = FastMath.cos(2 * DEG_TO_RAD * xlat);
glat = G_REF * (1. - .0026373 * c2);
rlat = 2. * glat / (3.085462e-6 + 2.27e-9 * c2) * 1.e-5;
// Convert latitude into radians
final double latr = DEG_TO_RAD * lat;
// Calculate legendre polynomials
final double c = FastMath.sin(latr);
final double s = FastMath.cos(latr);
plg[0][1] = c;
plg[0][2] = ( 3.0 * c * plg[0][1] - 1.0) / 2.0;
plg[0][3] = ( 5.0 * c * plg[0][2] - 2.0 * plg[0][1]) / 3.0;
plg[0][4] = ( 7.0 * c * plg[0][3] - 3.0 * plg[0][2]) / 4.0;
plg[0][5] = ( 9.0 * c * plg[0][4] - 4.0 * plg[0][3]) / 5.0;
plg[0][6] = (11.0 * c * plg[0][5] - 5.0 * plg[0][4]) / 6.0;
plg[1][1] = s;
plg[1][2] = 3.0 * c * plg[1][1];
plg[1][3] = ( 5.0 * c * plg[1][2] - 3.0 * plg[1][1]) / 2.0;
plg[1][4] = ( 7.0 * c * plg[1][3] - 4.0 * plg[1][2]) / 3.0;
plg[1][5] = ( 9.0 * c * plg[1][4] - 5.0 * plg[1][3]) / 4.0;
plg[1][6] = (11.0 * c * plg[1][5] - 6.0 * plg[1][4]) / 5.0;
plg[2][2] = 3.0 * s * plg[1][1];
plg[2][3] = 5.0 * c * plg[2][2];
plg[2][4] = ( 7.0 * c * plg[2][3] - 5.0 * plg[2][2]) / 2.0;
plg[2][5] = ( 9.0 * c * plg[2][4] - 6.0 * plg[2][3]) / 3.0;
plg[2][6] = (11.0 * c * plg[2][5] - 7.0 * plg[2][4]) / 4.0;
plg[2][7] = (13.0 * c * plg[2][6] - 8.0 * plg[2][5]) / 5.0;
plg[3][3] = 5.0 * s * plg[2][2];
plg[3][4] = 7.0 * c * plg[3][3];
plg[3][5] = ( 9.0 * c * plg[3][4] - 7.0 * plg[3][3]) / 2.0;
plg[3][6] = (11.0 * c * plg[3][5] - 8.0 * plg[3][4]) / 3.0;
// Calculate additional data
if (!(sw[7] == 0 && sw[8] == 0 && sw[14] == 0)) {
final double tloc = HOUR_TO_RAD * hl;
final double tlx2 = tloc + tloc;
final double tlx3 = tloc + tlx2;
stloc = FastMath.sin(tloc);
ctloc = FastMath.cos(tloc);
s2tloc = FastMath.sin(tlx2);
c2tloc = FastMath.cos(tlx2);
s3tloc = FastMath.sin(tlx3);
c3tloc = FastMath.cos(tlx3);
} else {
stloc = 0;
ctloc = 0;
s2tloc = 0;
c2tloc = 0;
s3tloc = 0;
c3tloc = 0;
}
}
/** Calculate temperatures and densities not including anomalous oxygen.
* <p>
* This method is the thermospheric portion of NRLMSISE-00 for alt > 72.5 km.
* </p>
* <p>NOTES ON INPUT VARIABLES:<br>
* Seconds, Local Time, and Longitude are used independently in the
* model and are not of equal importance for every situation.<br>
* For the most physically realistic calculation these three
* variables should be consistent (lst=sec/3600 + lon/15).<br>
* The Equation of Time departures from the above formula
* for apparent local time can be included if available but
* are of minor importance.<br><br>
*
* f107 and f107A values used to generate the model correspond
* to the 10.7 cm radio flux at the actual distance of the Earth
* from the Sun rather than the radio flux at 1 AU. The following
* site provides both classes of values:<br>
* ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
*
* f107, f107A, and ap effects are neither large nor well established below 80 km
* and these parameters should be set to 150., 150., and 4. respectively.
* </p>
* @param alt altitude (km)
*/
void gts7(final double alt) {
// Thermal diffusion coefficients for species
final double[] alpha = {-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
// Altitude limits for net density computation for species
final double[] altl = {200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
// N2 mixed density
final double xmm = PDM[2][4];
/**** Exospheric temperature ****/
double tinf = PTM[0] * PT[0];
// Tinf variations not important below ZA or ZN[0]
if (alt > ZN1[0]) {
tinf *= 1.0 + sw[16] * globe7(PT);
}
setTemperature(EXOSPHERIC, tinf);
// Gradient variations not important below ZN[4]
double g0 = PTM[3] * PS[0];
if (alt > ZN1[4]) {
g0 *= 1.0 + sw[19] * globe7(PS);
}
// Temperature at lower boundary
double tlb = PTM[1] * PD[3][0];
tlb *= 1.0 + sw[17] * globe7(PD[3]);
// Slope
final double s = g0 / (tinf - tlb);
// Lower thermosphere temp variations not significant for density above 300 km
meso_tn1[1] = PTM[6] * PTL[0][0];
meso_tn1[2] = PTM[2] * PTL[1][0];
meso_tn1[3] = PTM[7] * PTL[2][0];
meso_tn1[4] = PTM[4] * PTL[3][0];
meso_tgn1[1] = PTM[8] * PMA[8][0];
if (alt < 300.0) {
final double r = PTM[4] * PTL[3][0];
meso_tn1[1] /= 1.0 - sw[18] * glob7s(PTL[0]);
meso_tn1[2] /= 1.0 - sw[18] * glob7s(PTL[1]);
meso_tn1[3] /= 1.0 - sw[18] * glob7s(PTL[2]);
meso_tn1[4] /= 1.0 - sw[18] * sw[20] * glob7s(PTL[3]);
meso_tgn1[1] *= 1.0 + sw[18] * sw[20] * glob7s(PMA[8]);
meso_tgn1[1] *= meso_tn1[4] * meso_tn1[4] / (r * r);
}
/**** Temperature at altitude ****/
setTemperature(ALTITUDE, densu(alt, 1.0, tinf, tlb, 0.0, 0.0, PTM[5], s));
/**** N2 density ****/
/* Density variation factor at Zlb */
final double g28 = sw[21] * globe7(PD[2]);
/* Diffusive density at Zlb */
final double db28 = PDM[2][0] * FastMath.exp(g28) * PD[2][0];
/* Diffusive density at Alt */
double diffusiveDensity = densu(alt, db28, tinf, tlb, N2_MASS, alpha[2], PTM[5], s);
setDensity(MOLECULAR_NITROGEN, diffusiveDensity);
// Variation of turbopause height
final double zhf = PDL[1][24] * (1.0 + sw[5] * PDL[0][24] *
FastMath.sin(DEG_TO_RAD * lat) *
FastMath.cos(DAY_TO_RAD * (doy - PT[13])));
/* Turbopause */
final double zh28 = PDM[2][2] * zhf;
final double zhm28 = PDM[2][3] * PDL[1][5];
/* Mixed density at Zlb */
final double b28 = densu(zh28, db28, tinf, tlb, N2_MASS - xmm, alpha[2] - 1.0, PTM[5], s);
if (sw[15] != 0 && alt <= altl[2]) {
/* Mixed density at Alt */
dm28 = densu(alt, b28, tinf, tlb, xmm, alpha[2], PTM[5], s);
/* Net density at Alt */
setDensity(MOLECULAR_NITROGEN, dnet(diffusiveDensity, dm28, zhm28, xmm, N2_MASS));
}
/**** He density ****/
/* Density variation factor at Zlb */
final double g4 = sw[21] * globe7(PD[0]);
/* Diffusive density at Zlb */
final double db04 = PDM[0][0] * FastMath.exp(g4) * PD[0][0];
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db04, tinf, tlb, HE_MASS, alpha[0], PTM[5], s);
setDensity(HELIUM, diffusiveDensity);
if (sw[15] != 0 && alt < altl[0]) {
/* Turbopause */
final double zh04 = PDM[0][2];
/* Mixed density at Zlb */
final double b04 = densu(zh04, db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
/* Mixed density at Alt */
final double dm04 = densu(alt, b04, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm04 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm04, zhm04, xmm, HE_MASS);
/* Correction to specified mixing ratio at ground */
final double rl = FastMath.log(b28 * PDM[0][1] / b04);
final double zc04 = PDM[0][4] * PDL[1][0];
final double hc04 = PDM[0][5] * PDL[1][1];
/* Net density corrected at Alt */
setDensity(HELIUM, diffusiveDensity * ccor(alt, rl, hc04, zc04));
}
/**** O density ****/
/* Density variation factor at Zlb */
final double g16 = sw[21] * globe7(PD[1]);
/* Diffusive density at Zlb */
final double db16 = PDM[1][0] * FastMath.exp(g16) * PD[1][0];
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db16, tinf, tlb, O_MASS, alpha[1], PTM[5], s);
setDensity(ATOMIC_OXYGEN, diffusiveDensity);
if (sw[15] != 0 && alt < altl[1]) {
/* Turbopause */
final double zh16 = PDM[1][2];
/* Mixed density at Zlb */
final double b16 = densu(zh16, db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
/* Mixed density at Alt */
final double dm16 = densu(alt, b16, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm16 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm16, zhm16, xmm, O_MASS);
final double rl = PDM[1][1] * PDL[1][16] * (1.0 + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
final double hc16 = PDM[1][5] * PDL[1][3];
final double zc16 = PDM[1][4] * PDL[1][2];
final double hc216 = PDM[1][5] * PDL[1][4];
diffusiveDensity *= ccor2(alt, rl, hc16, zc16, hc216);
/* Chemistry correction */
final double hcc16 = PDM[1][7] * PDL[1][13];
final double zcc16 = PDM[1][6] * PDL[1][12];
final double rc16 = PDM[1][3] * PDL[1][14];
/* Net density corrected at Alt */
setDensity(ATOMIC_OXYGEN, diffusiveDensity * ccor(alt, rc16, hcc16, zcc16));
}
/**** O2 density ****/
/* Density variation factor at Zlb */
final double g32 = sw[21] * globe7(PD[4]);
/* Diffusive density at Zlb */
final double db32 = PDM[3][0] * FastMath.exp(g32) * PD[4][0];
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db32, tinf, tlb, O2_MASS, alpha[3], PTM[5], s);
setDensity(MOLECULAR_OXYGEN, diffusiveDensity);
if (sw[15] != 0) {
if (alt <= altl[3]) {
/* Turbopause */
final double zh32 = PDM[3][2];
/* Mixed density at Zlb */
final double b32 = densu(zh32, db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
/* Mixed density at Alt */
final double dm32 = densu(alt, b32, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm32 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm32, zhm32, xmm, O2_MASS);
/* Correction to specified mixing ratio at ground */
final double rl = FastMath.log(b28 * PDM[3][1] / b32);
final double hc32 = PDM[3][5] * PDL[1][7];
final double zc32 = PDM[3][4] * PDL[1][6];
diffusiveDensity *= ccor(alt, rl, hc32, zc32);
}
/* Correction for general departure from diffusive equilibrium above Zlb */
final double hcc32 = PDM[3][7] * PDL[1][22];
final double hcc232 = PDM[3][7] * PDL[0][22];
final double zcc32 = PDM[3][6] * PDL[1][21];
final double rc32 = PDM[3][3] * PDL[1][23] * (1. + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
/* Net density corrected at Alt */
setDensity(MOLECULAR_OXYGEN, diffusiveDensity * ccor2(alt, rc32, hcc32, zcc32, hcc232));
}
/**** Ar density ****/
/* Density variation factor at Zlb */
final double g40 = sw[21] * globe7(PD[5]);
/* Diffusive density at Zlb */
final double db40 = PDM[4][0] * FastMath.exp(g40) * PD[5][0];
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db40, tinf, tlb, AR_MASS, alpha[4], PTM[5], s);
setDensity(ARGON, diffusiveDensity);
if (sw[15] != 0 && alt <= altl[4]) {
/* Turbopause */
final double zh40 = PDM[4][2];
/* Mixed density at Zlb */
final double b40 = densu(zh40, db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
/* Mixed density at Alt */
final double dm40 = densu(alt, b40, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm40 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm40, zhm40, xmm, AR_MASS);
/* Correction to specified mixing ratio at ground */
final double rl = FastMath.log(b28 * PDM[4][1] / b40);
final double hc40 = PDM[4][5] * PDL[1][9];
final double zc40 = PDM[4][4] * PDL[1][8];
/* Net density corrected at Alt */
setDensity(ARGON, diffusiveDensity * ccor(alt, rl, hc40, zc40));
}
/**** H density ****/
/* Density variation factor at Zlb */
final double g1 = sw[21] * globe7(PD[6]);
/* Diffusive density at Zlb */
final double db01 = PDM[5][0] * FastMath.exp(g1) * PD[6][0];
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db01, tinf, tlb, H_MASS, alpha[6], PTM[5], s);
setDensity(HYDROGEN, diffusiveDensity);
if (sw[15] != 0 && alt <= altl[6]) {
/* Turbopause */
final double zh01 = PDM[5][2];
/* Mixed density at Zlb */
final double b01 = densu(zh01, db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
/* Mixed density at Alt */
final double dm01 = densu(alt, b01, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm01 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm01, zhm01, xmm, H_MASS);
/* Correction to specified mixing ratio at ground */
final double rl = FastMath.log(b28 * PDM[5][1] * FastMath.sqrt(PDL[1][17] * PDL[1][17]) / b01);
final double hc01 = PDM[5][5] * PDL[1][11];
final double zc01 = PDM[5][4] * PDL[1][10];
diffusiveDensity *= ccor(alt, rl, hc01, zc01);
/* Chemistry correction */
final double hcc01 = PDM[5][7] * PDL[1][19];
final double zcc01 = PDM[5][6] * PDL[1][18];
final double rc01 = PDM[5][3] * PDL[1][20];
/* Net density corrected at Alt */
setDensity(HYDROGEN, diffusiveDensity * ccor(alt, rc01, hcc01, zcc01));
}
/**** N density ****/
/* Density variation factor at Zlb */
final double g14 = sw[21] * globe7(PD[7]);
/* Diffusive density at Zlb */
final double db14 = PDM[6][0] * FastMath.exp(g14) * PD[7][0];
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db14, tinf, tlb, N_MASS, alpha[7], PTM[5], s);
setDensity(ATOMIC_NITROGEN, diffusiveDensity);
if (sw[15] != 0 && alt <= altl[7]) {
/* Turbopause */
final double zh14 = PDM[6][2];
/* Mixed density at Zlb */
final double b14 = densu(zh14, db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
/* Mixed density at Alt */
final double dm14 = densu(alt, b14, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm14 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm14, zhm14, xmm, N_MASS);
/* Correction to specified mixing ratio at ground */
final double rl = FastMath.log(b28 * PDM[6][1] * PDL[0][2] / b14);
final double hc14 = PDM[6][5] * PDL[0][1];
final double zc14 = PDM[6][4] * PDL[0][0];
diffusiveDensity *= ccor(alt, rl, hc14, zc14);
/* Chemistry correction */
final double hcc14 = PDM[6][7] * PDL[0][4];
final double zcc14 = PDM[6][6] * PDL[0][3];
final double rc14 = PDM[6][3] * PDL[0][5];
/* Net density corrected at Alt */
setDensity(ATOMIC_NITROGEN, diffusiveDensity * ccor(alt, rc14, hcc14, zcc14));
}
/**** Anomalous O density ****/
final double g16h = sw[21] * globe7(PD[8]);
final double db16h = PDM[7][0] * FastMath.exp(g16h) * PD[8][0];
final double tho = PDM[7][9] * PDL[0][6];
diffusiveDensity = densu(alt, db16h, tho, tho, O_MASS, alpha[8], PTM[5], s);
final double zsht = PDM[7][5];
final double zmho = PDM[7][4];
final double zsho = scalh(zmho, O_MASS, tho);
diffusiveDensity *= FastMath.exp(-zsht / zsho * (FastMath.exp((zmho - alt ) / zsht) - 1.));
setDensity(ANOMALOUS_OXYGEN, diffusiveDensity);
// Convert densities from cm-3 to m-3
for (int i = 0; i < 9; i++) {
setDensity(i, getDensity(i) * 1.0e+06);
}
/**** Total mass density ****/
final double tmd = AMU * (HE_MASS * getDensity(HELIUM) +
O_MASS * getDensity(ATOMIC_OXYGEN) +
N2_MASS * getDensity(MOLECULAR_NITROGEN) +
O2_MASS * getDensity(MOLECULAR_OXYGEN) +
AR_MASS * getDensity(ARGON) +
H_MASS * getDensity(HYDROGEN) +
N_MASS * getDensity(ATOMIC_NITROGEN));
setDensity(TOTAL_MASS, tmd);
}
/** Calculate temperatures and densities not including anomalous oxygen.
* <p>NOTES ON INPUT VARIABLES:<br>
* Seconds, Local Time, and Longitude are used independently in the
* model and are not of equal importance for every situation.<br>
* For the most physically realistic calculation these three
* variables should be consistent (lst=sec/3600 + lon/15).<br>
* The Equation of Time departures from the above formula
* for apparent local time can be included if available but
* are of minor importance.<br><br>
*
* f107 and f107A values used to generate the model correspond
* to the 10.7 cm radio flux at the actual distance of the Earth
* from the Sun rather than the radio flux at 1 AU. The following
* site provides both classes of values:<br>
* ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
*
* f107, f107A, and ap effects are neither large nor well established below 80 km
* and these parameters should be set to 150., 150., and 4. respectively.
* </p>
* @param alt altitude (km)
*/
void gtd7(final double alt) {
// Calculates for thermosphere/mesosphere (above ZN2[0])
final double altt = (alt > ZN2[0]) ? alt : ZN2[0];
gts7(altt);
if (alt >= ZN2[0]) {
return;
}
// Calculates for lower mesosphere/upper stratosphere (between ZN2[0] and ZN3[0]):
// Temperature at nodes and gradients at end nodes
// Inverse temperature a linear function of spherical harmonics
final double r = PMA[2][0] * PAVGM[2];
meso_tgn2[0] = meso_tgn1[1];
meso_tn2[0] = meso_tn1[4];
meso_tn2[1] = PMA[0][0] * PAVGM[0] / (1.0 - sw[20] * glob7s(PMA[0]));
meso_tn2[2] = PMA[1][0] * PAVGM[1] / (1.0 - sw[20] * glob7s(PMA[1]));
meso_tn2[3] = PMA[2][0] * PAVGM[2] / (1.0 - sw[20] * sw[22] * glob7s(PMA[2]));
meso_tgn2[1] = PMA[9][0] * PAVGM[8] * (1.0 + sw[20] * sw[22] * glob7s(PMA[9])) *
meso_tn2[3] * meso_tn2[3] / (r * r);
meso_tn3[0] = meso_tn2[3];
// Calculates for lower stratosphere and troposphere (below ZN3[0])
// Temperature at nodes and gradients at end nodes
// Inverse temperature a linear function of spherical harmonics
if (alt < ZN3[0]) {
final double q = PMA[6][0] * PAVGM[6];
meso_tgn3[0] = meso_tgn2[1];
meso_tn3[1] = PMA[3][0] * PAVGM[3] / (1.0 - sw[22] * glob7s(PMA[3]));
meso_tn3[2] = PMA[4][0] * PAVGM[4] / (1.0 - sw[22] * glob7s(PMA[4]));
meso_tn3[3] = PMA[5][0] * PAVGM[5] / (1.0 - sw[22] * glob7s(PMA[5]));
meso_tn3[4] = PMA[6][0] * PAVGM[6] / (1.0 - sw[22] * glob7s(PMA[6]));
meso_tgn3[1] = PMA[7][0] * PAVGM[7] * (1.0 + sw[22] * glob7s(PMA[7])) *
meso_tn3[4] * meso_tn3[4] / (q * q);
}
// Linear transition to full mixing below ZN2[0]
final double dmc = (alt > ZMIX) ? 1.0 - (ZN2[0] - alt) / (ZN2[0] - ZMIX) : 0.;
final double dz28 = getDensity(MOLECULAR_NITROGEN);
// N2 density
final double dm28m = dm28 * 1.0e+06;
double dmr = dz28 / dm28m - 1.0;
double dst = densm(alt, dm28m, PDM[2][4]) * (1.0 + dmr * dmc);
setDensity(MOLECULAR_NITROGEN, dst);
// HE density
dmr = getDensity(HELIUM) / (dz28 * PDM[0][1]) - 1.0;
dst = getDensity(MOLECULAR_NITROGEN) * PDM[0][1] * (1.0 + dmr * dmc);
setDensity(HELIUM, dst);
// O density
setDensity(ATOMIC_OXYGEN, 0.);
setDensity(ANOMALOUS_OXYGEN, 0.);
// O2 density
dmr = getDensity(MOLECULAR_OXYGEN) / (dz28 * PDM[3][1]) - 1.0;
dst = getDensity(MOLECULAR_NITROGEN) * PDM[3][1] * (1.0 + dmr * dmc);
setDensity(MOLECULAR_OXYGEN, dst);
// AR density
dmr = getDensity(ARGON) / (dz28 * PDM[4][1]) - 1.0;
dst = getDensity(MOLECULAR_NITROGEN) * PDM[4][1] * (1.0 + dmr * dmc);
setDensity(ARGON, dst);
// H density
setDensity(HYDROGEN, 0.);
// N density
setDensity(ATOMIC_NITROGEN, 0.);
// Total mass density
final double tmd = AMU * (HE_MASS * getDensity(HELIUM) +
O_MASS * getDensity(ATOMIC_OXYGEN) +
N2_MASS * getDensity(MOLECULAR_NITROGEN) +
O2_MASS * getDensity(MOLECULAR_OXYGEN) +
AR_MASS * getDensity(ARGON) +
H_MASS * getDensity(HYDROGEN) +
N_MASS * getDensity(ATOMIC_NITROGEN));
setDensity(TOTAL_MASS, tmd);
// Temperature at altitude
setTemperature(ALTITUDE, densm(alt, 1.0, 0));
}
/** Calculate temperatures and densities including anomalous oxygen.
* <p></p>
* <p>NOTES ON INPUT VARIABLES:<br>
* Seconds, Local Time, and Longitude are used independently in the
* model and are not of equal importance for every situation.<br>
* For the most physically realistic calculation these three
* variables should be consistent (lst=sec/3600 + lon/15).<br>
* The Equation of Time departures from the above formula
* for apparent local time can be included if available but
* are of minor importance.<br>
* <br>
* f107 and f107A values used to generate the model correspond
* to the 10.7 cm radio flux at the actual distance of the Earth
* from the Sun rather than the radio flux at 1 AU. The following
* site provides both classes of values:<br>
* ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br>
* <br>
* f107, f107A, and ap effects are neither large nor well established below 80 km
* and these parameters should be set to 150., 150., and 4. respectively.
* </p>
* @param alt altitude (km)
*/
void gtd7d(final double alt) {
// Compute densities and temperatures
gtd7(alt);
// Update the total mass density with anomalous oxygen contribution
final double dTot = getDensity(TOTAL_MASS) + AMU * O_MASS * getDensity(ANOMALOUS_OXYGEN);
setDensity(TOTAL_MASS, dTot);
}
/** Set one density.
* @param index one of the nine elements :
* <ul>
* <li>{@link #HELIUM}</li>
* <li>{@link #ATOMIC_OXYGEN}</li>
* <li>{@link #MOLECULAR_NITROGEN}</li>
* <li>{@link #MOLECULAR_OXYGEN}</li>
* <li>{@link #ARGON}</li>
* <li>{@link #TOTAL_MASS}</li>
* <li>{@link #HYDROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* </ul>
* @param d the value of density to set
*/
void setDensity(final int index, final double d) {
densities[index] = d;
}
/** Set one temperature.
* @param index one of the two elements :
* <ul>
* <li>{@link #EXOSPHERIC}</li>
* <li>{@link #ALTITUDE}</li>
* </ul>
* @param t the value of temperature to set
*/
void setTemperature(final int index, final double t) {
temperatures[index] = t;
}
/** Get one of the stored densities.
* @param index one of the nine elements :
* <ul>
* <li>{@link #HELIUM}</li>
* <li>{@link #ATOMIC_OXYGEN}</li>
* <li>{@link #MOLECULAR_NITROGEN}</li>
* <li>{@link #MOLECULAR_OXYGEN}</li>
* <li>{@link #ARGON}</li>
* <li>{@link #TOTAL_MASS}</li>
* <li>{@link #HYDROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* </ul>
* @return the requested density
*/
public double getDensity(final int index) {
return densities[index];
}
/** Calculate G(L) function with upper thermosphere parameters.
* @param p array of parameters
* @return G(L) value
*/
private double globe7(final double[] p) {
final double[] t = new double[14];
final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
// F10.7 effect
final double df = f107 - f107a;
final double dfa = f107a - FLUX_REF;
t[0] = p[19] * df * (1.0 + p[59] * dfa) + p[20] * df * df + p[21] * dfa + p[29] * dfa * dfa;
final double f1 = 1.0 + (p[47] * dfa + p[19] * df + p[20] * df * df) * swc[1];
final double f2 = 1.0 + (p[49] * dfa + p[19] * df + p[20] * df * df) * swc[1];
// Time independent
t[1] = (p[1] * plg[0][2] + p[2] * plg[0][4] + p[22] * plg[0][6]) +
(p[14] * plg[0][2]) * dfa * swc[1] + p[26] * plg[0][1];
// Symmetrical annual
t[2] = p[18] * cd32;
// Symmetrical semiannual
t[3] = (p[15] + p[16] * plg[0][2]) * cd18;
// Asymmetrical annual
t[4] = f1 * (p[9] * plg[0][1] + p[10] * plg[0][3]) * cd14;
// Asymmetrical semiannual
t[5] = p[37] * plg[0][1] * cd39;
// Diurnal
if (sw[7] != 0) {
final double t71 = (p[11] * plg[1][2]) * cd14 * swc[5];
final double t72 = (p[12] * plg[1][2]) * cd14 * swc[5];
t[6] = f2 * ((p[3] * plg[1][1] + p[4] * plg[1][3] + p[27] * plg[1][5] + t71) * ctloc +
(p[6] * plg[1][1] + p[7] * plg[1][3] + p[28] * plg[1][5] + t72) * stloc);
}
// Semidiurnal
if (sw[8] != 0) {
final double t81 = (p[23] * plg[2][3] + p[35] * plg[2][5]) * cd14 * swc[5];
final double t82 = (p[33] * plg[2][3] + p[36] * plg[2][5]) * cd14 * swc[5];
t[7] = f2 * ((p[5] * plg[2][2] + p[41] * plg[2][4] + t81) * c2tloc +
(p[8] * plg[2][2] + p[42] * plg[2][4] + t82) * s2tloc);
}
// Terdiurnal
if (sw[14] != 0) {
t[13] = f2 * ((p[39] * plg[3][3] + (p[93] * plg[3][4] + p[46] * plg[3][6]) * cd14 * swc[5]) * s3tloc +
(p[40] * plg[3][3] + (p[94] * plg[3][4] + p[48] * plg[3][6]) * cd14 * swc[5]) * c3tloc);
}
// magnetic activity based on daily ap
if (sw[9] == -1) {
if (p[51] != 0) {
final double exp1 = FastMath.exp(-10800.0 * FastMath.abs(p[51]) /
(1.0 + p[138] * (LAT_REF - FastMath.abs(lat))));
final double p24 = FastMath.max(p[24], 1.0e-4);
apt = sg0(FastMath.min(exp1, 0.99999), p24, p[25]);
t[8] = apt * (p[50] + p[96] * plg[0][2] + p[54] * plg[0][4] +
(p[125] * plg[0][1] + p[126] * plg[0][3] + p[127] * plg[0][5]) * cd14 * swc[5] +
(p[128] * plg[1][1] + p[129] * plg[1][3] + p[130] * plg[1][5]) * swc[7] *
FastMath.cos(HOUR_TO_RAD * (hl - p[131])));
}
} else {
final double apd = ap[0] - 4.0;
final double p44 = (p[43] < 0.) ? 1.0E-5 : p[43];
final double p45 = p[44];
apdf = apd + (p45 - 1.0) * (apd + (FastMath.exp(-p44 * apd) - 1.0) / p44);
if (sw[9] != 0) {
t[8] = apdf * (p[32] + p[45] * plg[0][2] + p[34] * plg[0][4] +
(p[100] * plg[0][1] + p[101] * plg[0][3] + p[102] * plg[0][5]) * cd14 * swc[5] +
(p[121] * plg[1][1] + p[122] * plg[1][3] + p[123] * plg[1][5]) * swc[7] *
FastMath.cos(HOUR_TO_RAD * (hl - p[124])));
}
}
if (sw[10] != 0) {
final double lonr = DEG_TO_RAD * lon;
// Longitudinal
if (sw[11] != 0) {
t[10] = (1.0 + p[80] * dfa * swc[1]) *
((p[64] * plg[1][2] + p[65] * plg[1][4] + p[66] * plg[1][6] +
p[103] * plg[1][1] + p[104] * plg[1][3] + p[105] * plg[1][5] +
(p[109] * plg[1][1] + p[110] * plg[1][3] + p[111] * plg[1][5]) * swc[5] * cd14) *
FastMath.cos(lonr) +
(p[90] * plg[1][2] + p[91] * plg[1][4] + p[92] * plg[1][6] +
p[106] * plg[1][1] + p[107] * plg[1][3] + p[108] * plg[1][5] +
(p[112] * plg[1][1] + p[113] * plg[1][3] + p[114] * plg[1][5]) * swc[5] * cd14) *
FastMath.sin(lonr));
}
// ut and mixed ut, longitude
if (sw[12] != 0) {
t[11] = (1.0 + p[95] * plg[0][1]) * (1.0 + p[81] * dfa * swc[1]) *
(1.0 + p[119] * plg[0][1] * swc[5] * cd14) *
(p[68] * plg[0][1] + p[69] * plg[0][3] + p[70] * plg[0][5]) *
FastMath.cos(SEC_TO_RAD * (sec - p[71]));
t[11] += swc[11] * (1.0 + p[137] * dfa * swc[1]) *
(p[76] * plg[2][3] + p[77] * plg[2][5] + p[78] * plg[2][7]) *
FastMath.cos(SEC_TO_RAD * (sec - p[79]) + 2.0 * lonr);
}
/* ut, longitude magnetic activity */
if (sw[13] != 0) {
if (sw[9] == -1) {
if (p[51] != 0.) {
t[12] = apt * swc[11] * (1. + p[132] * plg[0][1]) *
(p[52] * plg[1][2] + p[98] * plg[1][4] + p[67] * plg[1][6]) *
FastMath.cos(DEG_TO_RAD * (lon - p[97])) +
apt * swc[11] * swc[5] * cd14 *
(p[133] * plg[1][1] + p[134] * plg[1][3] + p[135] * plg[1][5]) *
FastMath.cos(DEG_TO_RAD * (lon - p[136])) +
apt * swc[12] *
(p[55] * plg[0][1] + p[56] * plg[0][3] + p[57] * plg[0][5]) *
FastMath.cos(SEC_TO_RAD * (sec - p[58]));
}
} else {
t[12] = apdf * swc[11] * (1.0 + p[120] * plg[0][1]) *
((p[60] * plg[1][2] + p[61] * plg[1][4] + p[62] * plg[1][6]) *
FastMath.cos(DEG_TO_RAD * (lon - p[63]))) +
apdf * swc[11] * swc[5] * cd14 *
(p[115] * plg[1][1] + p[116] * plg[1][3] + p[117] * plg[1][5]) *
FastMath.cos(DEG_TO_RAD * (lon - p[118])) +
apdf * swc[12] *
(p[83] * plg[0][1] + p[84] * plg[0][3] + p[85] * plg[0][5]) *
FastMath.cos(SEC_TO_RAD * (sec - p[75]));
}
}
}
// Sum all effects (params not used: 82, 89, 99, 139-149)
double tinf = p[30];
for (int i = 0; i < 14; i++) {
tinf += FastMath.abs(sw[i + 1]) * t[i];
}
// Return G(L)
return tinf;
}
/** Calculate G(L) function with lower atmosphere parameters.
* @param p array of parameters
* @return G(L) value
*/
private double glob7s(final double[] p) {
final double[] t = new double[14];
final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
// F10.7 effect
t[0] = p[21] * (f107a - FLUX_REF);
// Time independent
t[1] = p[1] * plg[0][2] + p[2] * plg[0][4] + p[22] * plg[0][6] +
p[26] * plg[0][1] + p[14] * plg[0][3] + p[59] * plg[0][5];
// Symmetrical annual
t[2] = (p[18] + p[47] * plg[0][2] + p[29] * plg[0][4]) * cd32;
// Symmetrical semiannual
t[3] = (p[15] + p[16] * plg[0][2] + p[30] * plg[0][4]) * cd18;
// Asymmetrical annual
t[4] = (p[9] * plg[0][1] + p[10] * plg[0][3] + p[20] * plg[0][5]) * cd14;
// Asymmetrical semiannual
t[5] = (p[37] * plg[0][1]) * cd39;
// Diurnal
if (sw[7] != 0) {
final double t71 = p[11] * plg[1][2] * cd14 * swc[5];
final double t72 = p[12] * plg[1][2] * cd14 * swc[5];
t[6] = (p[3] * plg[1][1] + p[4] * plg[1][3] + t71) * ctloc +
(p[6] * plg[1][1] + p[7] * plg[1][3] + t72) * stloc;
}
// Semidiurnal
if (sw[8] != 0) {
final double t81 = (p[23] * plg[2][3] + p[35] * plg[2][5]) * cd14 * swc[5];
final double t82 = (p[33] * plg[2][3] + p[36] * plg[2][5]) * cd14 * swc[5];
t[7] = (p[5] * plg[2][2] + p[41] * plg[2][4] + t81) * c2tloc +
(p[8] * plg[2][2] + p[42] * plg[2][4] + t82) * s2tloc;
}
// Terdiurnal
if (sw[14] != 0) {
t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
}
// Magnetic activity
if (sw[9] == 1) {
t[8] = apdf * (p[32] + p[45] * plg[0][2] * swc[2]);
} else if (sw[9] == -1) {
t[8] = apt * (p[50] + p[96] * plg[0][2] * swc[2]);
}
// Longitudinal
if (!(sw[10] == 0 || sw[11] == 0)) {
final double lonr = DEG_TO_RAD * lon;
t[10] = (1.0 + plg[0][1] * (p[80] * swc[5] * FastMath.cos(DAY_TO_RAD * (doy - p[81])) +
p[85] * swc[6] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[86]))) +
p[83] * swc[3] * FastMath.cos(DAY_TO_RAD * (doy - p[84])) +
p[87] * swc[4] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[88]))) *
((p[64] * plg[1][2] + p[65] * plg[1][4] + p[66] * plg[1][6] +
p[74] * plg[1][1] + p[75] * plg[1][3] + p[76] * plg[1][5]) * FastMath.cos(lonr) +
(p[90] * plg[1][2] + p[91] * plg[1][4] + p[92] * plg[1][6] +
p[77] * plg[1][1] + p[78] * plg[1][3] + p[79] * plg[1][5]) * FastMath.sin(lonr));
}
// Sum all effects
double gl = 0;
for (int i = 0; i < 14; i++) {
gl += FastMath.abs(sw[i + 1]) * t[i];
}
// Return G(L)
return gl;
}
/** Implements sg0 function (Eq. A24a).
* @param ex ex
* @param p24 abs(p[24])
* @param p25 p[25]
* @return sg0
*/
private double sg0(final double ex, final double p24, final double p25) {
final double g01 = g0(ap[1], p24, p25);
final double g02 = g0(ap[2], p24, p25);
final double g03 = g0(ap[3], p24, p25);
final double g04 = g0(ap[4], p24, p25);
final double g05 = g0(ap[5], p24, p25);
final double g06 = g0(ap[6], p24, p25);
final double ex2 = ex * ex;
final double ex3 = ex * ex2;
final double ex4 = ex2 * ex2;
final double ex8 = ex4 * ex4;
final double ex12 = ex4 * ex8;
final double g234 = g02 * ex + g03 * ex2 + g04 * ex3;
final double g56 = g05 * ex4 + g06 * ex12;
final double ex19 = ex3 * ex4 * ex12;
final double omex = 1.0 - ex;
final double sumex = 1.0 + (1.0 - ex19) / omex * FastMath.sqrt(ex);
return (g01 + (g234 + g56 * (1.0 - ex8) / omex)) / sumex;
}
/** Implements go function (Eq. A24d).
* @param apI 3 hrs ap
* @param p24 abs(p[24])
* @param p25 p[25]
* @return go
*/
private double g0(final double apI, final double p24, final double p25) {
final double am4 = apI - 4.0;
return am4 + (p25 - 1.0) * (am4 + (FastMath.exp(-p24 * am4) - 1.0) / p24);
}
/** Calculates chemistry/dissociation correction for MSIS models.
* @param alt altitude
* @param r target ratio
* @param h1 transition scale length
* @param zh altitude of 1/2 R
* @return correction
*/
private double ccor(final double alt, final double r, final double h1, final double zh) {
final double e = (alt - zh) / h1;
if (e > 70.) {
return 1.;
} else if (e < -70.) {
return FastMath.exp(r);
} else {
return FastMath.exp(r / (1.0 + FastMath.exp(e)));
}
}
/** Calculates O & O2 chemistry/dissociation correction for MSIS models.
* @param alt altitude
* @param r target ratio
* @param h1 transition scale length
* @param zh altitude of 1/2 R
* @param h2 transition scale length
* @return correction
*/
private double ccor2(final double alt, final double r,
final double h1, final double zh, final double h2) {
final double e1 = (alt - zh) / h1;
final double e2 = (alt - zh) / h2;
if ((e1 > 70.) || (e2 > 70.)) {
return 1.;
} else if ((e1 < -70.) && (e2 < -70.)) {
return FastMath.exp(r);
} else {
final double ex1 = FastMath.exp(e1);
final double ex2 = FastMath.exp(e2);
return FastMath.exp(r / (1.0 + 0.5 * (ex1 + ex2)));
}
}
/** Calculates scale height.
* @param alt altitude
* @param xm species molecular weight
* @param temp temperature
* @return scale height (km)
*/
private double scalh(final double alt, final double xm, final double temp) {
// Gravity at altitude
final double denom = 1.0 + alt / rlat;
final double galt = glat / (denom * denom);
return R_GAS * temp / (galt * xm);
}
/** Calculates turbopause correction for MSIS models.
* @param dd diffusive density
* @param dm full mixed density
* @param zhm transition scale length
* @param xmm full mixed molecular weight
* @param xm species molecular weight
* @return combined density
*/
private double dnet(final double dd, final double dm,
final double zhm, final double xmm, final double xm) {
if (!(dm > 0 && dd > 0)) {
double ddd = dd;
if (dd == 0 && dm == 0) {
ddd = 1;
}
if (dm == 0) {
return ddd;
}
if (dd == 0) {
return dm;
}
}
final double a = zhm / (xmm - xm);
final double ylog = a * FastMath.log(dm / dd);
if (ylog < -10.) {
return dd;
} else if (ylog > 10.) {
return dm;
} else {
return dd * FastMath.pow(1.0 + FastMath.exp(ylog), 1.0 / a);
}
}
/** Integrate cubic spline function from xa[0] to x.
* <p>ADAPTED FROM NUMERICAL RECIPES</p>
* @param xa array of abscissas in ascending order
* @param ya array of ordinates in ascending order by xa
* @param y2a array of second derivatives in ascending order by xa
* @param x abscissa end point
* @return integral value
*/
private double splini(final double[] xa, final double[] ya, final double[] y2a, final double x) {
final int n = xa.length;
double yi = 0;
int klo = 0;
int khi = 1;
while (x > xa[klo] && khi < n) {
double xx = x;
if (khi < n - 1) {
xx = (x < xa[khi]) ? x : xa[khi];
}
final double h = xa[khi] - xa[klo];
final double a = (xa[khi] - xx) / h;
final double b = (xx - xa[klo]) / h;
final double a2 = a * a;
final double b2 = b * b;
yi += ((1.0 - a2) * ya[klo] / 2.0 + b2 * ya[khi] / 2.0 +
((-(1.0 + a2 * a2) / 4.0 + a2 / 2.0) * y2a[klo] +
(b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi]) * h * h / 6.0) * h;
klo++;
khi++;
}
return yi;
}
/** Calculate cubic spline interpolated value.
* <p>ADAPTED FROM NUMERICAL RECIPES</p>
* @param xa array of abscissas in ascending order
* @param ya array of ordinates in ascending order by xa
* @param y2a array of second derivatives in ascending order by xa
* @param x abscissa for interpolation
* @return interpolated value
*/
private double splint(final double[] xa, final double[] ya, final double[] y2a, final double x) {
final int n = xa.length;
int klo = 0;
int khi = n - 1;
while (khi - klo > 1) {
final int k = (khi + klo) >>> 1;
if (xa[k] > x) {
khi = k;
} else {
klo = k;
}
}
final double h = xa[khi] - xa[klo];
final double a = (xa[khi] - x) / h;
final double b = (x - xa[klo]) / h;
return a * ya[klo] + b * ya[khi] +
((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * h * h / 6.0;
}
/** Calculate 2nd derivatives of cubic spline interpolation function.
* <p>ADAPTED FROM NUMERICAL RECIPES</p>
* @param x array of abscissas in ascending order
* @param y array of ordinates in ascending order by x
* @param yp1 derivative at x[0] (2nd derivatives null if > 1E30)
* @param ypn derivative at x[n-1] (2nd derivatives null if > 1E30)
* @return array of second derivatives
*/
private double[] spline(final double[] x, final double[] y, final double yp1, final double ypn) {
final int n = x.length;
final double[] y2 = new double[n];
final double[] u = new double[n];
if (yp1 < 1e+30) {
y2[0] = -0.5;
u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
}
for (int i = 1; i < n - 1; i++) {
final double sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
final double p = sig * y2[i - 1] + 2.0;
y2[i] = (sig - 1.0) / p;
u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) /
(x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
}
double qn = 0;
double un = 0;
if (ypn < 1e+30) {
qn = 0.5;
un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
}
y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
for (int k = n - 2; k >= 0; k--) {
y2[k] = y2[k] * y2[k + 1] + u[k];
}
return y2;
}
/** Calculate Temperature and Density Profiles for lower atmosphere.
* @param alt altitude
* @param d0 density
* @param xm mixed density
* @return temperature or density profile
*/
private double densm(final double alt, final double d0, final double xm) {
double densm = d0;
// stratosphere/mesosphere temperature
int mn = ZN2.length;
double z = (alt > ZN2[mn - 1]) ? alt : ZN2[mn - 1];
double z1 = ZN2[0];
double z2 = ZN2[mn - 1];
double t1 = meso_tn2[0];
double t2 = meso_tn2[mn - 1];
double zg = zeta(z, z1);
double zgdif = zeta(z2, z1);
/* set up spline nodes */
double[] xs = new double[mn];
double[] ys = new double[mn];
for (int k = 0; k < mn; k++) {
xs[k] = zeta(ZN2[k], z1) / zgdif;
ys[k] = 1.0 / meso_tn2[k];
}
final double qSM = (rlat + z2) / (rlat + z1);
double yd1 = -meso_tgn2[0] / (t1 * t1) * zgdif;
double yd2 = -meso_tgn2[1] / (t2 * t2) * zgdif * qSM * qSM;
/* calculate spline coefficients */
double[] y2out = spline(xs, ys, yd1, yd2);
double x = zg / zgdif;
double y = splint(xs, ys, y2out, x);
/* temperature at altitude */
double tz = 1.0 / y;
if (xm != 0.0) {
/* calculate stratosphere / mesospehere density */
final double glb = galt(z1);
final double gamm = xm * glb * zgdif / R_GAS;
/* Integrate temperature profile */
final double yi = splini(xs, ys, y2out, x);
final double expl = FastMath.min(50., gamm * yi);
/* Density at altitude */
densm *= (t1 / tz) * FastMath.exp(-expl);
}
if (alt > ZN3[0]) {
return (xm == 0.0) ? tz : densm;
}
// troposhere/stratosphere temperature
z = alt;
mn = ZN3.length;
z1 = ZN3[0];
z2 = ZN3[mn - 1];
t1 = meso_tn3[0];
t2 = meso_tn3[mn - 1];
zg = zeta(z, z1);
zgdif = zeta(z2, z1);
/* set up spline nodes */
xs = new double[mn];
ys = new double[mn];
for (int k = 0; k < mn; k++) {
xs[k] = zeta(ZN3[k], z1) / zgdif;
ys[k] = 1.0 / meso_tn3[k];
}
final double qTS = (rlat + z2) / (rlat + z1);
yd1 = -meso_tgn3[0] / (t1 * t1) * zgdif;
yd2 = -meso_tgn3[1] / (t2 * t2) * zgdif * qTS * qTS;
/* calculate spline coefficients */
y2out = spline(xs, ys, yd1, yd2);
x = zg / zgdif;
y = splint(xs, ys, y2out, x);
/* temperature at altitude */
tz = 1.0 / y;
if (xm != 0.0) {
/* calculate tropospheric / stratosphere density */
final double glb = galt(z1);
final double gamm = xm * glb * zgdif / R_GAS;
/* Integrate temperature profile */
final double yi = splini(xs, ys, y2out, x);
final double expl = FastMath.min(50., gamm * yi);
/* Density at altitude */
densm *= (t1 / tz) * FastMath.exp(-expl);
}
return (xm == 0.0) ? tz : densm;
}
/** Calculate temperature and density profiles according to new lower thermo polynomial.
* @param alt altitude
* @param dlb density at lower boundary
* @param tinf exospheric temperature
* @param tlb temperature at lower boundary
* @param xm species molecular weight
* @param alpha thermal diffusion coefficient
* @param zlb altitude of the lower boundary
* @param s2 slope
* @return temperature or density profile
*/
private double densu(final double alt, final double dlb, final double tinf,
final double tlb, final double xm, final double alpha,
final double zlb, final double s2) {
/* joining altitudes of Bates and spline */
double z = (alt > ZN1[0]) ? alt : ZN1[0];
/* geopotential altitude difference from ZLB */
final double zg2 = zeta(z, zlb);
/* Bates temperature */
final double tt = tinf - (tinf - tlb) * FastMath.exp(-s2 * zg2);
final double ta = tt;
double tz = tt;
final int mn = ZN1.length;
final double[] xs = new double[mn];
final double[] ys = new double[mn];
double x = 0.;
double[] y2out = new double[mn];
double zgdif = 0.;
if (alt < ZN1[0]) {
/* calculate temperature below ZA
* temperature gradient at ZA from Bates profile */
final double p = (rlat + zlb) / (rlat + ZN1[0]);
final double dta = (tinf - ta) * s2 * p * p;
meso_tgn1[0] = dta;
meso_tn1[0] = ta;
z = (alt > ZN1[mn - 1]) ? alt : ZN1[mn - 1];
final double t1 = meso_tn1[0];
final double t2 = meso_tn1[mn - 1];
/* geopotental difference from z1 */
final double zg = zeta(z, ZN1[0]);
zgdif = zeta(ZN1[mn - 1], ZN1[0]);
/* set up spline nodes */
for (int k = 0; k < mn; k++) {
xs[k] = zeta(ZN1[k], ZN1[0]) / zgdif;
ys[k] = 1.0 / meso_tn1[k];
}
/* end node derivatives */
final double q = (rlat + ZN1[mn - 1]) / (rlat + ZN1[0]);
final double yd1 = -meso_tgn1[0] / (t1 * t1) * zgdif;
final double yd2 = -meso_tgn1[1] / (t2 * t2) * zgdif * q * q;
/* calculate spline coefficients */
y2out = spline(xs, ys, yd1, yd2);
x = zg / zgdif;
final double y = splint(xs, ys, y2out, x);
/* temperature at altitude */
tz = 1.0 / y;
}
if (xm == 0) {
return tz;
}
/* calculate density above za */
double glb = galt(zlb);
double gamma = xm * glb / (R_GAS * s2 * tinf);
double expl = (tt <= 0) ? 50. : FastMath.min(50., FastMath.exp(-s2 * gamma * zg2));
double densu = dlb * FastMath.pow(tlb / tt, 1.0 + alpha + gamma) * expl;
/* calculate density below za */
if (alt < ZN1[0]) {
glb = galt(ZN1[0]);
gamma = xm * glb * zgdif / R_GAS;
/* integrate spline temperatures */
expl = (tz <= 0) ? 50.0 : FastMath.min(50., gamma * splini(xs, ys, y2out, x));
/* correct density at altitude */
densu *= FastMath.pow(meso_tn1[0] / tz, 1.0 + alpha) * FastMath.exp(-expl);
}
/* Return density at altitude */
return densu;
}
/** Calculate gravity at altitude.
* @param alt altitude (km)
* @return gravity at altitude (cm/s2)
*/
private double galt(final double alt) {
final double r = 1.0 + alt / rlat;
return glat / (r * r);
}
/** Calculate zeta function.
* @param zz zz value
* @param zl zl value
* @return value of zeta function
*/
private double zeta(final double zz, final double zl) {
return (zz - zl) * (rlat + zl) / (rlat + zz);
}
}
/**
* This class is a placeholder for the computed densities and temperatures.
* <p>
* Densities are provided as an array d such as:
* <ul>
* <li>d[0] = He number density (1/m³)</li>
* <li>d[1] = O number density (1/m³)</li>
* <li>d[2] = N2 number density (1/m³)</li>
* <li>d[3] = O2 number density (1/m³)</li>
* <li>d[4] = Ar number density (1/m³)</li>
* <li>d[5] = total mass density (kg/m³) (*)</li>
* <li>d[6] = H number density (1/m³)</li>
* <li>d[7] = N number density (1/m³)</li>
* <li>d[8] = anomalous oxygen number density (1/m³)
* </ul>
* Total mass density, d[5], is NOT the same for methods gtd7 and gtd7d:
* <ul>
* <li>For gtd7: d[5] is the sum of the mass densities of the species
* He, O, N2, O2, Ar, H and N but does NOT include anomalous oxygen.</li>
* <li>For gtd7d: d[5] is the "effective total mass density for drag" and is the sum
* of the mass densities of all species in this model, INCLUDING anomalous oxygen.</li>
* </ul>
* O, H, and N are set to zero below 72.5 km.
* </p>
* <p>
* Temperatures are provided as an array t such as:
* <ul>
* <li>t[0] = exospheric temperature (K)</li>
* <li>t[1] = temperature at altitude (K)</li>
* </ul>
* t[0] is set to global average for altitudes below 120 km.<br>
* The 120 km gradient is left at global average value for altitudes below 72 km.
* </p>
* @param <T> type of the field elements
* @since 9.0
*/
public class FieldOutput<T extends RealFieldElement<T>> {
/** Type of the field elements. */
private final Field<T> field;
/** Zero for the field. */
private final T zero;
/** Day of year (from 1 to 365 or 366). */
private final int doy;
/** Seconds in day (UT scale). */
private final T sec;
/** Geodetic latitude (°). */
private final T lat;
/** Geodetic longitude (°). */
private final T lon;
/** Local apparent solar time (hours). */
private final T hl;
/** 81 day average of F10.7 flux (centered on day). */
private final double f107a;
/** Daily F10.7 flux for previous day. */
private final double f107;
/** Array containing:
* <ul>
* <li>0: daily Ap</li>
* <li>1: 3 hr ap index for current time</li>
* <li>2: 3 hr ap index for 3 hrs before current time</li>
* <li>3: 3 hr ap index for 6 hrs before current time</li>
* <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
* <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
* <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
* </ul>. */
private final double[] ap;
/** Gravity at latitude (cm/s2). */
private final T glat;
/** Effective Earth radius at latitude (km). */
private final T rlat;
/** N2 mixed density at alt. */
private T dm28;
/** Legendre polynomials. */
private final T[][] plg;
/** Cosinus of local solar time. */
private final T ctloc;
/** Sinus of local solar time. */
private final T stloc;
/** Square of ctloc. */
private final T c2tloc;
/** Square of stloc. */
private final T s2tloc;
/** Cube of ctloc. */
private final T c3tloc;
/** Cube of stloc. */
private final T s3tloc;
/** Magnetic activity based on daily ap. */
private double apdf;
/** Magnetic activity based on daily ap. */
private T apt;
/** Temperature at nodes for ZN1 scale. */
private final T[] meso_tn1;
/** Temperature at nodes for ZN2 scale. */
private final T[] meso_tn2;
/** Temperature at nodes for ZN3 scale. */
private final T[] meso_tn3;
/** Temperature gradients at end nodes for ZN1 scale. */
private final T[] meso_tgn1;
/** Temperature gradients at end nodes for ZN2 scale. */
private final T[] meso_tgn2;
/** Temperature gradients at end nodes for ZN3 scale. */
private final T[] meso_tgn3;
/** Densities. */
private final T[] densities;
/** Temperatures. */
private final T[] temperatures;
/** Simple constructor.
* @param doy day of year (from 1 to 365 or 366)
* @param sec seconds in day (UT scale)
* @param lat geodetic latitude (°)
* @param lon geodetic longitude (°)
* @param hl local apparent solar time (hours)
* @param f107a 81 day average of F10.7 flux (centered on day)
* @param f107 daily F10.7 flux for previous day
* @param ap array containing:
* <ul>
* <li>0: daily Ap</li>
* <li>1: 3 hr ap index for current time</li>
* <li>2: 3 hr ap index for 3 hrs before current time</li>
* <li>3: 3 hr ap index for 6 hrs before current time</li>
* <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
* <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
* <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
* </ul>
*/
FieldOutput(final int doy, final T sec,
final T lat, final T lon, final T hl,
final double f107a, final double f107, final double[] ap) {
this.field = sec.getField();
this.zero = field.getZero();
this.doy = doy;
this.sec = sec;
this.lat = lat;
this.lon = lon;
this.hl = hl;
this.f107a = f107a;
this.f107 = f107;
this.ap = ap.clone();
this.plg = MathArrays.buildArray(field, 4, 8);
this.meso_tn1 = MathArrays.buildArray(field, ZN1.length);
this.meso_tn2 = MathArrays.buildArray(field, ZN2.length);
this.meso_tn3 = MathArrays.buildArray(field, ZN3.length);
this.meso_tgn1 = MathArrays.buildArray(field, 2);
this.meso_tgn2 = MathArrays.buildArray(field, 2);
this.meso_tgn3 = MathArrays.buildArray(field, 2);
densities = MathArrays.buildArray(field, 9);
temperatures = MathArrays.buildArray(field, 2);
// Calculates latitude variable gravity and effective radius
final T xlat = (sw[2] == 0) ? zero.add(LAT_REF) : lat;
final T c2 = xlat.multiply(2 * DEG_TO_RAD).cos();
glat = c2.multiply(-0.0026373).add(1).multiply(G_REF);
rlat = glat.multiply(2).divide(c2.multiply(2.27e-9).add(3.085462e-6)).multiply(1.e-5);
// Convert latitude into radians
final T latr = lat.multiply(DEG_TO_RAD);
// Calculate legendre polynomials
final T c = latr.sin();
final T s = latr.cos();
plg[0][1] = c;
plg[0][2] = c.multiply( 3.0).multiply(plg[0][1]).subtract(1.0).divide(2.0);
plg[0][3] = c.multiply( 5.0).multiply(plg[0][2]).subtract(plg[0][1].multiply(2.0)).divide(3.0);
plg[0][4] = c.multiply( 7.0).multiply(plg[0][3]).subtract(plg[0][2].multiply(3.0)).divide(4.0);
plg[0][5] = c.multiply( 9.0).multiply(plg[0][4]).subtract(plg[0][3].multiply(4.0)).divide(5.0);
plg[0][6] = c.multiply(11.0).multiply(plg[0][5]).subtract(plg[0][4].multiply(5.0)).divide(6.0);
plg[1][1] = s;
plg[1][2] = c.multiply( 3.0).multiply(plg[1][1]);
plg[1][3] = c.multiply( 5.0).multiply(plg[1][2]).subtract(plg[1][1].multiply(3.0)).divide(2.0);
plg[1][4] = c.multiply( 7.0).multiply(plg[1][3]).subtract(plg[1][2].multiply(4.0)).divide(3.0);
plg[1][5] = c.multiply( 9.0).multiply(plg[1][4]).subtract(plg[1][3].multiply(5.0)).divide(4.0);
plg[1][6] = c.multiply(11.0).multiply(plg[1][5]).subtract(plg[1][4].multiply(6.0)).divide(5.0);
plg[2][2] = s.multiply( 3.0).multiply(plg[1][1]);
plg[2][3] = c.multiply( 5.0).multiply(plg[2][2]);
plg[2][4] = c.multiply( 7.0).multiply(plg[2][3]).subtract(plg[2][2].multiply(5.0)).divide(2.0);
plg[2][5] = c.multiply( 9.0).multiply(plg[2][4]).subtract(plg[2][3].multiply(6.0)).divide(3.0);
plg[2][6] = c.multiply(11.0).multiply(plg[2][5]).subtract(plg[2][4].multiply(7.0)).divide(4.0);
plg[2][7] = c.multiply(13.0).multiply(plg[2][6]).subtract(plg[2][5].multiply(8.0)).divide(5.0);
plg[3][3] = s.multiply( 5.0).multiply(plg[2][2]);
plg[3][4] = c.multiply( 7.0).multiply(plg[3][3]);
plg[3][5] = c.multiply( 9.0).multiply(plg[3][4]).subtract(plg[3][3].multiply(7.0)).divide(2.0);
plg[3][6] = c.multiply(11.0).multiply(plg[3][5]).subtract(plg[3][4].multiply(8.0)).divide(3.0);
// Calculate additional data
if (!(sw[7] == 0 && sw[8] == 0 && sw[14] == 0)) {
final T tloc = hl.multiply(HOUR_TO_RAD);
final T tlx2 = tloc.add(tloc);
final T tlx3 = tloc.add(tlx2);
stloc = tloc.sin();
ctloc = tloc.cos();
s2tloc = tlx2.sin();
c2tloc = tlx2.cos();
s3tloc = tlx3.sin();
c3tloc = tlx3.cos();
} else {
stloc = zero;
ctloc = zero;
s2tloc = zero;
c2tloc = zero;
s3tloc = zero;
c3tloc = zero;
}
}
/** Calculate temperatures and densities not including anomalous oxygen.
* <p>
* This method is the thermospheric portion of NRLMSISE-00 for alt > 72.5 km.
* </p>
* <p>NOTES ON INPUT VARIABLES:<br>
* Seconds, Local Time, and Longitude are used independently in the
* model and are not of equal importance for every situation.<br>
* For the most physically realistic calculation these three
* variables should be consistent (lst=sec/3600 + lon/15).<br>
* The Equation of Time departures from the above formula
* for apparent local time can be included if available but
* are of minor importance.<br><br>
*
* f107 and f107A values used to generate the model correspond
* to the 10.7 cm radio flux at the actual distance of the Earth
* from the Sun rather than the radio flux at 1 AU. The following
* site provides both classes of values:<br>
* ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
*
* f107, f107A, and ap effects are neither large nor well established below 80 km
* and these parameters should be set to 150., 150., and 4. respectively.
* </p>
* @param alt altitude (km)
*/
void gts7(final T alt) {
// Thermal diffusion coefficients for species
final double[] alpha = {-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
// Altitude limits for net density computation for species
final double[] altl = {200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
// N2 mixed density
final double xmm = PDM[2][4];
/**** Exospheric temperature ****/
T tinf = zero.add(PTM[0] * PT[0]);
// Tinf variations not important below ZA or ZN[0]
if (alt.getReal() > ZN1[0]) {
tinf = tinf.multiply(globe7(PT).multiply(sw[16]).add(1));
}
setTemperature(EXOSPHERIC, tinf);
// Gradient variations not important below ZN[4]
T g0 = zero.add(PTM[3] * PS[0]);
if (alt.getReal() > ZN1[4]) {
g0 = g0.multiply(globe7(PS).multiply(sw[19]).add(1));
}
// Temperature at lower boundary
T tlb = zero.add(PTM[1] * PD[3][0]);
tlb = tlb.multiply(globe7(PD[3]).multiply(sw[17]).add(1));
// Slope
final T s = g0.divide(tinf.subtract(tlb));
// Lower thermosphere temp variations not significant for density above 300 km
meso_tn1[1] = zero.add(PTM[6] * PTL[0][0]);
meso_tn1[2] = zero.add(PTM[2] * PTL[1][0]);
meso_tn1[3] = zero.add(PTM[7] * PTL[2][0]);
meso_tn1[4] = zero.add(PTM[4] * PTL[3][0]);
meso_tgn1[1] = zero.add(PTM[8] * PMA[8][0]);
if (alt.getReal() < 300.0) {
final double r = PTM[4] * PTL[3][0];
meso_tn1[1] = meso_tn1[1].divide(glob7s(PTL[0]).multiply(sw[18] ).negate().add(1));
meso_tn1[2] = meso_tn1[2].divide(glob7s(PTL[1]).multiply(sw[18] ).negate().add(1));
meso_tn1[3] = meso_tn1[3].divide(glob7s(PTL[2]).multiply(sw[18] ).negate().add(1));
meso_tn1[4] = meso_tn1[4].divide(glob7s(PTL[3]).multiply(sw[18] * sw[20]).negate().add(1));
meso_tgn1[1] = meso_tgn1[1].multiply(glob7s(PMA[8]).multiply(sw[18] * sw[20]).add(1));
meso_tgn1[1] = meso_tgn1[1].multiply(meso_tn1[4].multiply(meso_tn1[4]).divide(r * r));
}
/**** Temperature at altitude ****/
setTemperature(ALTITUDE, densu(alt, zero.add(1.0), tinf, tlb, 0, 0, PTM[5], s));
/**** N2 density ****/
/* Density variation factor at Zlb */
final T g28 = globe7(PD[2]).multiply(sw[21]);
/* Diffusive density at Zlb */
final T db28 = g28.exp().multiply(PDM[2][0] * PD[2][0]);
/* Diffusive density at Alt */
T diffusiveDensity = densu(alt, db28, tinf, tlb, N2_MASS, alpha[2], PTM[5], s);
setDensity(MOLECULAR_NITROGEN, diffusiveDensity);
// Variation of turbopause height
final T zhf = lat.multiply(DEG_TO_RAD).sin().
multiply(sw[5] * PDL[0][24] * FastMath.cos(DAY_TO_RAD * (doy - PT[13]))).
add(1).
multiply(PDL[1][24]);
/* Turbopause */
final T zh28 = zhf.multiply(PDM[2][2]);
final double zhm28 = PDM[2][3] * PDL[1][5];
/* Mixed density at Zlb */
final T b28 = densu(zh28, db28, tinf, tlb, N2_MASS - xmm, alpha[2] - 1.0, PTM[5], s);
if (sw[15] != 0 && alt.getReal() <= altl[2]) {
/* Mixed density at Alt */
dm28 = densu(alt, b28, tinf, tlb, xmm, alpha[2], PTM[5], s);
/* Net density at Alt */
setDensity(MOLECULAR_NITROGEN, dnet(diffusiveDensity, dm28, zhm28, xmm, N2_MASS));
} else {
dm28 = zero;
}
/**** He density ****/
/* Density variation factor at Zlb */
final T g4 = globe7(PD[0]).multiply(sw[21]);
/* Diffusive density at Zlb */
final T db04 = g4.exp().multiply(PDM[0][0] * PD[0][0]);
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db04, tinf, tlb, HE_MASS, alpha[0], PTM[5], s);
setDensity(HELIUM, diffusiveDensity);
if (sw[15] != 0 && alt.getReal() < altl[0]) {
/* Turbopause */
final double zh04 = PDM[0][2];
/* Mixed density at Zlb */
final T b04 = densu(zero.add(zh04), db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm04 = densu(alt, b04, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm04 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm04, zhm04, xmm, HE_MASS);
/* Correction to specified mixing ratio at ground */
final T rl = b28.multiply(PDM[0][1]).divide(b04).log();
final double zc04 = PDM[0][4] * PDL[1][0];
final double hc04 = PDM[0][5] * PDL[1][1];
/* Net density corrected at Alt */
setDensity(HELIUM, diffusiveDensity.multiply(ccor(alt, rl, hc04, zc04)));
}
/**** O density ****/
/* Density variation factor at Zlb */
final T g16 = globe7(PD[1]).multiply(sw[21]);
/* Diffusive density at Zlb */
final T db16 = g16.exp().multiply(PDM[1][0] * PD[1][0]);
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db16, tinf, tlb, O_MASS, alpha[1], PTM[5], s);
setDensity(ATOMIC_OXYGEN, diffusiveDensity);
if (sw[15] != 0 && alt.getReal() < altl[1]) {
/* Turbopause */
final double zh16 = PDM[1][2];
/* Mixed density at Zlb */
final T b16 = densu(zero.add(zh16), db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
/* Mixed density at Alt */
final T dm16 = densu(alt, b16, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm16 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm16, zhm16, xmm, O_MASS);
final double rl = PDM[1][1] * PDL[1][16] * (1.0 + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
final double hc16 = PDM[1][5] * PDL[1][3];
final double zc16 = PDM[1][4] * PDL[1][2];
final double hc216 = PDM[1][5] * PDL[1][4];
diffusiveDensity = diffusiveDensity.multiply(ccor2(alt, rl, hc16, zc16, hc216));
/* Chemistry correction */
final double hcc16 = PDM[1][7] * PDL[1][13];
final double zcc16 = PDM[1][6] * PDL[1][12];
final double rc16 = PDM[1][3] * PDL[1][14];
/* Net density corrected at Alt */
setDensity(ATOMIC_OXYGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc16), hcc16, zcc16)));
}
/**** O2 density ****/
/* Density variation factor at Zlb */
final T g32 = globe7(PD[4]).multiply(sw[21]);
/* Diffusive density at Zlb */
final T db32 = g32.exp().multiply(PDM[3][0] * PD[4][0]);
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db32, tinf, tlb, O2_MASS, alpha[3], PTM[5], s);
setDensity(MOLECULAR_OXYGEN, diffusiveDensity);
if (sw[15] != 0) {
if (alt.getReal() <= altl[3]) {
/* Turbopause */
final double zh32 = PDM[3][2];
/* Mixed density at Zlb */
final T b32 = densu(zero.add(zh32), db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm32 = densu(alt, b32, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm32 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm32, zhm32, xmm, O2_MASS);
/* Correction to specified mixing ratio at ground */
final T rl = b28.multiply(PDM[3][1]).divide(b32).log();
final double hc32 = PDM[3][5] * PDL[1][7];
final double zc32 = PDM[3][4] * PDL[1][6];
diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc32, zc32));
}
/* Correction for general departure from diffusive equilibrium above Zlb */
final double hcc32 = PDM[3][7] * PDL[1][22];
final double hcc232 = PDM[3][7] * PDL[0][22];
final double zcc32 = PDM[3][6] * PDL[1][21];
final double rc32 = PDM[3][3] * PDL[1][23] * (1. + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
/* Net density corrected at Alt */
setDensity(MOLECULAR_OXYGEN, diffusiveDensity.multiply(ccor2(alt, rc32, hcc32, zcc32, hcc232)));
}
/**** Ar density ****/
/* Density variation factor at Zlb */
final T g40 = globe7(PD[5]).multiply(sw[21]);
/* Diffusive density at Zlb */
final T db40 = g40.exp().multiply(PDM[4][0] * PD[5][0]);
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db40, tinf, tlb, AR_MASS, alpha[4], PTM[5], s);
setDensity(ARGON, diffusiveDensity);
if (sw[15] != 0 && alt.getReal() <= altl[4]) {
/* Turbopause */
final double zh40 = PDM[4][2];
/* Mixed density at Zlb */
final T b40 = densu(zero.add(zh40), db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm40 = densu(alt, b40, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm40 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm40, zhm40, xmm, AR_MASS);
/* Correction to specified mixing ratio at ground */
final T rl = b28.multiply(PDM[4][1]).divide(b40).log();
final double hc40 = PDM[4][5] * PDL[1][9];
final double zc40 = PDM[4][4] * PDL[1][8];
/* Net density corrected at Alt */
setDensity(ARGON, diffusiveDensity.multiply(ccor(alt, rl, hc40, zc40)));
}
/**** H density ****/
/* Density variation factor at Zlb */
final T g1 = globe7(PD[6]).multiply(sw[21]);
/* Diffusive density at Zlb */
final T db01 = g1.exp().multiply(PDM[5][0] * PD[6][0]);
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db01, tinf, tlb, H_MASS, alpha[6], PTM[5], s);
setDensity(HYDROGEN, diffusiveDensity);
if (sw[15] != 0 && alt.getReal() <= altl[6]) {
/* Turbopause */
final double zh01 = PDM[5][2];
/* Mixed density at Zlb */
final T b01 = densu(zero.add(zh01), db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm01 = densu(alt, b01, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm01 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm01, zhm01, xmm, H_MASS);
/* Correction to specified mixing ratio at ground */
final T rl = b28.multiply(PDM[5][1] * FastMath.sqrt(PDL[1][17] * PDL[1][17])).divide(b01).log();
final double hc01 = PDM[5][5] * PDL[1][11];
final double zc01 = PDM[5][4] * PDL[1][10];
diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc01, zc01));
/* Chemistry correction */
final double hcc01 = PDM[5][7] * PDL[1][19];
final double zcc01 = PDM[5][6] * PDL[1][18];
final double rc01 = PDM[5][3] * PDL[1][20];
/* Net density corrected at Alt */
setDensity(HYDROGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc01), hcc01, zcc01)));
}
/**** N density ****/
/* Density variation factor at Zlb */
final T g14 = globe7(PD[7]).multiply(sw[21]);
/* Diffusive density at Zlb */
final T db14 = g14.exp().multiply(PDM[6][0] * PD[7][0]);
/* Diffusive density at Alt */
diffusiveDensity = densu(alt, db14, tinf, tlb, N_MASS, alpha[7], PTM[5], s);
setDensity(ATOMIC_NITROGEN, diffusiveDensity);
if (sw[15] != 0 && alt.getReal() <= altl[7]) {
/* Turbopause */
final double zh14 = PDM[6][2];
/* Mixed density at Zlb */
final T b14 = densu(zero.add(zh14), db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
/* Mixed density at Alt */
final T dm14 = densu(alt, b14, tinf, tlb, xmm, 0., PTM[5], s);
final double zhm14 = zhm28;
/* Net density at Alt */
diffusiveDensity = dnet(diffusiveDensity, dm14, zhm14, xmm, N_MASS);
/* Correction to specified mixing ratio at ground */
final T rl = b28.multiply(PDM[6][1] * PDL[0][2]).divide(b14).log();
final double hc14 = PDM[6][5] * PDL[0][1];
final double zc14 = PDM[6][4] * PDL[0][0];
diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc14, zc14));
/* Chemistry correction */
final double hcc14 = PDM[6][7] * PDL[0][4];
final double zcc14 = PDM[6][6] * PDL[0][3];
final double rc14 = PDM[6][3] * PDL[0][5];
/* Net density corrected at Alt */
setDensity(ATOMIC_NITROGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc14), hcc14, zcc14)));
}
/**** Anomalous O density ****/
final T g16h = globe7(PD[8]).multiply(sw[21]);
final T db16h = g16h.exp().multiply(PDM[7][0] * PD[8][0]);
final double tho = PDM[7][9] * PDL[0][6];
diffusiveDensity = densu(alt, db16h, zero.add(tho), zero.add(tho), O_MASS, alpha[8], PTM[5], s);
final double zsht = PDM[7][5];
final double zmho = PDM[7][4];
final T zsho = scalh(zmho, O_MASS, tho);
diffusiveDensity = diffusiveDensity.multiply(alt.negate().add(zmho).divide(zsht).exp().subtract(1).multiply(-zsht).divide(zsho).exp());
setDensity(ANOMALOUS_OXYGEN, diffusiveDensity);
// Convert densities from cm-3 to m-3
for (int i = 0; i < 9; i++) {
setDensity(i, getDensity(i).multiply(1.0e+06));
}
/**** Total mass density ****/
final T tmd = getDensity(HELIUM) .multiply(HE_MASS).
add(getDensity(ATOMIC_OXYGEN) .multiply( O_MASS)).
add(getDensity(MOLECULAR_NITROGEN).multiply(N2_MASS)).
add(getDensity(MOLECULAR_OXYGEN) .multiply(O2_MASS)).
add(getDensity(ARGON) .multiply(AR_MASS)).
add(getDensity(HYDROGEN) .multiply( H_MASS)).
add(getDensity(ATOMIC_NITROGEN) .multiply( N_MASS)).
multiply(AMU);
setDensity(TOTAL_MASS, tmd);
}
/** Calculate temperatures and densities not including anomalous oxygen.
* <p>NOTES ON INPUT VARIABLES:<br>
* Seconds, Local Time, and Longitude are used independently in the
* model and are not of equal importance for every situation.<br>
* For the most physically realistic calculation these three
* variables should be consistent (lst=sec/3600 + lon/15).<br>
* The Equation of Time departures from the above formula
* for apparent local time can be included if available but
* are of minor importance.<br><br>
*
* f107 and f107A values used to generate the model correspond
* to the 10.7 cm radio flux at the actual distance of the Earth
* from the Sun rather than the radio flux at 1 AU. The following
* site provides both classes of values:<br>
* ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
*
* f107, f107A, and ap effects are neither large nor well established below 80 km
* and these parameters should be set to 150., 150., and 4. respectively.
* </p>
* @param alt altitude (km)
*/
void gtd7(final T alt) {
// Calculates for thermosphere/mesosphere (above ZN2[0])
final T altt = (alt.getReal() > ZN2[0]) ? alt : zero.add(ZN2[0]);
gts7(altt);
if (alt.getReal() >= ZN2[0]) {
return;
}
// Calculates for lower mesosphere/upper stratosphere (between ZN2[0] and ZN3[0]):
// Temperature at nodes and gradients at end nodes
// Inverse temperature a linear function of spherical harmonics
final double r = PMA[2][0] * PAVGM[2];
meso_tgn2[0] = meso_tgn1[1];
meso_tn2[0] = meso_tn1[4];
meso_tn2[1] = glob7s(PMA[0]).multiply(sw[20] ).negate().add(1).reciprocal().multiply(PMA[0][0] * PAVGM[0]);
meso_tn2[2] = glob7s(PMA[1]).multiply(sw[20] ).negate().add(1).reciprocal().multiply(PMA[1][0] * PAVGM[1]);
meso_tn2[3] = glob7s(PMA[2]).multiply(sw[20] * sw[22]).negate().add(1).reciprocal().multiply(PMA[2][0] * PAVGM[2]);
meso_tgn2[1] = glob7s(PMA[9]).multiply(sw[20] * sw[22]).add(1).multiply(PMA[9][0] * PAVGM[8]).
multiply(meso_tn2[3]).multiply(meso_tn2[3]).divide(r * r);
meso_tn3[0] = meso_tn2[3];
// Calculates for lower stratosphere and troposphere (below ZN3[0])
// Temperature at nodes and gradients at end nodes
// Inverse temperature a linear function of spherical harmonics
if (alt.getReal() < ZN3[0]) {
final double q = PMA[6][0] * PAVGM[6];
meso_tgn3[0] = meso_tgn2[1];
meso_tn3[1] = glob7s(PMA[3]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[3][0] * PAVGM[3]);
meso_tn3[2] = glob7s(PMA[4]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[4][0] * PAVGM[4]);
meso_tn3[3] = glob7s(PMA[5]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[5][0] * PAVGM[5]);
meso_tn3[4] = glob7s(PMA[6]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[6][0] * PAVGM[6]);
meso_tgn3[1] = glob7s(PMA[7]).multiply(sw[22]) .add(1).multiply(PMA[7][0] * PAVGM[7]).
multiply(meso_tn3[4]).multiply(meso_tn3[4]).divide(q * q);
}
// Linear transition to full mixing below ZN2[0]
final T dmc = (alt.getReal() > ZMIX) ?
alt.subtract(ZN2[0]).divide(ZN2[0] - ZMIX).add(1) :
zero;
final T dz28 = getDensity(MOLECULAR_NITROGEN);
// N2 density
final T dm28m = dm28.multiply(1.0e+06);
T dmr = dz28.divide(dm28m).subtract(1);
T dst = densm(alt, dm28m, PDM[2][4]).multiply(dmr.multiply(dmc).add(1));
setDensity(MOLECULAR_NITROGEN, dst);
// HE density
dmr = getDensity(HELIUM).divide(dz28.multiply(PDM[0][1])).subtract(1);
dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[0][1]).multiply(dmr.multiply(dmc).add(1));
setDensity(HELIUM, dst);
// O density
setDensity(ATOMIC_OXYGEN, zero);
setDensity(ANOMALOUS_OXYGEN, zero);
// O2 density
dmr = getDensity(MOLECULAR_OXYGEN).divide(dz28.multiply(PDM[3][1])).subtract(1);
dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[3][1]).multiply(dmr.multiply(dmc).add(1));
setDensity(MOLECULAR_OXYGEN, dst);
// AR density
dmr = getDensity(ARGON).divide(dz28.multiply(PDM[4][1])).subtract(1);
dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[4][1]).multiply(dmr.multiply(dmc).add(1));
setDensity(ARGON, dst);
// H density
setDensity(HYDROGEN, zero);
// N density
setDensity(ATOMIC_NITROGEN, zero);
// Total mass density
final T tmd = getDensity(HELIUM) .multiply(HE_MASS).
add(getDensity(ATOMIC_OXYGEN) .multiply( O_MASS)).
add(getDensity(MOLECULAR_NITROGEN).multiply(N2_MASS)).
add(getDensity(MOLECULAR_OXYGEN) .multiply(O2_MASS)).
add(getDensity(ARGON) .multiply(AR_MASS)).
add(getDensity(HYDROGEN) .multiply( H_MASS)).
add(getDensity(ATOMIC_NITROGEN) .multiply( N_MASS)).
multiply(AMU);
setDensity(TOTAL_MASS, tmd);
// Temperature at altitude
setTemperature(ALTITUDE, densm(alt, field.getOne(), 0));
}
/** Calculate temperatures and densities including anomalous oxygen.
* <p></p>
* <p>NOTES ON INPUT VARIABLES:<br>
* Seconds, Local Time, and Longitude are used independently in the
* model and are not of equal importance for every situation.<br>
* For the most physically realistic calculation these three
* variables should be consistent (lst=sec/3600 + lon/15).<br>
* The Equation of Time departures from the above formula
* for apparent local time can be included if available but
* are of minor importance.<br>
* <br>
* f107 and f107A values used to generate the model correspond
* to the 10.7 cm radio flux at the actual distance of the Earth
* from the Sun rather than the radio flux at 1 AU. The following
* site provides both classes of values:<br>
* ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br>
* <br>
* f107, f107A, and ap effects are neither large nor well established below 80 km
* and these parameters should be set to 150., 150., and 4. respectively.
* </p>
* @param alt altitude (km)
*/
void gtd7d(final T alt) {
// Compute densities and temperatures
gtd7(alt);
// Update the total mass density with anomalous oxygen contribution
final T dTot = getDensity(TOTAL_MASS).add(getDensity(ANOMALOUS_OXYGEN).multiply( AMU * O_MASS));
setDensity(TOTAL_MASS, dTot);
}
/** Set one density.
* @param index one of the nine elements :
* <ul>
* <li>{@link #HELIUM}</li>
* <li>{@link #ATOMIC_OXYGEN}</li>
* <li>{@link #MOLECULAR_NITROGEN}</li>
* <li>{@link #MOLECULAR_OXYGEN}</li>
* <li>{@link #ARGON}</li>
* <li>{@link #TOTAL_MASS}</li>
* <li>{@link #HYDROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* </ul>
* @param d the value of density to set
*/
void setDensity(final int index, final T d) {
densities[index] = d;
}
/** Set one temperature.
* @param index one of the two elements :
* <ul>
* <li>{@link #EXOSPHERIC}</li>
* <li>{@link #ALTITUDE}</li>
* </ul>
* @param t the value of temperature to set
*/
void setTemperature(final int index, final T t) {
temperatures[index] = t;
}
/** Get one of the stored densities.
* @param index one of the nine elements :
* <ul>
* <li>{@link #HELIUM}</li>
* <li>{@link #ATOMIC_OXYGEN}</li>
* <li>{@link #MOLECULAR_NITROGEN}</li>
* <li>{@link #MOLECULAR_OXYGEN}</li>
* <li>{@link #ARGON}</li>
* <li>{@link #TOTAL_MASS}</li>
* <li>{@link #HYDROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* <li>{@link #ATOMIC_NITROGEN}</li>
* </ul>
* @return the requested density
*/
public T getDensity(final int index) {
return densities[index];
}
/** Calculate G(L) function with upper thermosphere parameters.
* @param p array of parameters
* @return G(L) value
*/
private T globe7(final double[] p) {
final T[] t = MathArrays.buildArray(field, 14);
final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
// F10.7 effect
final double df = f107 - f107a;
final double dfa = f107a - FLUX_REF;
t[0] = zero.add(p[19] * df * (1.0 + p[59] * dfa) +
p[20] * df * df +
p[21] * dfa +
p[29] * dfa * dfa);
final double f1 = 1.0 + (p[47] * dfa + p[19] * df + p[20] * df * df) * swc[1];
final double f2 = 1.0 + (p[49] * dfa + p[19] * df + p[20] * df * df) * swc[1];
// Time independent
t[1] = plg[0][2].multiply(p[ 1]).
add(plg[0][4].multiply(p[ 2])).
add(plg[0][6].multiply(p[22])).
add(plg[0][2].multiply(p[14] * dfa * swc[1])).
add(plg[0][1].multiply(p[26]));
// Symmetrical annual
t[2] = zero.add(p[18] * cd32);
// Symmetrical semiannual
t[3] = plg[0][2].multiply(p[16]).add(p[15]).multiply(cd18);
// Asymmetrical annual
t[4] = plg[0][1].multiply(p[9]).add(plg[0][3].multiply(p[10])).multiply(f1 * cd14);
// Asymmetrical semiannual
t[5] = plg[0][1].multiply(p[37] * cd39);
// Diurnal
if (sw[7] != 0) {
final T t71 = plg[1][2].multiply(p[11] * cd14 * swc[5]);
final T t72 = plg[1][2].multiply(p[12] * cd14 * swc[5]);
t[6] = plg[1][1].multiply(p[3]).add(plg[1][3].multiply(p[4])).add(plg[1][5].multiply(p[27])).add(t71).multiply(ctloc).
add(plg[1][1].multiply(p[6]).add(plg[1][3].multiply(p[7])).add(plg[1][5].multiply(p[28])).add(t72).multiply(stloc)).
multiply(f2);
}
// Semidiurnal
if (sw[8] != 0) {
final T t81 = plg[2][3].multiply(p[23]).add(plg[2][5].multiply(p[35])).multiply(cd14 * swc[5]);
final T t82 = plg[2][3].multiply(p[33]).add(plg[2][5].multiply(p[36])).multiply(cd14 * swc[5]);
t[7] = plg[2][2].multiply(p[5]).add(plg[2][4].multiply(p[41])).add(t81).multiply(c2tloc).
add(plg[2][2].multiply(p[8]).add(plg[2][4].multiply(p[42])).add(t82).multiply(s2tloc)).
multiply(f2);
}
// Terdiurnal
if (sw[14] != 0) {
t[13] = plg[3][3].multiply(p[39]).add(plg[3][4].multiply(p[93]).add(plg[3][6].multiply(p[46])).multiply(cd14 * swc[5])).multiply(s3tloc).
add(plg[3][3].multiply(p[40]).add(plg[3][4].multiply(p[94]).add(plg[3][6].multiply(p[48])).multiply(cd14 * swc[5])).multiply(c3tloc)).
multiply(f2);
}
// magnetic activity based on daily ap
if (sw[9] == -1) {
if (p[51] != 0) {
final T exp1 = lat.abs().negate().add(LAT_REF).multiply(p[138]).add(1).
reciprocal().multiply(-10800.0 * FastMath.abs(p[51])).
exp();
final double p24 = FastMath.max(p[24], 1.0e-4);
apt = sg0(min(0.99999, exp1), p24, p[25]);
t[8] = plg[0][2].multiply(p[96]).add(plg[0][4].multiply(p[54])).add(p[50]).
add((plg[0][1].multiply(p[125]).add(plg[0][3].multiply(p[126])).add(plg[0][5].multiply(p[127]))).multiply(cd14 * swc[5])).
add((plg[1][1].multiply(p[128]).add(plg[1][3].multiply(p[129])).add(plg[1][5].multiply(p[130]))).multiply(swc[7]).multiply(hl.subtract(p[131]).multiply(HOUR_TO_RAD).cos())).
multiply(apt);
}
} else {
final double apd = ap[0] - 4.0;
final double p44 = (p[43] < 0.) ? 1.0E-5 : p[43];
final double p45 = p[44];
apdf = apd + (p45 - 1.0) * (apd + (FastMath.exp(-p44 * apd) - 1.0) / p44);
if (sw[9] != 0) {
t[8] = plg[0][2].multiply(p[45]).add(plg[0][4].multiply(p[34])).add(p[32]).
add((plg[0][1].multiply(p[100]).add(plg[0][3].multiply(p[101])).add(plg[0][5].multiply(p[102]))).multiply(cd14 * swc[5])).
add((plg[1][1].multiply(p[121]).add(plg[1][3].multiply(p[122])).add(plg[1][5].multiply(p[123]))).multiply(swc[7]).multiply(hl.subtract(p[124]).multiply(HOUR_TO_RAD).cos())).
multiply(apdf);
}
}
if (sw[10] != 0) {
final T lonr = lon.multiply(DEG_TO_RAD);
// Longitudinal
if (sw[11] != 0) {
t[10] = plg[1][2].multiply(p[ 64]) .add(plg[1][4].multiply(p[ 65])).add(plg[1][6].multiply(p[ 66])).
add(plg[1][1].multiply(p[103])).add(plg[1][3].multiply(p[104])).add(plg[1][5].multiply(p[105])).
add((plg[1][1].multiply(p[109])).add(plg[1][3].multiply(p[110])).add(plg[1][5].multiply(p[111])).multiply(swc[5] * cd14)).
multiply(lonr.cos()).
add( plg[1][2].multiply(p[ 90]) .add(plg[1][4].multiply(p[ 91])).add(plg[1][6].multiply(p[ 92])).
add(plg[1][1].multiply(p[106])).add(plg[1][3].multiply(p[107])).add(plg[1][5].multiply(p[108])).
add((plg[1][1].multiply(p[112])).add(plg[1][3].multiply(p[113])).add(plg[1][5].multiply(p[114])).multiply(swc[5] * cd14)).
multiply(lonr.sin())).
multiply(1.0 + p[80] * dfa * swc[1]);
}
// ut and mixed ut, longitude
if (sw[12] != 0) {
t[11] = plg[0][1].multiply(p[95]).add(1).multiply(1.0 + p[81] * dfa * swc[1]).
multiply(plg[0][1].multiply(p[119] * swc[5] * cd14).add(1)).
multiply(plg[0][1].multiply(p[68]).add(plg[0][3].multiply(p[69])).add(plg[0][5].multiply(p[70]))).
multiply(sec.subtract(p[71]).multiply(SEC_TO_RAD).cos());
t[11] = t[11].
add(plg[2][3].multiply(p[76]).add(plg[2][5].multiply(p[77])).add(plg[2][7].multiply(p[78])).
multiply(swc[11] * (1.0 + p[137] * dfa * swc[1])).
multiply(sec.subtract(p[79]).multiply(SEC_TO_RAD).add(lonr.multiply(2)).cos()));
}
/* ut, longitude magnetic activity */
if (sw[13] != 0) {
if (sw[9] == -1) {
if (p[51] != 0.) {
t[12] = apt.multiply(swc[11]).multiply(plg[0][1].multiply(p[132]).add(1)).
multiply(plg[1][2].multiply(p[52]).add(plg[1][4].multiply(p[98])).add(plg[1][6].multiply(p[67]))).
multiply(lon.subtract(p[97]).multiply(DEG_TO_RAD).cos()).
add(apt.multiply(swc[11] * swc[5] * cd14).
multiply(plg[1][1].multiply(p[133]).add(plg[1][3].multiply(p[134])).add(plg[1][5].multiply(p[135]))).
multiply(lon.subtract(p[136]).multiply(DEG_TO_RAD).cos())).
add(apt.multiply(swc[12]).
multiply(plg[0][1].multiply(p[55]).add(plg[0][3].multiply(p[56])).add(plg[0][5].multiply(p[57]))).
multiply(sec.subtract(p[58]).multiply(SEC_TO_RAD).cos()));
}
} else {
t[12] = plg[0][1].multiply(p[120]).add(1).multiply(apdf * swc[11]).
multiply(plg[1][2].multiply(p[60]).add(plg[1][4].multiply(p[61])).add(plg[1][6].multiply(p[62]))).
multiply(lon.subtract(p[63]).multiply(DEG_TO_RAD).cos()).
add(plg[1][1].multiply(p[115]).add(plg[1][3].multiply(p[116])).add(plg[1][5].multiply(p[117])).
multiply(apdf * swc[11] * swc[5] * cd14).
multiply(lon.subtract(p[118]).multiply(DEG_TO_RAD).cos())).
add(plg[0][1].multiply(p[83]).add(plg[0][3].multiply(p[84])).add(plg[0][5].multiply(p[85])).
multiply(apdf * swc[12]).
multiply(sec.subtract(p[75]).multiply(SEC_TO_RAD).cos()));
}
}
}
// Sum all effects (params not used: 82, 89, 99, 139-149)
T tinf = zero.add(p[30]);
for (int i = 0; i < 14; i++) {
tinf = tinf.add(t[i].multiply(FastMath.abs(sw[i + 1])));
}
// Return G(L)
return tinf;
}
/** Calculate G(L) function with lower atmosphere parameters.
* @param p array of parameters
* @return G(L) value
*/
private T glob7s(final double[] p) {
final T[] t = MathArrays.buildArray(field, 14);
final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));
// F10.7 effect
t[0] = zero.add(p[21] * (f107a - FLUX_REF));
// Time independent
t[1] = plg[0][2].multiply(p[1]).
add(plg[0][4].multiply(p[2])).
add(plg[0][6].multiply(p[22])).
add(plg[0][1].multiply(p[26])).
add(plg[0][3].multiply(p[14])).
add(plg[0][5].multiply(p[59]));
// Symmetrical annual
t[2] = plg[0][2].multiply(p[47]).add(plg[0][4].multiply(p[29])).add(p[18]).multiply(cd32);
// Symmetrical semiannual
t[3] = plg[0][2].multiply(p[16]).add(plg[0][4].multiply(p[30])).add(p[15]).multiply(cd18);
// Asymmetrical annual
t[4] = plg[0][1].multiply(p[9]).add(plg[0][3].multiply(p[10])).add(plg[0][5].multiply(p[20])).multiply(cd14);
// Asymmetrical semiannual
t[5] = plg[0][1].multiply(p[37]).multiply(cd39);
// Diurnal
if (sw[7] != 0) {
final T t71 = plg[1][2].multiply(p[11]).multiply(cd14 * swc[5]);
final T t72 = plg[1][2].multiply(p[12]).multiply(cd14 * swc[5]);
t[6] = plg[1][1].multiply(p[3]).add(plg[1][3].multiply(p[4])).add(t71).multiply(ctloc).
add(plg[1][1].multiply(p[6]).add(plg[1][3].multiply(p[7])).add(t72).multiply(stloc));
}
// Semidiurnal
if (sw[8] != 0) {
final T t81 = plg[2][3].multiply(p[23]).add(plg[2][5].multiply(p[35])).multiply(cd14 * swc[5]);
final T t82 = plg[2][3].multiply(p[33]).add(plg[2][5].multiply(p[36])).multiply(cd14 * swc[5]);
t[7] = plg[2][2].multiply(p[5]).add(plg[2][4].multiply(p[41])).add(t81).multiply(c2tloc).
add(plg[2][2].multiply(p[8]).add(plg[2][4].multiply(p[42])).add(t82).multiply(s2tloc));
}
// Terdiurnal
if (sw[14] != 0) {
t[13] = plg[3][3].multiply(p[39]).multiply(s3tloc).add(plg[3][3].multiply(p[40]).multiply(c3tloc));
}
// Magnetic activity
if (sw[9] == 1) {
t[8] = plg[0][2].multiply(p[45] * swc[2]).add(p[32]).multiply(apdf);
} else if (sw[9] == -1) {
t[8] = plg[0][2].multiply(p[96] * swc[2]).add(p[50]).multiply(apt);
}
// Longitudinal
if (!(sw[10] == 0 || sw[11] == 0)) {
final T lonr = lon.multiply(DEG_TO_RAD);
t[10] = plg[0][1].multiply(p[80] * swc[5] * FastMath.cos(DAY_TO_RAD * (doy - p[81])) +
p[85] * swc[6] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[86]))).
add(1.0 +
p[83] * swc[3] * FastMath.cos(DAY_TO_RAD * (doy - p[84])) +
p[87] * swc[4] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[88]))).
multiply( plg[1][2].multiply(p[64]).
add(plg[1][4].multiply(p[65])).
add(plg[1][6].multiply(p[66])).
add(plg[1][1].multiply(p[74])).
add(plg[1][3].multiply(p[75])).
add(plg[1][5].multiply(p[76])).multiply(lonr.cos()).
add( plg[1][2].multiply(p[90]).
add(plg[1][4].multiply(p[91])).
add(plg[1][6].multiply(p[92])).
add(plg[1][1].multiply(p[77])).
add(plg[1][3].multiply(p[78])).
add(plg[1][5].multiply(p[79])).multiply(lonr.sin())));
}
// Sum all effects
T gl = zero;
for (int i = 0; i < 14; i++) {
gl = gl.add(t[i].multiply(FastMath.abs(sw[i + 1])));
}
// Return G(L)
return gl;
}
/** Implements sg0 function (Eq. A24a).
* @param ex ex
* @param p24 abs(p[24])
* @param p25 p[25]
* @return sg0
*/
private T sg0(final T ex, final double p24, final double p25) {
final double g01 = g0(ap[1], p24, p25);
final double g02 = g0(ap[2], p24, p25);
final double g03 = g0(ap[3], p24, p25);
final double g04 = g0(ap[4], p24, p25);
final double g05 = g0(ap[5], p24, p25);
final double g06 = g0(ap[6], p24, p25);
final T ex2 = ex.multiply(ex);
final T ex3 = ex.multiply(ex2);
final T ex4 = ex2.multiply(ex2);
final T ex8 = ex4.multiply(ex4);
final T ex12 = ex4.multiply(ex8);
final T g234 = ex.multiply(g02).add(ex2.multiply(g03)).add(ex3.multiply(g04));
final T g56 = ex4.multiply(g05).add(ex12.multiply(g06));
final T ex19 = ex3.multiply(ex4).multiply(ex12);
final T omex = ex.negate().add(1);
final T sumex = ex19.negate().add(1).divide(omex).multiply(ex.sqrt()).add(1);
return ex8.negate().add(1).multiply(g56).divide(omex).add(g234).add(g01).divide(sumex);
}
/** Implements go function (Eq. A24d).
* @param apI 3 hrs ap
* @param p24 abs(p[24])
* @param p25 p[25]
* @return go
*/
private double g0(final double apI, final double p24, final double p25) {
final double am4 = apI - 4.0;
return am4 + (p25 - 1.0) * (am4 + (FastMath.exp(-p24 * am4) - 1.0) / p24);
}
/** Calculates chemistry/dissociation correction for MSIS models.
* @param alt altitude
* @param r target ratio
* @param h1 transition scale length
* @param zh altitude of 1/2 R
* @return correction
*/
private T ccor(final T alt, final T r, final double h1, final double zh) {
final T e = alt.subtract(zh).divide(h1);
if (e.getReal() > 70.) {
return field.getOne();
} else if (e.getReal() < -70.) {
return r.exp();
} else {
return r.divide(e.exp().add(1)).exp();
}
}
/** Calculates O & O2 chemistry/dissociation correction for MSIS models.
* @param alt altitude
* @param r target ratio
* @param h1 transition scale length
* @param zh altitude of 1/2 R
* @param h2 transition scale length
* @return correction
*/
private T ccor2(final T alt, final double r, final double h1, final double zh, final double h2) {
final T e1 = alt.subtract(zh).divide(h1);
final T e2 = alt.subtract(zh).divide(h2);
if ((e1.getReal() > 70.) || (e2.getReal() > 70.)) {
return field.getOne();
} else if ((e1.getReal() < -70.) && (e2.getReal() < -70.)) {
return zero.add(FastMath.exp(r));
} else {
final T ex1 = e1.exp();
final T ex2 = e2.exp();
return ex1.add(ex2).multiply(0.5).add(1).reciprocal().multiply(r).exp();
}
}
/** Calculates scale height.
* @param alt altitude
* @param xm species molecular weight
* @param temp temperature
* @return scale height (km)
*/
private T scalh(final double alt, final double xm, final double temp) {
// Gravity at altitude
final T denom = rlat.reciprocal().multiply(alt).add(1);
final T galt = glat.divide(denom.multiply(denom));
return galt.reciprocal().multiply(R_GAS * temp / xm);
}
/** Calculates turbopause correction for MSIS models.
* @param dd diffusive density
* @param dm full mixed density
* @param zhm transition scale length
* @param xmm full mixed molecular weight
* @param xm species molecular weight
* @return combined density
*/
private T dnet(final T dd, final T dm, final double zhm, final double xmm, final double xm) {
if (!(dm.getReal() > 0 && dd.getReal() > 0)) {
T ddd = dd;
if (dd.getReal() == 0 && dm.getReal() == 0) {
ddd = field.getOne();
}
if (dm.getReal() == 0) {
return ddd;
}
if (dd.getReal() == 0) {
return dm;
}
}
final double a = zhm / (xmm - xm);
final T ylog = dm.divide(dd).log().multiply(a);
if (ylog.getReal() < -10.) {
return dd;
} else if (ylog.getReal() > 10.) {
return dm;
} else {
return ylog.exp().add(1).pow(1.0 / a).multiply(dd);
}
}
/** Integrate cubic spline function from xa[0] to x.
* <p>ADAPTED FROM NUMERICAL RECIPES</p>
* @param xa array of abscissas in ascending order
* @param ya array of ordinates in ascending order by xa
* @param y2a array of second derivatives in ascending order by xa
* @param x abscissa end point
* @return integral value
*/
private T splini(final T[] xa, final T[] ya, final T[] y2a, final T x) {
final int n = xa.length;
T yi = zero;
int klo = 0;
int khi = 1;
while (x.getReal() > xa[klo].getReal() && khi < n) {
T xx = x;
if (khi < n - 1) {
xx = (x.getReal() < xa[khi].getReal()) ? x : xa[khi];
}
final T h = xa[khi].subtract(xa[klo]);
final T a = xa[khi].subtract(xx).divide(h);
final T b = xx.subtract(xa[klo]).divide(h);
final T a2 = a.multiply(a);
final T b2 = b.multiply(b);
final T z =
a2.divide(2).subtract(a2.multiply(a2).add(1).divide(4)).multiply(y2a[klo]).
add(b2.multiply(b2).divide(4).subtract(b2.divide(2)).multiply(y2a[khi]));
yi = yi.add( a2.negate().add(1).multiply(ya[klo]).divide(2).
add(b2.multiply(ya[khi]).divide(2)).
add(z.multiply(h).multiply(h).divide(6)).
multiply(h));
klo++;
khi++;
}
return yi;
}
/** Calculate cubic spline interpolated value.
* <p>ADAPTED FROM NUMERICAL RECIPES</p>
* @param xa array of abscissas in ascending order
* @param ya array of ordinates in ascending order by xa
* @param y2a array of second derivatives in ascending order by xa
* @param x abscissa for interpolation
* @return interpolated value
*/
private T splint(final T[] xa, final T[] ya, final T[] y2a, final T x) {
final int n = xa.length;
int klo = 0;
int khi = n - 1;
while (khi - klo > 1) {
final int k = (khi + klo) >>> 1;
if (xa[k].getReal() > x.getReal()) {
khi = k;
} else {
klo = k;
}
}
final T h = xa[khi].subtract(xa[klo]);
final T a = xa[khi].subtract(x).divide(h);
final T b = x.subtract(xa[klo]).divide(h);
return a.multiply(ya[klo]).add(b.multiply(ya[khi])).
add(( a.multiply(a).multiply(a).subtract(a).multiply(y2a[klo]).
add(b.multiply(b).multiply(b).subtract(b).multiply(y2a[khi]))
).multiply(h).multiply(h).divide(6));
}
/** Calculate 2nd derivatives of cubic spline interpolation function.
* <p>ADAPTED FROM NUMERICAL RECIPES</p>
* @param x array of abscissas in ascending order
* @param y array of ordinates in ascending order by x
* @param yp1 derivative at x[0] (2nd derivatives null if > 1E30)
* @param ypn derivative at x[n-1] (2nd derivatives null if > 1E30)
* @return array of second derivatives
*/
private T[] spline(final T[] x, final T[] y, final T yp1, final T ypn) {
final int n = x.length;
final T[] y2 = MathArrays.buildArray(field, n);
final T[] u = MathArrays.buildArray(field, n);
if (yp1.getReal() < 1e+30) {
y2[0] = zero.add(-0.5);
final T dx = x[1].subtract(x[0]);
final T dy = y[1].subtract(y[0]);
u[0] = dx.reciprocal().multiply(3.0).multiply(dy.divide(dx).subtract(yp1));
}
for (int i = 1; i < n - 1; i++) {
final T dx0m = x[i].subtract(x[i - 1]);
final T dy0m = y[i].subtract(y[i - 1]);
final T dxpm = x[i + 1].subtract(x[i - 1]);
final T dxp0 = x[i + 1].subtract(x[i]);
final T dyp0 = y[i + 1].subtract(y[i]);
final T sig = dx0m.divide(dxpm);
final T p = sig.multiply(y2[i - 1]).add(2.0);
y2[i] = sig.subtract(1.0).divide(p);
u[i] = dyp0.divide(dxp0).subtract(dy0m.divide(dx0m)).multiply(6).divide(dxpm).subtract(sig.multiply(u[i - 1])).divide(p);
}
double qn = 0;
T un = zero;
if (ypn.getReal() < 1e+30) {
final T dx12 = x[n - 1].subtract(x[n - 2]);
final T dy12 = y[n - 1].subtract(y[n - 2]);
qn = 0.5;
un = dx12.reciprocal().multiply(3.0).multiply(ypn.subtract(dy12.divide(dx12)));
}
y2[n - 1] = un.subtract(u[n - 2].multiply(qn)).divide(y2[n - 2].multiply(qn).add(1.0));
for (int k = n - 2; k >= 0; k--) {
y2[k] = y2[k].multiply(y2[k + 1]).add(u[k]);
}
return y2;
}
/** Calculate Temperature and Density Profiles for lower atmosphere.
* @param alt altitude
* @param d0 density
* @param xm mixed density
* @return temperature or density profile
*/
private T densm(final T alt, final T d0, final double xm) {
T densm = d0;
// stratosphere/mesosphere temperature
int mn = ZN2.length;
T z = (alt.getReal() > ZN2[mn - 1]) ? alt : zero.add(ZN2[mn - 1]);
double z1 = ZN2[0];
double z2 = ZN2[mn - 1];
T t1 = meso_tn2[0];
T t2 = meso_tn2[mn - 1];
T zg = zeta(z, z1);
T zgdif = zeta(zero.add(z2), z1);
/* set up spline nodes */
T[] xs = MathArrays.buildArray(field, mn);
T[] ys = MathArrays.buildArray(field, mn);
for (int k = 0; k < mn; k++) {
xs[k] = zeta(zero.add(ZN2[k]), z1).divide(zgdif);
ys[k] = meso_tn2[k].reciprocal();
}
final T qSM = rlat.add(z2).divide(rlat.add(z1));
T yd1 = meso_tgn2[0].negate().divide(t1.multiply(t1)).multiply(zgdif);
T yd2 = meso_tgn2[1].negate().divide(t2.multiply(t2)).multiply(zgdif).multiply(qSM).multiply(qSM);
/* calculate spline coefficients */
T[] y2out = spline(xs, ys, yd1, yd2);
T x = zg.divide(zgdif);
T y = splint(xs, ys, y2out, x);
/* temperature at altitude */
T tz = y.reciprocal();
if (xm != 0.0) {
/* calculate stratosphere / mesospehere density */
final T glb = galt(zero.add(z1));
final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);
/* Integrate temperature profile */
final T yi = splini(xs, ys, y2out, x);
final T expl = min(50., gamm.multiply(yi));
/* Density at altitude */
densm = densm.multiply(t1.divide(tz).multiply(expl.negate().exp()));
}
if (alt.getReal() > ZN3[0]) {
return (xm == 0.0) ? tz : densm;
}
// troposhere/stratosphere temperature
z = alt;
mn = ZN3.length;
z1 = ZN3[0];
z2 = ZN3[mn - 1];
t1 = meso_tn3[0];
t2 = meso_tn3[mn - 1];
zg = zeta(z, z1);
zgdif = zeta(zero.add(z2), z1);
/* set up spline nodes */
xs = MathArrays.buildArray(field, mn);
ys = MathArrays.buildArray(field, mn);
for (int k = 0; k < mn; k++) {
xs[k] = zeta(zero.add(ZN3[k]), z1).divide(zgdif);
ys[k] = meso_tn3[k].reciprocal();
}
final T qTS = rlat.add(z2) .divide(rlat.add(z1));
yd1 = meso_tgn3[0].negate().divide(t1.multiply(t1)).multiply(zgdif);
yd2 = meso_tgn3[1].negate().divide(t2.multiply(t2)).multiply(zgdif).multiply(qTS).multiply(qTS);
/* calculate spline coefficients */
y2out = spline(xs, ys, yd1, yd2);
x = zg.divide(zgdif);
y = splint(xs, ys, y2out, x);
/* temperature at altitude */
tz = y.reciprocal();
if (xm != 0.0) {
/* calculate tropospheric / stratosphere density */
final T glb = galt(zero.add(z1));
final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);
/* Integrate temperature profile */
final T yi = splini(xs, ys, y2out, x);
final T expl = min(50., gamm.multiply(yi));
/* Density at altitude */
densm = densm.multiply(t1.divide(tz).multiply(expl.negate().exp()));
}
return (xm == 0.0) ? tz : densm;
}
/** Calculate temperature and density profiles according to new lower thermo polynomial.
* @param alt altitude
* @param dlb density at lower boundary
* @param tinf exospheric temperature
* @param tlb temperature at lower boundary
* @param xm species molecular weight
* @param alpha thermal diffusion coefficient
* @param zlb altitude of the lower boundary
* @param s2 slope
* @return temperature or density profile
*/
private T densu(final T alt, final T dlb, final T tinf,
final T tlb, final double xm, final double alpha,
final double zlb, final T s2) {
/* joining altitudes of Bates and spline */
T z = (alt.getReal() > ZN1[0]) ? alt : zero.add(ZN1[0]);
/* geopotential altitude difference from ZLB */
final T zg2 = zeta(z, zlb);
/* Bates temperature */
final T tt = tinf.subtract(tinf.subtract(tlb).multiply(s2.negate().multiply(zg2).exp()));
final T ta = tt;
T tz = tt;
final int mn = ZN1.length;
final T[] xs = MathArrays.buildArray(field, mn);
final T[] ys = MathArrays.buildArray(field, mn);
T x = zero;
T[] y2out = MathArrays.buildArray(field, mn);
T zgdif = zero;
if (alt.getReal() < ZN1[0]) {
/* calculate temperature below ZA
* temperature gradient at ZA from Bates profile */
final T p = rlat.add(zlb).divide(rlat.add(ZN1[0]));
final T dta = tinf.subtract(ta).multiply(s2).multiply(p.multiply(p));
meso_tgn1[0] = dta;
meso_tn1[0] = ta;
final T tzn1mn1 = zero.add(ZN1[mn - 1]);
z = (alt.getReal() > ZN1[mn - 1]) ? alt : tzn1mn1;
final T t1 = meso_tn1[0];
final T t2 = meso_tn1[mn - 1];
/* geopotental difference from z1 */
final T zg = zeta(z, ZN1[0]);
zgdif = zeta(tzn1mn1, ZN1[0]);
/* set up spline nodes */
for (int k = 0; k < mn; k++) {
xs[k] = zeta(zero.add(ZN1[k]), ZN1[0]).divide(zgdif);
ys[k] = meso_tn1[k].reciprocal();
}
/* end node derivatives */
final T q = rlat.add(ZN1[mn - 1]).divide(rlat.add(ZN1[0]));
final T yd1 = meso_tgn1[0].negate().divide(t1.multiply(t1)).multiply(zgdif);
final T yd2 = meso_tgn1[1].negate().divide(t2.multiply(t2)).multiply(zgdif).multiply(q.multiply(q));
/* calculate spline coefficients */
y2out = spline(xs, ys, yd1, yd2);
x = zg.divide(zgdif);
final T y = splint(xs, ys, y2out, x);
/* temperature at altitude */
tz = y.reciprocal();
}
if (xm == 0) {
return tz;
}
/* calculate density above za */
T glb = galt(zero.add(zlb));
T gamma = glb.divide(s2.multiply(tinf)).multiply(xm / R_GAS);
T expl = tt.getReal() <= 0 ?
zero.add(50) :
min(50.0, s2.negate().multiply(gamma).multiply(zg2).exp());
T densu = dlb.multiply(tlb.divide(tt).pow(gamma.add(alpha + 1))).multiply(expl);
/* calculate density below za */
if (alt.getReal() < ZN1[0]) {
glb = galt(zero.add(ZN1[0]));
gamma = glb.multiply(zgdif).multiply(xm / R_GAS);
/* integrate spline temperatures */
expl = tz.getReal() <= 0 ?
zero.add(50.0) :
min(50.0, gamma.multiply(splini(xs, ys, y2out, x)));
/* correct density at altitude */
densu = densu.multiply(meso_tn1[0].divide(tz).pow(alpha + 1).multiply(expl.negate().exp()));
}
/* Return density at altitude */
return densu;
}
/** Compute min of two values, one double and one field element.
* @param d double value
* @param f field element
* @return min value
*/
private T min(final double d, final T f) {
return (f.getReal() > d) ? zero.add(d) : f;
}
/** Calculate gravity at altitude.
* @param alt altitude (km)
* @return gravity at altitude (cm/s2)
*/
private T galt(final T alt) {
final T r = alt.divide(rlat).add(1);
return glat.divide(r.multiply(r));
}
/** Calculate zeta function.
* @param zz zz value
* @param zl zl value
* @return value of zeta function
*/
private T zeta(final T zz, final double zl) {
return zz.subtract(zl).multiply(rlat.add(zl)).divide(rlat.add(zz));
}
}
}