/* 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.RealFieldElement; import org.hipparchus.geometry.euclidean.threed.FieldVector3D; import org.hipparchus.geometry.euclidean.threed.Vector3D; import org.hipparchus.util.FastMath; import org.orekit.bodies.BodyShape; 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.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> * * @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; /** Error message when index is wrong. */ private static final String WRONG_INDEX = "Index is expected to be between 1 and 23"; // 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; /** Input switches. */ private int[] switches = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; /** Switches for main effects. */ private int[] sw = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; /** Switches for cross effects. */ private int[] swc = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; /** Gravity at latitude (cm/s2). */ private double glat; /** Effective Earth radius at latitude (km). */ private double rlat; /** N2 mixed density at alt. */ private double dm28; /** Legendre polynomials. */ private double[][] plg = new double[4][8]; /** Cosinus of local solar time. */ private double ctloc; /** Sinus of local solar time. */ private double stloc; /** Square of ctloc. */ private double c2tloc; /** Square of stloc. */ private double s2tloc; /** Cube of ctloc. */ private double c3tloc; /** Cube of stloc. */ private 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 double[] meso_tn1 = new double[ZN1.length]; /** Temperature at nodes for ZN2 scale. */ private double[] meso_tn2 = new double[ZN2.length]; /** Temperature at nodes for ZN3 scale. */ private double[] meso_tn3 = new double[ZN3.length]; /** Temperature gradients at end nodes for ZN1 scale. */ private double[] meso_tgn1 = new double[2]; /** Temperature gradients at end nodes for ZN2 scale. */ private double[] meso_tgn2 = new double[2]; /** Temperature gradients at end nodes for ZN3 scale. */ private double[] meso_tgn3 = new double[2]; /** 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.inputParams = parameters; this.sun = sun; this.earth = earth; } /** {@inheritDoc} */ public Frame getFrame() { return earth.getBodyFrame(); } /** {@inheritDoc} */ 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 = gtd7d(doy, sec, alt, lat, lon, lst, inputParams.getAverageFlux(date), inputParams.getDailyFlux(date), inputParams.getAp(date)); // return the local density return out.getDensity(Output.TOTAL_MASS); } @Override public <T extends RealFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position, final Frame frame) throws OrekitException { // TODO: field implementation throw new UnsupportedOperationException(); } /** 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 the current switches. * @return the array of switches */ public int[] getSwitches() { return Arrays.copyOfRange(switches, 1, switches.length); } /** Get a specific switch. * @param number the number in the array of switches (between 1 and 23) * @return the switch at the given index */ public int getSwitch(final int number) { if (number < 1 || number > 23) { throw new IllegalArgumentException(WRONG_INDEX); } return switches[number]; } /** Set all the switches. * @param switches the array of 23 switches */ public void setSwitches(final int[] switches) { if (switches.length != 23) { throw new IllegalArgumentException("The array of switches is expected to have 23 elements"); } for (int i = 0; i < switches.length; i++) { setSwitch(i + 1, switches[i]); } } /** Set a value a specific switch. * <p> * For switches from #1 to #8 and #10 to #23, any value not equal to 1 or 2 is set to 0.<br> * For switch #9, any value not equal to -1 or 1 is set to 0. * </p> * @param number the number of the switch to be set (between 1 and 23) * @param value the value to set */ public void setSwitch(final int number, final int value) { if (number < 1 || number > 23) { throw new IllegalArgumentException(WRONG_INDEX); } this.switches[number] = value; if (number != 9) { this.sw[number] = (value == 1) ? 1 : 0; this.swc[number] = (value > 0) ? 1 : 0; } else { if (value == -1 || value == 1) { this.sw[number] = value; } else { this.sw[number] = 0; } this.swc[number] = this.sw[number]; } } /** 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 doy day of year (from 1 to 365 or 366) * @param sec seconds in day (UT scale) * @param alt altitude (km) * @param lat geodetic latitude (°) * @param lon geodetic longitude (°) * @param hl local apparent solar time (hours) * @param f107 daily F10.7 flux for previous day * @param f107a 81 day average of F10.7 flux (centered on 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> * @return {@link Output output data} */ public Output gtd7d(final int doy, final double sec, final double alt, final double lat, final double lon, final double hl, final double f107a, final double f107, final double[] ap) { // Compute densities and temperatures final Output out = gtd7(doy, sec, alt, lat, lon, hl, f107a, f107, ap); // Update the total mass density with anomalous oxygen contribution final double dTot = out.getDensity(Output.TOTAL_MASS) + AMU * O_MASS * out.getDensity(Output.ANOMALOUS_OXYGEN); out.setDensity(Output.TOTAL_MASS, dTot); // Return the updated output return out; } /** 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 doy day of year (from 1 to 365 or 366) * @param sec seconds in day (UT scale) * @param alt altitude (km) * @param lat geodetic latitude (°) * @param lon geodetic longitude (°) * @param hl local apparent solar time (hours) * @param f107 daily F10.7 flux for previous day * @param f107a 81 day average of F10.7 flux (centered on 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> * @return {@link Output output data} */ public Output gtd7(final int doy, final double sec, final double alt, final double lat, final double lon, final double hl, final double f107a, final double f107, final double[] ap) { // Calculates latitude variable gravity and effective radius final double xlat = (sw[2] == 0) ? LAT_REF : lat; glatf(xlat); // Calculates for thermosphere/mesosphere (above ZN2[0]) final double altt = (alt > ZN2[0]) ? alt : ZN2[0]; final Output output = gts7(doy, sec, altt, lat, lon, hl, f107a, f107, ap); if (alt >= ZN2[0]) { return output; } // 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 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], doy, lon, f107a)); meso_tn2[2] = PMA[1][0] * PAVGM[1] / (1.0 - sw[20] * glob7s(PMA[1], doy, lon, f107a)); meso_tn2[3] = PMA[2][0] * PAVGM[2] / (1.0 - sw[20] * sw[22] * glob7s(PMA[2], doy, lon, f107a)); meso_tgn2[1] = PMA[9][0] * PAVGM[8] * (1.0 + sw[20] * sw[22] * glob7s(PMA[9], doy, lon, f107a)) * meso_tn2[3] * meso_tn2[3] / FastMath.pow(PMA[2][0] * PAVGM[2], 2); 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]) { meso_tgn3[0] = meso_tgn2[1]; meso_tn3[1] = PMA[3][0] * PAVGM[3] / (1.0 - sw[22] * glob7s(PMA[3], doy, lon, f107a)); meso_tn3[2] = PMA[4][0] * PAVGM[4] / (1.0 - sw[22] * glob7s(PMA[4], doy, lon, f107a)); meso_tn3[3] = PMA[5][0] * PAVGM[5] / (1.0 - sw[22] * glob7s(PMA[5], doy, lon, f107a)); meso_tn3[4] = PMA[6][0] * PAVGM[6] / (1.0 - sw[22] * glob7s(PMA[6], doy, lon, f107a)); meso_tgn3[1] = PMA[7][0] * PAVGM[7] * (1.0 + sw[22] * glob7s(PMA[7], doy, lon, f107a)) * meso_tn3[4] * meso_tn3[4] / FastMath.pow(PMA[6][0] * PAVGM[6], 2); } // 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 = output.getDensity(Output.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); output.setDensity(Output.MOLECULAR_NITROGEN, dst); // HE density dmr = output.getDensity(Output.HELIUM) / (dz28 * PDM[0][1]) - 1.0; dst = output.getDensity(Output.MOLECULAR_NITROGEN) * PDM[0][1] * (1.0 + dmr * dmc); output.setDensity(Output.HELIUM, dst); // O density output.setDensity(Output.ATOMIC_OXYGEN, 0.); output.setDensity(Output.ANOMALOUS_OXYGEN, 0.); // O2 density dmr = output.getDensity(Output.MOLECULAR_OXYGEN) / (dz28 * PDM[3][1]) - 1.0; dst = output.getDensity(Output.MOLECULAR_NITROGEN) * PDM[3][1] * (1.0 + dmr * dmc); output.setDensity(Output.MOLECULAR_OXYGEN, dst); // AR density dmr = output.getDensity(Output.ARGON) / (dz28 * PDM[4][1]) - 1.0; dst = output.getDensity(Output.MOLECULAR_NITROGEN) * PDM[4][1] * (1.0 + dmr * dmc); output.setDensity(Output.ARGON, dst); // H density output.setDensity(Output.HYDROGEN, 0.); // N density output.setDensity(Output.ATOMIC_NITROGEN, 0.); // Total mass density final double tmd = AMU * (HE_MASS * output.getDensity(Output.HELIUM) + O_MASS * output.getDensity(Output.ATOMIC_OXYGEN) + N2_MASS * output.getDensity(Output.MOLECULAR_NITROGEN) + O2_MASS * output.getDensity(Output.MOLECULAR_OXYGEN) + AR_MASS * output.getDensity(Output.ARGON) + H_MASS * output.getDensity(Output.HYDROGEN) + N_MASS * output.getDensity(Output.ATOMIC_NITROGEN)); output.setDensity(Output.TOTAL_MASS, tmd); // Temperature at altitude output.setTemperature(Output.ALTITUDE, densm(alt, 1.0, 0)); return output; } /** Calculates latitude variable gravity and effective radius. * @param lat latitude (°) */ private void glatf(final double lat) { final double latr = DEG_TO_RAD * lat; final double c2 = FastMath.cos(2 * latr); glat = G_REF * (1. - .0026373 * c2); rlat = 2. * glat / (3.085462e-6 + 2.27e-9 * c2) * 1.e-5; } /** 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 doy day of year (from 1 to 365 or 366) * @param sec seconds in day (UT scale) * @param alt altitude (km) * @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> * @return {@link Output output data} */ public Output gts7(final int doy, final double sec, final double alt, final double lat, final double lon, final double hl, final double f107a, final double f107, final double[] ap) { // 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]; // Create output final Output output = new Output(); // Calculate Legendre polynomials and related data initialize(lat, hl); /**** 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, doy, sec, alt, lat, lon, hl, f107a, f107, ap); } output.setTemperature(Output.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, doy, sec, alt, lat, lon, hl, f107a, f107, ap); } // Temperature at lower boundary double tlb = PTM[1] * PD[3][0]; tlb *= 1.0 + sw[17] * globe7(PD[3], doy, sec, alt, lat, lon, hl, f107a, f107, ap); // 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) { meso_tn1[1] /= 1.0 - sw[18] * glob7s(PTL[0], doy, lon, f107a); meso_tn1[2] /= 1.0 - sw[18] * glob7s(PTL[1], doy, lon, f107a); meso_tn1[3] /= 1.0 - sw[18] * glob7s(PTL[2], doy, lon, f107a); meso_tn1[4] /= 1.0 - sw[18] * sw[20] * glob7s(PTL[3], doy, lon, f107a); meso_tgn1[1] *= 1.0 + sw[18] * sw[20] * glob7s(PMA[8], doy, lon, f107a); meso_tgn1[1] *= meso_tn1[4] * meso_tn1[4] / FastMath.pow(PTM[4] * PTL[3][0], 2); } /**** Temperature at altitude ****/ output.setTemperature(Output.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], doy, sec, alt, lat, lon, hl, f107a, f107, ap); /* Diffusive density at Zlb */ final double db28 = PDM[2][0] * FastMath.exp(g28) * PD[2][0]; /* Diffusive density at Alt */ double dd = densu(alt, db28, tinf, tlb, N2_MASS, alpha[2], PTM[5], s); output.setDensity(Output.MOLECULAR_NITROGEN, dd); // 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 */ output.setDensity(Output.MOLECULAR_NITROGEN, dnet(dd, dm28, zhm28, xmm, N2_MASS)); } /**** He density ****/ /* Density variation factor at Zlb */ final double g4 = sw[21] * globe7(PD[0], doy, sec, alt, lat, lon, hl, f107a, f107, ap); /* Diffusive density at Zlb */ final double db04 = PDM[0][0] * FastMath.exp(g4) * PD[0][0]; /* Diffusive density at Alt */ dd = densu(alt, db04, tinf, tlb, HE_MASS, alpha[0], PTM[5], s); output.setDensity(Output.HELIUM, dd); 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 */ dd = dnet(dd, 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 */ output.setDensity(Output.HELIUM, dd * ccor(alt, rl, hc04, zc04)); } /**** O density ****/ /* Density variation factor at Zlb */ final double g16 = sw[21] * globe7(PD[1], doy, sec, alt, lat, lon, hl, f107a, f107, ap); /* Diffusive density at Zlb */ final double db16 = PDM[1][0] * FastMath.exp(g16) * PD[1][0]; /* Diffusive density at Alt */ dd = densu(alt, db16, tinf, tlb, O_MASS, alpha[1], PTM[5], s); output.setDensity(Output.ATOMIC_OXYGEN, dd); 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 */ dd = dnet(dd, 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]; dd *= 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 */ output.setDensity(Output.ATOMIC_OXYGEN, dd * ccor(alt, rc16, hcc16, zcc16)); } /**** O2 density ****/ /* Density variation factor at Zlb */ final double g32 = sw[21] * globe7(PD[4], doy, sec, alt, lat, lon, hl, f107a, f107, ap); /* Diffusive density at Zlb */ final double db32 = PDM[3][0] * FastMath.exp(g32) * PD[4][0]; /* Diffusive density at Alt */ dd = densu(alt, db32, tinf, tlb, O2_MASS, alpha[3], PTM[5], s); output.setDensity(Output.MOLECULAR_OXYGEN, dd); 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 */ dd = dnet(dd, 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]; dd *= 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 */ output.setDensity(Output.MOLECULAR_OXYGEN, dd * ccor2(alt, rc32, hcc32, zcc32, hcc232)); } /**** Ar density ****/ /* Density variation factor at Zlb */ final double g40 = sw[21] * globe7(PD[5], doy, sec, alt, lat, lon, hl, f107a, f107, ap); /* Diffusive density at Zlb */ final double db40 = PDM[4][0] * FastMath.exp(g40) * PD[5][0]; /* Diffusive density at Alt */ dd = densu(alt, db40, tinf, tlb, AR_MASS, alpha[4], PTM[5], s); output.setDensity(Output.ARGON, dd); 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 */ dd = dnet(dd, 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 */ output.setDensity(Output.ARGON, dd * ccor(alt, rl, hc40, zc40)); } /**** H density ****/ /* Density variation factor at Zlb */ final double g1 = sw[21] * globe7(PD[6], doy, sec, alt, lat, lon, hl, f107a, f107, ap); /* Diffusive density at Zlb */ final double db01 = PDM[5][0] * FastMath.exp(g1) * PD[6][0]; /* Diffusive density at Alt */ dd = densu(alt, db01, tinf, tlb, H_MASS, alpha[6], PTM[5], s); output.setDensity(Output.HYDROGEN, dd); 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 */ dd = dnet(dd, 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]; dd *= 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 */ output.setDensity(Output.HYDROGEN, dd * ccor(alt, rc01, hcc01, zcc01)); } /**** N density ****/ /* Density variation factor at Zlb */ final double g14 = sw[21] * globe7(PD[7], doy, sec, alt, lat, lon, hl, f107a, f107, ap); /* Diffusive density at Zlb */ final double db14 = PDM[6][0] * FastMath.exp(g14) * PD[7][0]; /* Diffusive density at Alt */ dd = densu(alt, db14, tinf, tlb, N_MASS, alpha[7], PTM[5], s); output.setDensity(Output.ATOMIC_NITROGEN, dd); 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 */ dd = dnet(dd, 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]; dd *= 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 */ output.setDensity(Output.ATOMIC_NITROGEN, dd * ccor(alt, rc14, hcc14, zcc14)); } /**** Anomalous O density ****/ final double g16h = sw[21] * globe7(PD[8], doy, sec, alt, lat, lon, hl, f107a, f107, ap); final double db16h = PDM[7][0] * FastMath.exp(g16h) * PD[8][0]; final double tho = PDM[7][9] * PDL[0][6]; dd = 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); dd *= FastMath.exp(-zsht / zsho * (FastMath.exp((zmho - alt ) / zsht) - 1.)); output.setDensity(Output.ANOMALOUS_OXYGEN, dd); // Convert densities from cm-3 to m-3 for (int i = 0; i < 9; i++) { output.setDensity(i, output.getDensity(i) * 1.0e+06); } /**** Total mass density ****/ final double tmd = AMU * (HE_MASS * output.getDensity(Output.HELIUM) + O_MASS * output.getDensity(Output.ATOMIC_OXYGEN) + N2_MASS * output.getDensity(Output.MOLECULAR_NITROGEN) + O2_MASS * output.getDensity(Output.MOLECULAR_OXYGEN) + AR_MASS * output.getDensity(Output.ARGON) + H_MASS * output.getDensity(Output.HYDROGEN) + N_MASS * output.getDensity(Output.ATOMIC_NITROGEN)); output.setDensity(Output.TOTAL_MASS, tmd); // Return output data return output; } /** Calculate Legendre polynomials and related data. * @param lat latitude (°) * @param hl local apparent solar time (hours) */ private void initialize(final double lat, final double hl) { // 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); } } /** Calculate G(L) function with upper thermosphere parameters. * @param p array of parameters * @param doy day of year (from 1 to 365 or 366) * @param sec seconds in day (UT scale) * @param alt altitude (km) * @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 of ap indices * @return G(L) value */ private double globe7(final double[] p, final int doy, final double sec, final double alt, final double lat, final double lon, final double hl, final double f107a, final double f107, final double[] ap) { 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], ap); 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 * @param doy day of year (from 1 to 365 or 366) * @param lon geodetic longitude (°) * @param f107a 81 day average of F10.7 flux (centered on day) * @return G(L) value */ private double glob7s(final double[] p, final int doy, final double lon, final double f107a) { 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 tt = 0;; for (int i = 0; i < 14; i++) { tt += FastMath.abs(sw[i + 1]) * t[i]; } // Return G(L) return tt; } /** Implements sg0 function (Eq. A24a). * @param ex ex * @param p24 abs(p[24]) * @param p25 p[25] * @param ap array of ap indices * @return sg0 */ private double sg0(final double ex, final double p24, final double p25, final double[] ap) { 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 ap 3 hrs ap * @param p24 abs(p[24]) * @param p25 p[25] * @return go */ private double g0(final double ap, final double p24, final double p25) { final double am4 = ap - 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]; } double yd1 = -meso_tgn2[0] / (t1 * t1) * zgdif; double yd2 = -meso_tgn2[1] / (t2 * t2) * zgdif * FastMath.pow((rlat + z2) / (rlat + z1), 2); /* 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]; } yd1 = -meso_tgn3[0] / (t1 * t1) * zgdif; yd2 = -meso_tgn3[1] / (t2 * t2) * zgdif * FastMath.pow((rlat + z2) / (rlat + z1), 2); /* 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 dta = (tinf - ta) * s2 * FastMath.pow((rlat + zlb) / (rlat + ZN1[0]), 2); 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 yd1 = -meso_tgn1[0] / (t1 * t1) * zgdif; final double yd2 = -meso_tgn1[1] / (t2 * t2) * zgdif * FastMath.pow((rlat + ZN1[mn - 1]) / (rlat + ZN1[0]), 2); /* 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) { return glat / FastMath.pow(1.0 + alt / rlat, 2); } /** 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> */ public static class Output { // Constants /** Identifier for helium density. */ public static final int HELIUM = 0; /** Identifier for atomic oxygen density. */ public static final int ATOMIC_OXYGEN = 1; /** Identifier for molecular nitrogen density. */ public static final int MOLECULAR_NITROGEN = 2; /** Identifier for molecular oxygen density. */ public static final int MOLECULAR_OXYGEN = 3; /** Identifier for argon density. */ public static final int ARGON = 4; /** Identifier for atomic nitrogen density. */ public static final int TOTAL_MASS = 5; /** Identifier for hydrogen density. */ public static final int HYDROGEN = 6; /** Identifier for atomic nitrogen density. */ public static final int ATOMIC_NITROGEN = 7; /** Identifier for anomalous oxygen density. */ public static final int ANOMALOUS_OXYGEN = 8; /** Identifier for exospheric temperature. */ public static final int EXOSPHERIC = 0; /** Identifier for temperature at altitude. */ public static final int ALTITUDE = 1; /** Error message for wrong index. */ private static final String WRONG_DIMENSION = "Unexpected array length"; /** Error message for wrong index. */ private static final String WRONG_INDEX = "Unexpected index"; // Fields /** Densities. */ private double[] dd = new double[9]; /** Temperatures. */ private double[] tt = new double[2]; /** Simple constructor. */ public Output() { } /** Constructor. * @param d the densities as an array of 9 values * @param t the temperatures as an array of 2 values */ public Output(final double[] d, final double[] t) { dd = d.clone(); tt = t.clone(); } /** Set the densities. * @param d the densities as an array of 9 values */ public void setDensities(final double[] d) { if (d.length != 9) { throw new IllegalArgumentException(WRONG_DIMENSION); } dd = d.clone(); } /** 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 */ public void setDensity(final int index, final double d) { if (index < 0 || index > 8) { throw new IllegalArgumentException(WRONG_INDEX); } dd[index] = d; } /** Set the temperatures. * @param t the temperatures as an array of 2 values */ public void setTemperatures(final double[] t) { if (t.length != 2) { throw new IllegalArgumentException(WRONG_DIMENSION); } tt = t.clone(); } /** 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 */ public void setTemperature(final int index, final double t) { if (index < 0 || index > 1) { throw new IllegalArgumentException(WRONG_INDEX); } tt[index] = t; } /** Get the stored densities. * @return the densities as an array of 9 values */ public double[] getDensities() { return dd.clone(); } /** 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) { if (index < 0 || index > 8) { throw new IllegalArgumentException(WRONG_INDEX); } return dd[index]; } /** Get the stored temperatures. * @return the temperatures as an array of 2 values */ public double[] getTemperatures() { return tt.clone(); } /** Get one of the stored temperatures. * @param index one of the two elements : * <ul> * <li>{@link #EXOSPHERIC}</li> * <li>{@link #ALTITUDE}</li> * </ul> * @return the requested temperature */ public double getTemperature(final int index) { if (index < 0 || index > 1) { throw new IllegalArgumentException(WRONG_INDEX); } return tt[index]; } } }