NRLMSISE00.java

/* Copyright 2002-2018 CS Systèmes d'Information
 * Licensed to CS Systèmes d'Information (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.forces.drag.atmosphere;

import java.util.Arrays;

import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinatesProvider;


/** This class implements the mathematical representation of the 2001
 *  Naval Research Laboratory Mass Spectrometer and Incoherent Scatter
 *  Radar Exosphere (NRLMSISE-00) of the MSIS® class model.
 *  <p>
 *  NRLMSISE-00 calculates the neutral atmosphere empirical model from the surface
 *  to lower exosphere (0 to 1000 km) and provides:
 *  <ul>
 *  <li>Exospheric Temperature above Input Position (K)</li>
 *  <li>Local Temperature at Input Position (K)</li>
 *  <li>Total Mass-Density at Input Position (kg/m³)</li>
 *  <li>Partial Densities at Input Position (1/m³) for:
 *  <ul>
 *      <li>He,</li>
 *      <li>H,</li>
 *      <li>N,</li>
 *      <li>O,</li>
 *      <li>Ar,</li>
 *      <li>N2,</li>
 *      <li>O2,</li>
 *      <li>anomalous oxygen.</li>
 *  </ul>
 *  </li>
 *  </ul>
 *  </p>
 *  <p>
 *  The model needs geographical and time information to compute general values,
 *  but also needs space weather data:
 *  <ul>
 *  <li>mean and daily solar flux,</li>
 *  <li>geomagnetic indices.</li>
 *  </ul>
 *  </p>
 *  <p>
 *  Switches can be used to turn on and off particular variations:<br>
 *  0 is off, 1 is on, and 2 is main effects off but cross terms on.<br>
 *  The standard value is 1 for all the 23 available switches.<br>
 *  Function of each switch according to its number:
 *  <ul>
 *  <li>#1 - F10.7 effect on mean</li>
 *  <li>#2 - Independent of time</li>
 *  <li>#3 - Symmetrical annual</li>
 *  <li>#4 - Symmetrical semiannual</li>
 *  <li>#5 - Asymmetrical annual</li>
 *  <li>#6 - Asymmetrical semiannual</li>
 *  <li>#7 - Diurnal</li>
 *  <li>#8 - Semidiurnal</li>
 *  <li>#9 - Daily Ap [**]</li>
 *  <li>#10 - All UT, longitudinal effects</li>
 *  <li>#11 - Longitudinal</li>
 *  <li>#12 - UT and mixed UT, longitudinal</li>
 *  <li>#13 - Mixed AP, UT, longitudinal</li>
 *  <li>#14 - Terdiurnal</li>
 *  <li>#15 - Departures from diffusive equilibrium</li>
 *  <li>#16 - All exospheric temperature variations</li>
 *  <li>#17 - All variations from 120 km temperature (TLB)</li>
 *  <li>#18 - All lower thermosphere (TN1) temperature variations</li>
 *  <li>#19 - All 120 km gradient (S) variations</li>
 *  <li>#20 - All upper stratosphere (TN2) temperature variations</li>
 *  <li>#21 - All variations from 120 km values (ZLB)</li>
 *  <li>#22 - All lower mesosphere temperature (TN3) variations</li>
 *  <li>#23 - Turbopause scale height variations</li>
 *  </ul>
 *  [**] Switch #9 is a bit specific:
 *  <ul>
 *  <li>set to  1, the daily Ap only is used (first element of ap array),</li>
 *  <li>set to -1, the entire array of ap is used, including 3 hr ap indices.</li>
 *  </ul>
 *  </p>
 *  <p>
 *  The NRLMSISE-00 model was developed by Mike Picone, Alan Hedin, and Doug Drob.<br>
 *  They also wrote a NRLMSISE-00 distribution package in FORTRAN available at:<br>
 *  ftp://hanna.ccmc.gsfc.nasa.gov/pub/modelweb/atmospheric/msis/nrlmsise00/<br>
 *  <br>
 *  Dominik Brodowski implemented a C version of the NRLMSISE-00 model available at:<br>
 *  http://www.brodo.de/space/nrlmsise/index.html
 *  </p>
 *  <p>
 *  Instances of this class are immutable.
 *  </p>
 *
 *  @author Mike Picone & al (Naval Research Laboratory), 2001: FORTRAN routine
 *  @author Dominik Brodowski, 2004: C routine
 *  @author Pascal Parraud, 2016: Java translation
 *  @since 8.1
 */
public class NRLMSISE00 implements Atmosphere {

    /** Serializable UID. */
    private static final long serialVersionUID = -7923498628122574334L;

    // Constants

    /** Identifier for helium density. */
    private static final int HELIUM = 0;

    /** Identifier for atomic oxygen density. */
    private static final int ATOMIC_OXYGEN = 1;

    /** Identifier for molecular nitrogen density. */
    private static final int MOLECULAR_NITROGEN = 2;

    /** Identifier for molecular oxygen density. */
    private static final int MOLECULAR_OXYGEN = 3;

    /** Identifier for argon density. */
    private static final int ARGON = 4;

    /** Identifier for atomic nitrogen density. */
    private static final int TOTAL_MASS = 5;

    /** Identifier for hydrogen density. */
    private static final int HYDROGEN = 6;

    /** Identifier for atomic nitrogen density. */
    private static final int ATOMIC_NITROGEN = 7;

    /** Identifier for anomalous oxygen density. */
    private static final int ANOMALOUS_OXYGEN = 8;

    /** Identifier for exospheric temperature. */
    private static final int EXOSPHERIC = 0;

    /** Identifier for temperature at altitude. */
    private static final int ALTITUDE = 1;

    // CONVERSION CONSTANTS

    /** Conversion from degree to radian. */
    private static final double DEG_TO_RAD = 1.74533e-2;

    /** Conversion from day to radian. */
    private static final double DAY_TO_RAD = 1.72142e-2;

    /** Conversion from hour to radian. */
    private static final double HOUR_TO_RAD = 0.2618;

    /** Conversion from second to radian. */
    private static final double SEC_TO_RAD = 7.2722e-5;

    // EARTH GEOPHYSICAL CONSTANTS

    /** Reference latitude (°). */
    private static final double LAT_REF = 45.;

    /** Reference gravity on Earth surface at reference latitude (cm/s2). */
    private static final double G_REF = 980.616;

    // CHEMICAL CONSTANTS

    /** Unified atomic mass unit (kg). */
    private static final double AMU = 1.66e-27;

    /** Gas constant (inverse of). */
    private static final double R_GAS = 831.4;

    /** Hydrogen atomic mass. */
    private static final double H_MASS = 1.;

    /** Helium atomic mass. */
    private static final double HE_MASS = 4.;

    /** Nitrogen atomic mass. */
    private static final double N_MASS = 14.;

    /** N2 molecular mass. */
    private static final double N2_MASS = 2. * N_MASS;

    /** Oxygen atomic mass. */
    private static final double O_MASS = 16.;

    /** O2 molecular mass. */
    private static final double O2_MASS = 2. * O_MASS;

    /** Argon atomic mass. */
    private static final double AR_MASS = 40.;

    // NRL MSISE 2000 SPECIFIC CONSTANTS

    /** Reference average flux. */
    private static final double FLUX_REF = 150.;

    /** Array of altitudes #1. */
    private static final double[] ZN1 = {123.435, 110.0, 100.0, 90.0, 72.5};

    /** Array of altitudes #2. */
    private static final double[] ZN2 = {72.5, 55.0, 45.0, 32.5};

    /** Array of altitudes #3. */
    private static final double[] ZN3 = {32.5, 20.0, 15.0, 10.0, 0.0};

    /** Mix altitude (km). */
    private static final double ZMIX = 62.5;

    /** NRLMSISE-00 data: temperature pt[150]. */
    private static final double[] PT = {
        9.86573e-01, 1.62228e-02, 1.55270e-02, -1.04323e-01, -3.75801e-03,
        -1.18538e-03, -1.24043e-01, 4.56820e-03, 8.76018e-03, -1.36235e-01,
        -3.52427e-02, 8.84181e-03, -5.92127e-03, -8.61650e+00, 0.00000e+00,
        1.28492e-02, 0.00000e+00, 1.30096e+02, 1.04567e-02, 1.65686e-03,
        -5.53887e-06, 2.97810e-03, 0.00000e+00, 5.13122e-03, 8.66784e-02,
        1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.27026e-06,
        0.00000e+00, 6.74494e+00, 4.93933e-03, 2.21656e-03, 2.50802e-03,
        0.00000e+00, 0.00000e+00, -2.08841e-02, -1.79873e+00, 1.45103e-03,
        2.81769e-04, -1.44703e-03, -5.16394e-05, 8.47001e-02, 1.70147e-01,
        5.72562e-03, 5.07493e-05, 4.36148e-03, 1.17863e-04, 4.74364e-03,
        6.61278e-03, 4.34292e-05, 1.44373e-03, 2.41470e-05, 2.84426e-03,
        8.56560e-04, 2.04028e-03, 0.00000e+00, -3.15994e+03, -2.46423e-03,
        1.13843e-03, 4.20512e-04, 0.00000e+00, -9.77214e+01, 6.77794e-03,
        5.27499e-03, 1.14936e-03, 0.00000e+00, -6.61311e-03, -1.84255e-02,
        -1.96259e-02, 2.98618e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        6.44574e+02, 8.84668e-04, 5.05066e-04, 0.00000e+00, 4.02881e+03,
        -1.89503e-03, 0.00000e+00, 0.00000e+00, 8.21407e-04, 2.06780e-03,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        -1.20410e-02, -3.63963e-03, 9.92070e-05, -1.15284e-04, -6.33059e-05,
        -6.05545e-01, 8.34218e-03, -9.13036e+01, 3.71042e-04, 0.00000e+00,
        4.19000e-04, 2.70928e-03, 3.31507e-03, -4.44508e-03, -4.96334e-03,
        -1.60449e-03, 3.95119e-03, 2.48924e-03, 5.09815e-04, 4.05302e-03,
        2.24076e-03, 0.00000e+00, 6.84256e-03, 4.66354e-04, 0.00000e+00,
        -3.68328e-04, 0.00000e+00, 0.00000e+00, -1.46870e+02, 0.00000e+00,
        0.00000e+00, 1.09501e-03, 4.65156e-04, 5.62583e-04, 3.21596e+00,
        6.43168e-04, 3.14860e-03, 3.40738e-03, 1.78481e-03, 9.62532e-04,
        5.58171e-04, 3.43731e+00, -2.33195e-01, 5.10289e-04, 0.00000e+00,
        0.00000e+00, -9.25347e+04, 0.00000e+00, -1.99639e-03, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
    };

    /** NRLMSISE-00 data: density pd[9][150]. */
    private static final double[][] PD = {
        // HE DENSITY
        {
            1.09979e+00, -4.88060e-02, -1.97501e-01, -9.10280e-02, -6.96558e-03,
            2.42136e-02, 3.91333e-01, -7.20068e-03, -3.22718e-02, 1.41508e+00,
            1.68194e-01, 1.85282e-02, 1.09384e-01, -7.24282e+00, 0.00000e+00,
            2.96377e-01, -4.97210e-02, 1.04114e+02, -8.61108e-02, -7.29177e-04,
            1.48998e-06, 1.08629e-03, 0.00000e+00, 0.00000e+00, 8.31090e-02,
            1.12818e-01, -5.75005e-02, -1.29919e-02, -1.78849e-02, -2.86343e-06,
            0.00000e+00, -1.51187e+02, -6.65902e-03, 0.00000e+00, -2.02069e-03,
            0.00000e+00, 0.00000e+00, 4.32264e-02, -2.80444e+01, -3.26789e-03,
            2.47461e-03, 0.00000e+00, 0.00000e+00, 9.82100e-02, 1.22714e-01,
            -3.96450e-02, 0.00000e+00, -2.76489e-03, 0.00000e+00, 1.87723e-03,
            -8.09813e-03, 4.34428e-05, -7.70932e-03, 0.00000e+00, -2.28894e-03,
            -5.69070e-03, -5.22193e-03, 6.00692e-03, -7.80434e+03, -3.48336e-03,
            -6.38362e-03, -1.82190e-03, 0.00000e+00, -7.58976e+01, -2.17875e-02,
            -1.72524e-02, -9.06287e-03, 0.00000e+00, 2.44725e-02, 8.66040e-02,
            1.05712e-01, 3.02543e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            -6.01364e+03, -5.64668e-03, -2.54157e-03, 0.00000e+00, 3.15611e+02,
            -5.69158e-03, 0.00000e+00, 0.00000e+00, -4.47216e-03, -4.49523e-03,
            4.64428e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            4.51236e-02, 2.46520e-02, 6.17794e-03, 0.00000e+00, 0.00000e+00,
            -3.62944e-01, -4.80022e-02, -7.57230e+01, -1.99656e-03, 0.00000e+00,
            -5.18780e-03, -1.73990e-02, -9.03485e-03, 7.48465e-03, 1.53267e-02,
            1.06296e-02, 1.18655e-02, 2.55569e-03, 1.69020e-03, 3.51936e-02,
            -1.81242e-02, 0.00000e+00, -1.00529e-01, -5.10574e-03, 0.00000e+00,
            2.10228e-03, 0.00000e+00, 0.00000e+00, -1.73255e+02, 5.07833e-01,
            -2.41408e-01, 8.75414e-03, 2.77527e-03, -8.90353e-05, -5.25148e+00,
            -5.83899e-03, -2.09122e-02, -9.63530e-03, 9.77164e-03, 4.07051e-03,
            2.53555e-04, -5.52875e+00, -3.55993e-01, -2.49231e-03, 0.00000e+00,
            0.00000e+00, 2.86026e+01, 0.00000e+00, 3.42722e-04, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // O DENSITY
        {
            1.02315e+00, -1.59710e-01, -1.06630e-01, -1.77074e-02, -4.42726e-03,
            3.44803e-02, 4.45613e-02, -3.33751e-02, -5.73598e-02, 3.50360e-01,
            6.33053e-02, 2.16221e-02, 5.42577e-02, -5.74193e+00, 0.00000e+00,
            1.90891e-01, -1.39194e-02, 1.01102e+02, 8.16363e-02, 1.33717e-04,
            6.54403e-06, 3.10295e-03, 0.00000e+00, 0.00000e+00, 5.38205e-02,
            1.23910e-01, -1.39831e-02, 0.00000e+00, 0.00000e+00, -3.95915e-06,
            0.00000e+00, -7.14651e-01, -5.01027e-03, 0.00000e+00, -3.24756e-03,
            0.00000e+00, 0.00000e+00, 4.42173e-02, -1.31598e+01, -3.15626e-03,
            1.24574e-03, -1.47626e-03, -1.55461e-03, 6.40682e-02, 1.34898e-01,
            -2.42415e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 6.13666e-04,
            -5.40373e-03, 2.61635e-05, -3.33012e-03, 0.00000e+00, -3.08101e-03,
            -2.42679e-03, -3.36086e-03, 0.00000e+00, -1.18979e+03, -5.04738e-02,
            -2.61547e-03, -1.03132e-03, 1.91583e-04, -8.38132e+01, -1.40517e-02,
            -1.14167e-02, -4.08012e-03, 1.73522e-04, -1.39644e-02, -6.64128e-02,
            -6.85152e-02, -1.34414e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            6.07916e+02, -4.12220e-03, -2.20996e-03, 0.00000e+00, 1.70277e+03,
            -4.63015e-03, 0.00000e+00, 0.00000e+00, -2.25360e-03, -2.96204e-03,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            3.92786e-02, 1.31186e-02, -1.78086e-03, 0.00000e+00, 0.00000e+00,
            -3.90083e-01, -2.84741e-02, -7.78400e+01, -1.02601e-03, 0.00000e+00,
            -7.26485e-04, -5.42181e-03, -5.59305e-03, 1.22825e-02, 1.23868e-02,
            6.68835e-03, -1.03303e-02, -9.51903e-03, 2.70021e-04, -2.57084e-02,
            -1.32430e-02, 0.00000e+00, -3.81000e-02, -3.16810e-03, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -9.05762e-04, -2.14590e-03, -1.17824e-03, 3.66732e+00,
            -3.79729e-04, -6.13966e-03, -5.09082e-03, -1.96332e-03, -3.08280e-03,
            -9.75222e-04, 4.03315e+00, -2.52710e-01, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // N2 DENSITY
        {
            1.16112e+00, 0.00000e+00, 0.00000e+00, 3.33725e-02, 0.00000e+00,
            3.48637e-02, -5.44368e-03, 0.00000e+00, -6.73940e-02, 1.74754e-01,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
            1.26733e-01, 0.00000e+00, 1.03154e+02, 5.52075e-02, 0.00000e+00,
            0.00000e+00, 8.13525e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -2.50482e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.48894e-03,
            6.16053e-04, -5.79716e-04, 2.95482e-03, 8.47001e-02, 1.70147e-01,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // TOTAL MASS
        {
            9.44846e-01, 0.00000e+00, 0.00000e+00, -3.08617e-02, 0.00000e+00,
            -2.44019e-02, 6.48607e-03, 0.00000e+00, 3.08181e-02, 4.59392e-02,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
            2.13260e-02, 0.00000e+00, -3.56958e+02, 0.00000e+00, 1.82278e-04,
            0.00000e+00, 3.07472e-04, 0.00000e+00, 0.00000e+00, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 3.83054e-03, 0.00000e+00, 0.00000e+00,
            -1.93065e-03, -1.45090e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -1.23493e-03, 1.36736e-03, 8.47001e-02, 1.70147e-01,
            3.71469e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            5.10250e-03, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 3.68756e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // O2 DENSITY
        {
            1.35580e+00, 1.44816e-01, 0.00000e+00, 6.07767e-02, 0.00000e+00,
            2.94777e-02, 7.46900e-02, 0.00000e+00, -9.23822e-02, 8.57342e-02,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 2.38636e+01, 0.00000e+00,
            7.71653e-02, 0.00000e+00, 8.18751e+01, 1.87736e-02, 0.00000e+00,
            0.00000e+00, 1.49667e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -3.67874e+02, 5.48158e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
            1.22631e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            8.17187e-03, 3.71617e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.10826e-03,
            -3.13640e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            -7.35742e-02, -5.00266e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 1.94965e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // AR DENSITY
        {
            1.04761e+00, 2.00165e-01, 2.37697e-01, 3.68552e-02, 0.00000e+00,
            3.57202e-02, -2.14075e-01, 0.00000e+00, -1.08018e-01, -3.73981e-01,
            0.00000e+00, 3.10022e-02, -1.16305e-03, -2.07596e+01, 0.00000e+00,
            8.64502e-02, 0.00000e+00, 9.74908e+01, 5.16707e-02, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 3.46193e+02, 1.34297e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.48509e-03,
            -1.54689e-04, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
            1.47753e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            1.89320e-02, 3.68181e-05, 1.32570e-02, 0.00000e+00, 0.00000e+00,
            3.59719e-03, 7.44328e-03, -1.00023e-03, -6.50528e+03, 0.00000e+00,
            1.03485e-02, -1.00983e-03, -4.06916e-03, -6.60864e+01, -1.71533e-02,
            1.10605e-02, 1.20300e-02, -5.20034e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            -2.62769e+03, 7.13755e-03, 4.17999e-03, 0.00000e+00, 1.25910e+04,
            0.00000e+00, 0.00000e+00, 0.00000e+00, -2.23595e-03, 4.60217e-03,
            5.71794e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            -3.18353e-02, -2.35526e-02, -1.36189e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 2.03522e-02, -6.67837e+01, -1.09724e-03, 0.00000e+00,
            -1.38821e-02, 1.60468e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.51574e-02,
            -5.44470e-04, 0.00000e+00, 7.28224e-02, 6.59413e-02, 0.00000e+00,
            -5.15692e-03, 0.00000e+00, 0.00000e+00, -3.70367e+03, 0.00000e+00,
            0.00000e+00, 1.36131e-02, 5.38153e-03, 0.00000e+00, 4.76285e+00,
            -1.75677e-02, 2.26301e-02, 0.00000e+00, 1.76631e-02, 4.77162e-03,
            0.00000e+00, 5.39354e+00, 0.00000e+00, -7.51710e-03, 0.00000e+00,
            0.00000e+00, -8.82736e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // H DENSITY
        {
            1.26376e+00, -2.14304e-01, -1.49984e-01, 2.30404e-01, 2.98237e-02,
            2.68673e-02, 2.96228e-01, 2.21900e-02, -2.07655e-02, 4.52506e-01,
            1.20105e-01, 3.24420e-02, 4.24816e-02, -9.14313e+00, 0.00000e+00,
            2.47178e-02, -2.88229e-02, 8.12805e+01, 5.10380e-02, -5.80611e-03,
            2.51236e-05, -1.24083e-02, 0.00000e+00, 0.00000e+00, 8.66784e-02,
            1.58727e-01, -3.48190e-02, 0.00000e+00, 0.00000e+00, 2.89885e-05,
            0.00000e+00, 1.53595e+02, -1.68604e-02, 0.00000e+00, 1.01015e-02,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.84552e-04,
            -1.22181e-03, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
            -1.04927e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, -5.91313e-03,
            -2.30501e-02, 3.14758e-05, 0.00000e+00, 0.00000e+00, 1.26956e-02,
            8.35489e-03, 3.10513e-04, 0.00000e+00, 3.42119e+03, -2.45017e-03,
            -4.27154e-04, 5.45152e-04, 1.89896e-03, 2.89121e+01, -6.49973e-03,
            -1.93855e-02, -1.48492e-02, 0.00000e+00, -5.10576e-02, 7.87306e-02,
            9.51981e-02, -1.49422e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            2.65503e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 6.37110e-03, 3.24789e-04,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            6.14274e-02, 1.00376e-02, -8.41083e-04, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -1.27099e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            -3.94077e-03, -1.28601e-02, -7.97616e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -6.71465e-03, -1.69799e-03, 1.93772e-03, 3.81140e+00,
            -7.79290e-03, -1.82589e-02, -1.25860e-02, -1.04311e-02, -3.02465e-03,
            2.43063e-03, 3.63237e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // N DENSITY
        {
            7.09557e+01, -3.26740e-01, 0.00000e+00, -5.16829e-01, -1.71664e-03,
            9.09310e-02, -6.71500e-01, -1.47771e-01, -9.27471e-02, -2.30862e-01,
            -1.56410e-01, 1.34455e-02, -1.19717e-01, 2.52151e+00, 0.00000e+00,
            -2.41582e-01, 5.92939e-02, 4.39756e+00, 9.15280e-02, 4.41292e-03,
            0.00000e+00, 8.66807e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,
            1.58727e-01, 9.74701e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 6.70217e+01, -1.31660e-03, 0.00000e+00, -1.65317e-02,
            0.00000e+00, 0.00000e+00, 8.50247e-02, 2.77428e+01, 4.98658e-03,
            6.15115e-03, 9.50156e-03, -2.12723e-02, 8.47001e-02, 1.70147e-01,
            -2.38645e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.37380e-03,
            -8.41918e-03, 2.80145e-05, 7.12383e-03, 0.00000e+00, -1.66209e-02,
            1.03533e-04, -1.68898e-02, 0.00000e+00, 3.64526e+03, 0.00000e+00,
            6.54077e-03, 3.69130e-04, 9.94419e-04, 8.42803e+01, -1.16124e-02,
            -7.74414e-03, -1.68844e-03, 1.42809e-03, -1.92955e-03, 1.17225e-01,
            -2.41512e-02, 1.50521e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            1.60261e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, -3.54403e-04, -1.87270e-02,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            2.76439e-02, 6.43207e-03, -3.54300e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -2.80221e-02, 8.11228e+01, -6.75255e-04, 0.00000e+00,
            -1.05162e-02, -3.48292e-03, -6.97321e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -1.45546e-03, -1.31970e-02, -3.57751e-03, -1.09021e+00,
            -1.50181e-02, -7.12841e-03, -6.64590e-03, -3.52610e-03, -1.87773e-02,
            -2.22432e-03, -3.93895e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        // HOT O DENSITY
        {
            6.04050e-02, 1.57034e+00, 2.99387e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.51018e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, -8.61650e+00, 1.26454e-02,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 5.50878e-03, 0.00000e+00, 0.00000e+00, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 6.23881e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 8.47001e-02, 1.70147e-01,
            -9.45934e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        }
    };

    /** NRLMSISE-00 data: ps[150]. */
    private static final double[] PS = {
        9.56827e-01, 6.20637e-02, 3.18433e-02, 0.00000e+00, 0.00000e+00,
        3.94900e-02, 0.00000e+00, 0.00000e+00, -9.24882e-03, -7.94023e-03,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 1.74712e+02, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 2.74677e-03, 0.00000e+00, 1.54951e-02, 8.66784e-02,
        1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, -6.99007e-04, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 1.24362e-02, -5.28756e-03, 8.47001e-02, 1.70147e-01,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 2.47425e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
        0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
    };

    /** NRLMSISE-00 data: TURBO pdl[2][25]. */
    private static final double[][] PDL = {
        {
            1.09930e+00, 3.90631e+00, 3.07165e+00, 9.86161e-01, 1.63536e+01,
            4.63830e+00, 1.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 1.28840e+00, 3.10302e-02, 1.18339e-01
        },
        {
            1.00000e+00, 7.00000e-01, 1.15020e+00, 3.44689e+00, 1.28840e+00,
            1.00000e+00, 1.08738e+00, 1.22947e+00, 1.10016e+00, 7.34129e-01,
            1.15241e+00, 2.22784e+00, 7.95046e-01, 4.01612e+00, 4.47749e+00,
            1.23435e+02, -7.60535e-02, 1.68986e-06, 7.44294e-01, 1.03604e+00,
            1.72783e+02, 1.15020e+00, 3.44689e+00, -7.46230e-01, 9.49154e-01
        }
    };

    /** NRLMSISE-00 data: LOWER BOUNDARY ptm[10]. */
    private static final double[] PTM = {
        1.04130e+03, 3.86000e+02, 1.95000e+02, 1.66728e+01, 2.13000e+02,
        1.20000e+02, 2.40000e+02, 1.87000e+02, -2.00000e+00, 0.00000e+00
    };

    /** NRLMSISE-00 data: pdm[8][10]. */
    private static final double[][] PDM = {
        {
            2.45600e+07, 6.71072e-06, 1.00000e+02, 0.00000e+00, 1.10000e+02,
            1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        {
            8.59400E+10, 1.00000e+00, 1.05000e+02, -8.00000e+00, 1.10000e+02,
            1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00
        },
        {
            2.81000E+11, 0.00000e+00, 1.05000e+02, 2.80000e+01, 2.89500e+01,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        {
            3.30000E+10, 2.68270e-01, 1.05000e+02, 1.00000e+00, 1.10000e+02,
            1.00000e+01, 1.10000e+02, -1.00000e+01, 0.00000e+00, 0.00000e+00
        },
        {
            1.33000e+09, 1.19615e-02, 1.05000e+02, 0.00000e+00, 1.10000e+02,
            1.00000e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00
        },
        {
            1.76100e+05, 1.00000e+00, 9.50000e+01, -8.00000e+00, 1.10000e+02,
            1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00,
        },
        {
            1.00000e+07, 1.00000e+00, 1.05000e+02, -8.00000e+00, 1.10000e+02,
            1.00000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 0.00000e+00
        },
        {
            1.00000e+06, 1.00000e+00, 1.05000e+02, -8.00000e+00, 5.50000e+02,
            7.60000e+01, 9.00000e+01, 2.00000e+00, 0.00000e+00, 4.00000e+03
        }
    };

    /** NRLMSISE-00 data: ptl[4][100]. */
    private static final double[][] PTL = {
        // TN1(2)
        {
            1.00858e+00, 4.56011e-02, -2.22972e-02, -5.44388e-02, 5.23136e-04,
            -1.88849e-02, 5.23707e-02, -9.43646e-03, 6.31707e-03, -7.80460e-02,
            -4.88430e-02, 0.00000e+00, 0.00000e+00, -7.60250e+00, 0.00000e+00,
            -1.44635e-02, -1.76843e-02, -1.21517e+02, 2.85647e-02, 0.00000e+00,
            0.00000e+00, 6.31792e-04, 0.00000e+00, 5.77197e-03, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -8.90272e+03, 3.30611e-03, 3.02172e-03, 0.00000e+00,
            -2.13673e-03, -3.20910e-04, 0.00000e+00, 0.00000e+00, 2.76034e-03,
            2.82487e-03, -2.97592e-04, -4.21534e-03, 8.47001e-02, 1.70147e-01,
            8.96456e-03, 0.00000e+00, -1.08596e-02, 0.00000e+00, 0.00000e+00,
            5.57917e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 9.65405e-03, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN1(3)
        {
            9.39664e-01, 8.56514e-02, -6.79989e-03, 2.65929e-02, -4.74283e-03,
            1.21855e-02, -2.14905e-02, 6.49651e-03, -2.05477e-02, -4.24952e-02,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 1.19148e+01, 0.00000e+00,
            1.18777e-02, -7.28230e-02, -8.15965e+01, 1.73887e-02, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -1.44691e-02, 2.80259e-04, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 2.16584e+02, 3.18713e-03, 7.37479e-03, 0.00000e+00,
            -2.55018e-03, -3.92806e-03, 0.00000e+00, 0.00000e+00, -2.89757e-03,
            -1.33549e-03, 1.02661e-03, 3.53775e-04, 8.47001e-02, 1.70147e-01,
            -9.17497e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            3.56082e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -1.00902e-02, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN1(4)
        {
            9.85982e-01, -4.55435e-02, 1.21106e-02, 2.04127e-02, -2.40836e-03,
            1.11383e-02, -4.51926e-02, 1.35074e-02, -6.54139e-03, 1.15275e-01,
            1.28247e-01, 0.00000e+00, 0.00000e+00, -5.30705e+00, 0.00000e+00,
            -3.79332e-02, -6.24741e-02, 7.71062e-01, 2.96315e-02, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 6.81051e-03, -4.34767e-03, 8.66784e-02,
            1.58727e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 1.07003e+01, -2.76907e-03, 4.32474e-04, 0.00000e+00,
            1.31497e-03, -6.47517e-04, 0.00000e+00, -2.20621e+01, -1.10804e-03,
            -8.09338e-04, 4.18184e-04, 4.29650e-03, 8.47001e-02, 1.70147e-01,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            -4.04337e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -9.52550e-04,
            8.56253e-04, 4.33114e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.21223e-03,
            2.38694e-04, 9.15245e-04, 1.28385e-03, 8.67668e-04, -5.61425e-06,
            1.04445e+00, 3.41112e+01, 0.00000e+00, -8.40704e-01, -2.39639e+02,
            7.06668e-01, -2.05873e+01, -3.63696e-01, 2.39245e+01, 0.00000e+00,
            -1.06657e-03, -7.67292e-04, 1.54534e-04, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN1(5) TN2(1)
        {
            1.00320e+00, 3.83501e-02, -2.38983e-03, 2.83950e-03, 4.20956e-03,
            5.86619e-04, 2.19054e-02, -1.00946e-02, -3.50259e-03, 4.17392e-02,
            -8.44404e-03, 0.00000e+00, 0.00000e+00, 4.96949e+00, 0.00000e+00,
            -7.06478e-03, -1.46494e-02, 3.13258e+01, -1.86493e-03, 0.00000e+00,
            -1.67499e-02, 0.00000e+00, 0.00000e+00, 5.12686e-04, 8.66784e-02,
            1.58727e-01, -4.64167e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            4.37353e-03, -1.99069e+02, 0.00000e+00, -5.34884e-03, 0.00000e+00,
            1.62458e-03, 2.93016e-03, 2.67926e-03, 5.90449e+02, 0.00000e+00,
            0.00000e+00, -1.17266e-03, -3.58890e-04, 8.47001e-02, 1.70147e-01,
            0.00000e+00, 0.00000e+00, 1.38673e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.60571e-03,
            6.28078e-04, 5.05469e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.57829e-03,
            -4.00855e-04, 5.04077e-05, -1.39001e-03, -2.33406e-03, -4.81197e-04,
            1.46758e+00, 6.20332e+00, 0.00000e+00, 3.66476e-01, -6.19760e+01,
            3.09198e-01, -1.98999e+01, 0.00000e+00, -3.29933e+02, 0.00000e+00,
            -1.10080e-03, -9.39310e-05, 1.39638e-04, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        }
    };

    /** NRLMSISE-00 data: pma[10][100]. */
    private static final double[][] PMA = {
        // TN2(2)
        {
            9.81637e-01, -1.41317e-03, 3.87323e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.58707e-02,
            -8.63658e-03, 0.00000e+00, 0.00000e+00, -2.02226e+00, 0.00000e+00,
            -8.69424e-03, -1.91397e-02, 8.76779e+01, 4.52188e-03, 0.00000e+00,
            2.23760e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -7.07572e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            -4.11210e-03, 3.50060e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -8.36657e-03, 1.61347e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -1.45130e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.24152e-03,
            6.43365e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.33255e-03,
            2.42657e-03, 1.60666e-03, -1.85728e-03, -1.46874e-03, -4.79163e-06,
            1.22464e+00, 3.53510e+01, 0.00000e+00, 4.49223e-01, -4.77466e+01,
            4.70681e-01, 8.41861e+00, -2.88198e-01, 1.67854e+02, 0.00000e+00,
            7.11493e-04, 6.05601e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN2(3)
        {
            1.00422e+00, -7.11212e-03, 5.24480e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -5.28914e-02,
            -2.41301e-02, 0.00000e+00, 0.00000e+00, -2.12219e+01, -1.03830e-02,
            -3.28077e-03, 1.65727e-02, 1.68564e+00, -6.68154e-03, 0.00000e+00,
            1.45155e-02, 0.00000e+00, 8.42365e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -4.34645e-03, 0.00000e+00, 0.00000e+00, 2.16780e-02,
            0.00000e+00, -1.38459e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 7.04573e-03, -4.73204e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 1.08767e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -8.08279e-03,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.21769e-04,
            -2.27387e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.26769e-03,
            3.16901e-03, 4.60316e-04, -1.01431e-04, 1.02131e-03, 9.96601e-04,
            1.25707e+00, 2.50114e+01, 0.00000e+00, 4.24472e-01, -2.77655e+01,
            3.44625e-01, 2.75412e+01, 0.00000e+00, 7.94251e+02, 0.00000e+00,
            2.45835e-03, 1.38871e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN2(4) TN3(1)
        {
            1.01890e+00, -2.46603e-02, 1.00078e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -6.70977e-02,
            -4.02286e-02, 0.00000e+00, 0.00000e+00, -2.29466e+01, -7.47019e-03,
            2.26580e-03, 2.63931e-02, 3.72625e+01, -6.39041e-03, 0.00000e+00,
            9.58383e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -1.85291e-03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 1.39717e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 9.19771e-03, -3.69121e+02, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -1.57067e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.07265e-03,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.92953e-03,
            -2.77739e-03, -4.40092e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.47280e-03,
            2.95035e-04, -1.81246e-03, 2.81945e-03, 4.27296e-03, 9.78863e-04,
            1.40545e+00, -6.19173e+00, 0.00000e+00, 0.00000e+00, -7.93632e+01,
            4.44643e-01, -4.03085e+02, 0.00000e+00, 1.15603e+01, 0.00000e+00,
            2.25068e-03, 8.48557e-04, -2.98493e-04, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN3(2)
        {
            9.75801e-01, 3.80680e-02, -3.05198e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.85575e-02,
            5.04057e-02, 0.00000e+00, 0.00000e+00, -1.76046e+02, 1.44594e-02,
            -1.48297e-03, -3.68560e-03, 3.02185e+01, -3.23338e-03, 0.00000e+00,
            1.53569e-02, 0.00000e+00, -1.15558e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 4.89620e-03, 0.00000e+00, 0.00000e+00, -1.00616e-02,
            -8.21324e-03, -1.57757e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 6.63564e-03, 4.58410e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -2.51280e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 9.91215e-03,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -8.73148e-04,
            -1.29648e-03, -7.32026e-05, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -4.68110e-03,
            -4.66003e-03, -1.31567e-03, -7.39390e-04, 6.32499e-04, -4.65588e-04,
            -1.29785e+00, -1.57139e+02, 0.00000e+00, 2.58350e-01, -3.69453e+01,
            4.10672e-01, 9.78196e+00, -1.52064e-01, -3.85084e+03, 0.00000e+00,
            -8.52706e-04, -1.40945e-03, -7.26786e-04, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN3(3)
        {
            9.60722e-01, 7.03757e-02, -3.00266e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.22671e-02,
            4.10423e-02, 0.00000e+00, 0.00000e+00, -1.63070e+02, 1.06073e-02,
            5.40747e-04, 7.79481e-03, 1.44908e+02, 1.51484e-04, 0.00000e+00,
            1.97547e-02, 0.00000e+00, -1.41844e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 5.77884e-03, 0.00000e+00, 0.00000e+00, 9.74319e-03,
            0.00000e+00, -2.88015e+03, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -4.44902e-03, -2.92760e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 2.34419e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.36685e-03,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -4.65325e-04,
            -5.50628e-04, 3.31465e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.06179e-03,
            -3.08575e-03, -7.93589e-04, -1.08629e-04, 5.95511e-04, -9.05050e-04,
            1.18997e+00, 4.15924e+01, 0.00000e+00, -4.72064e-01, -9.47150e+02,
            3.98723e-01, 1.98304e+01, 0.00000e+00, 3.73219e+03, 0.00000e+00,
            -1.50040e-03, -1.14933e-03, -1.56769e-04, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN3(4)
        {
            1.03123e+00, -7.05124e-02, 8.71615e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -3.82621e-02,
            -9.80975e-03, 0.00000e+00, 0.00000e+00, 2.89286e+01, 9.57341e-03,
            0.00000e+00, 0.00000e+00, 8.66153e+01, 7.91938e-04, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 4.68917e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 7.86638e-03, 0.00000e+00, 0.00000e+00, 9.90827e-03,
            0.00000e+00, 6.55573e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, -4.00200e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 7.07457e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.72268e-03,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.04970e-04,
            1.21560e-03, -8.05579e-06, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.49941e-03,
            -4.57256e-04, -1.59311e-04, 2.96481e-04, -1.77318e-03, -6.37918e-04,
            1.02395e+00, 1.28172e+01, 0.00000e+00, 1.49903e-01, -2.63818e+01,
            0.00000e+00, 4.70628e+01, -2.22139e-01, 4.82292e-02, 0.00000e+00,
            -8.67075e-04, -5.86479e-04, 5.32462e-04, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TN3(5) SURFACE TEMP TSL
        {
            1.00828e+00, -9.10404e-02, -2.26549e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -2.32420e-02,
            -9.08925e-03, 0.00000e+00, 0.00000e+00, 3.36105e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -1.24957e+01, -5.87939e-03, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 2.79765e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 2.01237e+03, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -1.75553e-02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.29699e-03,
            1.26659e-03, 2.68402e-04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 1.17894e-03,
            1.48746e-03, 1.06478e-04, 1.34743e-04, -2.20939e-03, -6.23523e-04,
            6.36539e-01, 1.13621e+01, 0.00000e+00, -3.93777e-01, 2.38687e+03,
            0.00000e+00, 6.61865e+02, -1.21434e-01, 9.27608e+00, 0.00000e+00,
            1.68478e-04, 1.24892e-03, 1.71345e-03, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TGN3(2) SURFACE GRAD TSLG
        {
            1.57293e+00, -6.78400e-01, 6.47500e-01, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -7.62974e-02,
            -3.60423e-01, 0.00000e+00, 0.00000e+00, 1.28358e+02, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 4.68038e+01, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -1.67898e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 2.90994e+04, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 3.15706e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TGN2(1) TGN1(2)
        {
            8.60028e-01, 3.77052e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.17570e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 7.77757e-03, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 1.01024e+02, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 6.54251e+02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, -1.56959e-02,
            1.91001e-02, 3.15971e-02, 1.00982e-02, -6.71565e-03, 2.57693e-03,
            1.38692e+00, 2.82132e-01, 0.00000e+00, 0.00000e+00, 3.81511e+02,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        },
        // TGN3(1) TGN2(2)
        {
            1.06029e+00, -5.25231e-02, 3.73034e-01, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.31072e-02,
            -3.88409e-01, 0.00000e+00, 0.00000e+00, -1.65295e+02, -2.13801e-01,
            -4.38916e-02, -3.22716e-01, -8.82393e+01, 1.18458e-01, 0.00000e+00,
            -4.35863e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, -1.19782e-01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 2.62229e+01, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, -5.37443e+01, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, -4.55788e-01, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 3.84009e-02,
            3.96733e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 5.05494e-02,
            7.39617e-02, 1.92200e-02, -8.46151e-03, -1.34244e-02, 1.96338e-02,
            1.50421e+00, 1.88368e+01, 0.00000e+00, 0.00000e+00, -5.13114e+01,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            5.11923e-02, 3.61225e-02, 0.00000e+00, 0.00000e+00, 0.00000e+00,
            0.00000e+00, 0.00000e+00, 0.00000e+00, 0.00000e+00, 2.00000e+00
        }
    };

    /**  NRLMSISE-00 data: MIDDLE ATMOSPHERE AVERAGES pavgm[10]. */
    private static final double[] PAVGM = {
        2.61000e+02, 2.64000e+02, 2.29000e+02, 2.17000e+02, 2.17000e+02,
        2.23000e+02, 2.86760e+02, -2.93940e+00, 2.50000e+00, 0.00000e+00
    };

    // Fields

    /** External data container. */
    private final NRLMSISE00InputParameters inputParams;

    /** Sun position. */
    private PVCoordinatesProvider sun;

    /** Earth body shape. */
    private final BodyShape earth;

    /** Switches for main effects. */
    private final int[] sw;

    /** Switches for cross effects. */
    private final int[] swc;

    /** Constructor.
     * <p>
     * The model is constructed with all switches set to 1.
     * </p>
     * <p>
     * Parameters are mandatory only for the
     * {@link #getDensity(AbsoluteDate, Vector3D, Frame) getDensity()} and
     * {@link #getVelocity(AbsoluteDate, Vector3D, Frame) getVelocity()} methods.
     * </p>
     * @param parameters the solar and magnetic activity data
     * @param sun the Sun position
     * @param earth the Earth body shape
     */
    public NRLMSISE00(final NRLMSISE00InputParameters parameters,
                      final PVCoordinatesProvider sun,
                      final BodyShape earth) {
        this(parameters, sun, earth, allOnes(), allOnes());
    }

    /** Constructor.
     * <p>
     * The model is constructed with all switches set to 1.
     * </p>
     * <p>
     * Parameters are mandatory only for the
     * {@link #getDensity(AbsoluteDate, Vector3D, Frame) getDensity()} and
     * {@link #getVelocity(AbsoluteDate, Vector3D, Frame) getVelocity()} methods.
     * </p>
     * @param parameters the solar and magnetic activity data
     * @param sun the Sun position
     * @param earth the Earth body shape
     * @param sw switches for main effects
     * @param swc switches for cross effects
     */
    private NRLMSISE00(final NRLMSISE00InputParameters parameters,
                      final PVCoordinatesProvider sun,
                      final BodyShape earth,
                      final int[] sw,
                      final int[] swc) {
        this.inputParams = parameters;
        this.sun         = sun;
        this.earth       = earth;
        this.sw          = sw;
        this.swc         = swc;
    }

    /** Change a switch.
     * <p>
     * This method creates a new instance, the current instance is
     * not changed at all!
     * </p>
     * @param number switch number between 1 and 23
     * @param value switch value
     * @return a <em>new</em> instance, with switch changed
     * @exception OrekitException if switch number is not between 1 and 23
     */
    public NRLMSISE00 withSwitch(final int number, final int value)
        throws OrekitException {
        if (number < 1 || number > 23) {
            throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, number, 1, 23);
        }

        final int[] newSw       = sw.clone();
        final int[] newSwc      = swc.clone();
        if (number != 9) {
            newSw[number]  = (value == 1) ? 1 : 0;
            newSwc[number] = (value > 0) ? 1 : 0;
        } else {
            if (value == -1 || value == 1) {
                newSw[number] = value;
            } else {
                newSw[number] = 0;
            }
            newSwc[number] = newSw[number];
        }

        return new NRLMSISE00(inputParams, sun, earth, newSwc, newSwc);

    }

    /** Create an array of switches set to 1.
     * @return array of switches
     */
    private static int[] allOnes() {
        final int[] array = new int[24];
        Arrays.fill(array, 1);
        return array;
    }

    /** {@inheritDoc} */
    @Override
    public Frame getFrame() {
        return earth.getBodyFrame();
    }

    /** {@inheritDoc} */
    @Override
    public double getDensity(final AbsoluteDate date,
                             final Vector3D position,
                             final Frame frame)
        throws OrekitException {

        // check if data are available :
        if ((date.compareTo(inputParams.getMaxDate()) > 0) ||
            (date.compareTo(inputParams.getMinDate()) < 0)) {
            throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
                                      date, inputParams.getMinDate(), inputParams.getMaxDate());
        }

        // compute day number in current year and the seconds within the day
        final DateTimeComponents dtc = date.getComponents(TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true));
        final int    doy = dtc.getDate().getDayOfYear();
        final double sec = dtc.getTime().getSecondsInLocalDay();

        // compute geodetic position (km and °)
        final GeodeticPoint inBody = earth.transform(position, frame, date);
        final double alt = inBody.getAltitude() / 1000.;
        final double lon = FastMath.toDegrees(inBody.getLongitude());
        final double lat = FastMath.toDegrees(inBody.getLatitude());

        // compute local solar time
        final double lst = localSolarTime(date, position, frame);

        // get solar activity data and compute
        final Output out = new Output(doy, sec, lat, lon, lst, inputParams.getAverageFlux(date),
                                      inputParams.getDailyFlux(date), inputParams.getAp(date));
        out.gtd7d(alt);

        // return the local density
        return out.getDensity(TOTAL_MASS);

    }

    /** {@inheritDoc} */
    @Override
    public <T extends RealFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
                                                        final FieldVector3D<T> position,
                                                        final Frame frame)
        throws OrekitException {
        // check if data are available :
        final AbsoluteDate dateD = date.toAbsoluteDate();
        if ((dateD.compareTo(inputParams.getMaxDate()) > 0) ||
            (dateD.compareTo(inputParams.getMinDate()) < 0)) {
            throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
                                      dateD, inputParams.getMinDate(), inputParams.getMaxDate());
        }

        // compute day number in current year and the seconds within the day
        final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
        final DateTimeComponents dtc = dateD.getComponents(ut1);
        final int    doy = dtc.getDate().getDayOfYear();
        final T sec = date.durationFrom(new AbsoluteDate(dtc.getDate(), TimeComponents.H00, ut1));

        // compute geodetic position (km and °)
        final FieldGeodeticPoint<T> inBody = earth.transform(position, frame, date);
        final T alt = inBody.getAltitude().divide(1000.);
        final T lon = inBody.getLongitude().multiply(180.0 / FastMath.PI);
        final T lat = inBody.getLatitude().multiply(180.0 / FastMath.PI);

        // compute local solar time
        final T lst = localSolarTime(dateD, position, frame);

        // get solar activity data and compute
        final FieldOutput<T> out = new FieldOutput<>(doy, sec, lat, lon, lst,
                                                     inputParams.getAverageFlux(dateD),
                                                     inputParams.getDailyFlux(dateD), inputParams.getAp(dateD));
        out.gtd7d(alt);

        // return the local density
        return out.getDensity(TOTAL_MASS);

    }

    /** Get local solar time.
     * @param date current date
     * @param position current position in frame
     * @param frame the frame in which is defined the position
     * @return the local solar time (hour in [0, 24[)
     * @throws OrekitException if position cannot be computed in given frame
     */
    private double localSolarTime(final AbsoluteDate date,
                                  final Vector3D position,
                                  final Frame frame) throws OrekitException {
        final Vector3D sunPos = sun.getPVCoordinates(date, frame).getPosition();
        final double lst = FastMath.PI + FastMath.atan2(
                sunPos.getX() * position.getY() - sunPos.getY() * position.getX(),
                sunPos.getX() * position.getX() + sunPos.getY() * position.getY());
        return lst * 12. / FastMath.PI;
    }

    /** Get local solar time.
     * @param date current date
     * @param position current position in frame
     * @param frame the frame in which is defined the position
     * @param <T> type of the filed elements
     * @return the local solar time (hour in [0, 24[)
     * @throws OrekitException if position cannot be computed in given frame
     */
    private <T extends RealFieldElement<T>> T localSolarTime(final AbsoluteDate date,
                                                             final FieldVector3D<T> position,
                                                             final Frame frame)
        throws OrekitException {
        final Vector3D sunPos = sun.getPVCoordinates(date, frame).getPosition();
        final T y  = position.getY().multiply(sunPos.getX()).subtract(position.getX().multiply(sunPos.getY()));
        final T x  = position.getX().multiply(sunPos.getX()).add(position.getY().multiply(sunPos.getY()));
        final T hl = y.atan2(x).add(FastMath.PI);

        return hl.multiply(12. / FastMath.PI);

    }

    /**
     * This class is a placeholder for the computed densities and temperatures.
     * <p>
     * Densities are provided as an array d such as:
     * <ul>
     * <li>d[0] = He number density (1/m³)</li>
     * <li>d[1] = O number density (1/m³)</li>
     * <li>d[2] = N2 number density (1/m³)</li>
     * <li>d[3] = O2 number density (1/m³)</li>
     * <li>d[4] = Ar number density (1/m³)</li>
     * <li>d[5] = total mass density (kg/m³) (*)</li>
     * <li>d[6] = H number density (1/m³)</li>
     * <li>d[7] = N number density (1/m³)</li>
     * <li>d[8] = anomalous oxygen number density (1/m³)
     * </ul>
     * Total mass density, d[5], is NOT the same for methods gtd7 and gtd7d:
     * <ul>
     * <li>For gtd7: d[5] is the sum of the mass densities of the species
     * He, O, N2, O2, Ar, H and N but does NOT include anomalous oxygen.</li>
     * <li>For gtd7d: d[5] is the "effective total mass density for drag" and is the sum
     * of the mass densities of all species in this model, INCLUDING anomalous oxygen.</li>
     * </ul>
     * O, H, and N are set to zero below 72.5 km.
     * </p>
     * <p>
     * Temperatures are provided as an array t such as:
     * <ul>
     * <li>t[0] = exospheric temperature (K)</li>
     * <li>t[1] = temperature at altitude (K)</li>
     * </ul>
     * t[0] is set to global average for altitudes below 120 km.<br>
     * The 120 km gradient is left at global average value for altitudes below 72 km.
     * </p>
     */
    private class Output {

        /** Day of year (from 1 to 365 or 366). */
        private final int doy;

        /** Seconds in day (UT scale). */
        private final double sec;

        /** Geodetic latitude (°). */
        private final double lat;

        /** Geodetic longitude (°). */
        private final double lon;

        /** Local apparent solar time (hours). */
        private final double hl;

        /** 81 day average of F10.7 flux (centered on day). */
        private final double f107a;

        /** Daily F10.7 flux for previous day. */
        private final double f107;

        /** Array containing:
        *  <ul>
        *  <li>0: daily Ap</li>
        *  <li>1: 3 hr ap index for current time</li>
        *  <li>2: 3 hr ap index for 3 hrs before current time</li>
        *  <li>3: 3 hr ap index for 6 hrs before current time</li>
        *  <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
        *  <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
        *  <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
        *  </ul>. */
        private final double[] ap;

        /** Gravity at latitude (cm/s2). */
        private final double glat;

        /** Effective Earth radius at latitude (km). */
        private final double rlat;

        /** N2 mixed density at alt. */
        private double dm28;

        /** Legendre polynomials. */
        private final double[][] plg;

        /** Cosinus of local solar time. */
        private final double ctloc;
        /** Sinus of local solar time. */
        private final double stloc;
        /** Square of ctloc. */
        private final double c2tloc;
        /** Square of stloc. */
        private final double s2tloc;
        /** Cube of ctloc. */
        private final double c3tloc;
        /** Cube of stloc. */
        private final double s3tloc;

        /** Magnetic activity based on daily ap. */
        private double apdf;

        /** Magnetic activity based on daily ap. */
        private double apt;

        /** Temperature at nodes for ZN1 scale. */
        private final double[] meso_tn1;

        /** Temperature at nodes for ZN2 scale. */
        private final double[] meso_tn2;

        /** Temperature at nodes for ZN3 scale. */
        private final double[] meso_tn3;

        /** Temperature gradients at end nodes for ZN1 scale. */
        private final double[] meso_tgn1;

        /** Temperature gradients at end nodes for ZN2 scale. */
        private final double[] meso_tgn2;

        /** Temperature gradients at end nodes for ZN3 scale. */
        private final double[] meso_tgn3;

        /** Densities. */
        private final double[] densities;

        /** Temperatures. */
        private final double[] temperatures;

        /** Simple constructor.
         *  @param doy day of year (from 1 to 365 or 366)
         *  @param sec seconds in day (UT scale)
         *  @param lat geodetic latitude (°)
         *  @param lon geodetic longitude (°)
         *  @param hl local apparent solar time (hours)
         *  @param f107a 81 day average of F10.7 flux (centered on day)
         *  @param f107 daily F10.7 flux for previous day
         *  @param ap array containing:
         *  <ul>
         *  <li>0: daily Ap</li>
         *  <li>1: 3 hr ap index for current time</li>
         *  <li>2: 3 hr ap index for 3 hrs before current time</li>
         *  <li>3: 3 hr ap index for 6 hrs before current time</li>
         *  <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
         *  <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
         *  <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
         *  </ul>
         */
        Output(final int doy, final double sec,
               final double lat, final double lon, final double hl,
               final double f107a, final double f107, final double[] ap) {

            this.doy   = doy;
            this.sec   = sec;
            this.lat   = lat;
            this.lon   = lon;
            this.hl    = hl;
            this.f107a = f107a;
            this.f107  = f107;
            this.ap    = ap.clone();

            this.plg       = new double[4][8];

            this.meso_tn1  = new double[ZN1.length];
            this.meso_tn2  = new double[ZN2.length];
            this.meso_tn3  = new double[ZN3.length];
            this.meso_tgn1 = new double[2];
            this.meso_tgn2 = new double[2];
            this.meso_tgn3 = new double[2];

            densities       = new double[9];
            temperatures    = new double[2];

            // Calculates latitude variable gravity and effective radius
            final double xlat = (sw[2] == 0) ? LAT_REF : lat;
            final double c2   = FastMath.cos(2 * DEG_TO_RAD * xlat);
            glat = G_REF * (1. - .0026373 * c2);
            rlat = 2. * glat / (3.085462e-6 + 2.27e-9 * c2) * 1.e-5;

            // Convert latitude into radians
            final double latr = DEG_TO_RAD * lat;

            // Calculate legendre polynomials
            final double c = FastMath.sin(latr);
            final double s = FastMath.cos(latr);

            plg[0][1] = c;
            plg[0][2] = ( 3.0 * c * plg[0][1] - 1.0) / 2.0;
            plg[0][3] = ( 5.0 * c * plg[0][2] - 2.0 * plg[0][1]) / 3.0;
            plg[0][4] = ( 7.0 * c * plg[0][3] - 3.0 * plg[0][2]) / 4.0;
            plg[0][5] = ( 9.0 * c * plg[0][4] - 4.0 * plg[0][3]) / 5.0;
            plg[0][6] = (11.0 * c * plg[0][5] - 5.0 * plg[0][4]) / 6.0;

            plg[1][1] = s;
            plg[1][2] =   3.0 * c * plg[1][1];
            plg[1][3] = ( 5.0 * c * plg[1][2] - 3.0 * plg[1][1]) / 2.0;
            plg[1][4] = ( 7.0 * c * plg[1][3] - 4.0 * plg[1][2]) / 3.0;
            plg[1][5] = ( 9.0 * c * plg[1][4] - 5.0 * plg[1][3]) / 4.0;
            plg[1][6] = (11.0 * c * plg[1][5] - 6.0 * plg[1][4]) / 5.0;

            plg[2][2] = 3.0 * s * plg[1][1];
            plg[2][3] =   5.0 * c * plg[2][2];
            plg[2][4] = ( 7.0 * c * plg[2][3] - 5.0 * plg[2][2]) / 2.0;
            plg[2][5] = ( 9.0 * c * plg[2][4] - 6.0 * plg[2][3]) / 3.0;
            plg[2][6] = (11.0 * c * plg[2][5] - 7.0 * plg[2][4]) / 4.0;
            plg[2][7] = (13.0 * c * plg[2][6] - 8.0 * plg[2][5]) / 5.0;

            plg[3][3] = 5.0 * s * plg[2][2];
            plg[3][4] =   7.0 * c * plg[3][3];
            plg[3][5] = ( 9.0 * c * plg[3][4] - 7.0 * plg[3][3]) / 2.0;
            plg[3][6] = (11.0 * c * plg[3][5] - 8.0 * plg[3][4]) / 3.0;

            // Calculate additional data
            if (!(sw[7] == 0 && sw[8] == 0 && sw[14] == 0)) {
                final double tloc = HOUR_TO_RAD * hl;
                final double tlx2 = tloc + tloc;
                final double tlx3 = tloc + tlx2;
                stloc  = FastMath.sin(tloc);
                ctloc  = FastMath.cos(tloc);
                s2tloc = FastMath.sin(tlx2);
                c2tloc = FastMath.cos(tlx2);
                s3tloc = FastMath.sin(tlx3);
                c3tloc = FastMath.cos(tlx3);
            } else {
                stloc  = 0;
                ctloc  = 0;
                s2tloc = 0;
                c2tloc = 0;
                s3tloc = 0;
                c3tloc = 0;
            }

        }

        /** Calculate temperatures and densities not including anomalous oxygen.
         *  <p>
         *  This method is the thermospheric portion of NRLMSISE-00 for alt > 72.5 km.
         *  </p>
         *  <p>NOTES ON INPUT VARIABLES:<br>
         *  Seconds, Local Time, and Longitude are used independently in the
         *  model and are not of equal importance for every situation.<br>
         *  For the most physically realistic calculation these three
         *  variables should be consistent (lst=sec/3600 + lon/15).<br>
         *  The Equation of Time departures from the above formula
         *  for apparent local time can be included if available but
         *  are of minor importance.<br><br>
         *
         *  f107 and f107A values used to generate the model correspond
         *  to the 10.7 cm radio flux at the actual distance of the Earth
         *  from the Sun rather than the radio flux at 1 AU. The following
         *  site provides both classes of values:<br>
         *  ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
         *
         *  f107, f107A, and ap effects are neither large nor well established below 80 km
         *  and these parameters should be set to 150., 150., and 4. respectively.
         *  </p>
         *  @param alt altitude (km)
         */
        void gts7(final double alt) {

            // Thermal diffusion coefficients for species
            final double[] alpha = {-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
            // Altitude limits for net density computation for species
            final double[] altl  = {200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
            // N2 mixed density
            final double xmm = PDM[2][4];

            /**** Exospheric temperature ****/
            double tinf = PTM[0] * PT[0];
            // Tinf variations not important below ZA or ZN[0]
            if (alt > ZN1[0]) {
                tinf *= 1.0 + sw[16] * globe7(PT);
            }
            setTemperature(EXOSPHERIC, tinf);

            // Gradient variations not important below ZN[4]
            double g0 = PTM[3] * PS[0];
            if (alt > ZN1[4]) {
                g0 *= 1.0 + sw[19] * globe7(PS);
            }

            // Temperature at lower boundary
            double tlb = PTM[1] * PD[3][0];
            tlb *= 1.0 + sw[17] * globe7(PD[3]);

            // Slope
            final double s = g0 / (tinf - tlb);

            // Lower thermosphere temp variations not significant for density above 300 km
            meso_tn1[1]  = PTM[6] * PTL[0][0];
            meso_tn1[2]  = PTM[2] * PTL[1][0];
            meso_tn1[3]  = PTM[7] * PTL[2][0];
            meso_tn1[4]  = PTM[4] * PTL[3][0];
            meso_tgn1[1] = PTM[8] * PMA[8][0];
            if (alt < 300.0) {
                final double r = PTM[4] * PTL[3][0];
                meso_tn1[1]  /= 1.0 - sw[18] * glob7s(PTL[0]);
                meso_tn1[2]  /= 1.0 - sw[18] * glob7s(PTL[1]);
                meso_tn1[3]  /= 1.0 - sw[18] * glob7s(PTL[2]);
                meso_tn1[4]  /= 1.0 - sw[18] * sw[20] * glob7s(PTL[3]);
                meso_tgn1[1] *= 1.0 + sw[18] * sw[20] * glob7s(PMA[8]);
                meso_tgn1[1] *= meso_tn1[4] * meso_tn1[4] / (r * r);
            }

            /**** Temperature at altitude ****/
            setTemperature(ALTITUDE, densu(alt, 1.0, tinf, tlb, 0.0, 0.0, PTM[5], s));

            /**** N2 density ****/
            /*   Density variation factor at Zlb */
            final double g28 = sw[21] * globe7(PD[2]);
            /* Diffusive density at Zlb */
            final double db28 = PDM[2][0] * FastMath.exp(g28) * PD[2][0];
            /* Diffusive density at Alt */
            double diffusiveDensity = densu(alt, db28, tinf, tlb, N2_MASS, alpha[2], PTM[5], s);
            setDensity(MOLECULAR_NITROGEN, diffusiveDensity);
            // Variation of turbopause height
            final double zhf = PDL[1][24] * (1.0 + sw[5] * PDL[0][24] *
                                       FastMath.sin(DEG_TO_RAD * lat) *
                                       FastMath.cos(DAY_TO_RAD * (doy - PT[13])));
            /* Turbopause */
            final double zh28  = PDM[2][2] * zhf;
            final double zhm28 = PDM[2][3] * PDL[1][5];
            /* Mixed density at Zlb */
            final double b28 = densu(zh28, db28, tinf, tlb, N2_MASS - xmm, alpha[2] - 1.0, PTM[5], s);
            if (sw[15] != 0 && alt <= altl[2]) {
                /*  Mixed density at Alt */
                dm28 = densu(alt, b28, tinf, tlb, xmm, alpha[2], PTM[5], s);
                /*  Net density at Alt */
                setDensity(MOLECULAR_NITROGEN, dnet(diffusiveDensity, dm28, zhm28, xmm, N2_MASS));
            }

            /**** He density ****/
            /*   Density variation factor at Zlb */
            final double g4 = sw[21] * globe7(PD[0]);
            /*  Diffusive density at Zlb */
            final double db04 = PDM[0][0] * FastMath.exp(g4) * PD[0][0];
            /*  Diffusive density at Alt */
            diffusiveDensity = densu(alt, db04, tinf, tlb, HE_MASS, alpha[0], PTM[5], s);
            setDensity(HELIUM, diffusiveDensity);
            if (sw[15] != 0 && alt < altl[0]) {
                /*  Turbopause */
                final double zh04 = PDM[0][2];
                /*  Mixed density at Zlb */
                final double b04 = densu(zh04, db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
                /*  Mixed density at Alt */
                final double dm04 = densu(alt, b04, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm04 = zhm28;
                /*  Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm04, zhm04, xmm, HE_MASS);
                /*  Correction to specified mixing ratio at ground */
                final double rl = FastMath.log(b28 * PDM[0][1] / b04);
                final double zc04 = PDM[0][4] * PDL[1][0];
                final double hc04 = PDM[0][5] * PDL[1][1];
                /*  Net density corrected at Alt */
                setDensity(HELIUM, diffusiveDensity * ccor(alt, rl, hc04, zc04));
            }

            /**** O density ****/
            /* Density variation factor at Zlb */
            final double g16 = sw[21] * globe7(PD[1]);
            /* Diffusive density at Zlb */
            final double db16 = PDM[1][0] * FastMath.exp(g16) * PD[1][0];
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db16, tinf, tlb, O_MASS, alpha[1], PTM[5], s);
            setDensity(ATOMIC_OXYGEN, diffusiveDensity);
            if (sw[15] != 0 && alt < altl[1]) {
                /* Turbopause */
                final double zh16 = PDM[1][2];
                /* Mixed density at Zlb */
                final double b16 = densu(zh16, db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
                /* Mixed density at Alt */
                final double dm16 = densu(alt, b16, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm16 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm16, zhm16, xmm, O_MASS);
                final double rl = PDM[1][1] * PDL[1][16] * (1.0 + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
                final double hc16 = PDM[1][5] * PDL[1][3];
                final double zc16 = PDM[1][4] * PDL[1][2];
                final double hc216 = PDM[1][5] * PDL[1][4];
                diffusiveDensity *= ccor2(alt, rl, hc16, zc16, hc216);
                /* Chemistry correction */
                final double hcc16 = PDM[1][7] * PDL[1][13];
                final double zcc16 = PDM[1][6] * PDL[1][12];
                final double rc16  = PDM[1][3] * PDL[1][14];
                /* Net density corrected at Alt */
                setDensity(ATOMIC_OXYGEN, diffusiveDensity * ccor(alt, rc16, hcc16, zcc16));
            }

            /**** O2 density ****/
            /* Density variation factor at Zlb */
            final double g32 = sw[21] * globe7(PD[4]);
            /* Diffusive density at Zlb */
            final double db32 = PDM[3][0] * FastMath.exp(g32) * PD[4][0];
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db32, tinf, tlb, O2_MASS, alpha[3], PTM[5], s);
            setDensity(MOLECULAR_OXYGEN, diffusiveDensity);
            if (sw[15] != 0) {
                if (alt <= altl[3]) {
                    /* Turbopause */
                    final double zh32 = PDM[3][2];
                    /* Mixed density at Zlb */
                    final double b32 = densu(zh32, db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
                    /* Mixed density at Alt */
                    final double dm32 = densu(alt, b32, tinf, tlb, xmm, 0., PTM[5], s);
                    final double zhm32 = zhm28;
                    /* Net density at Alt */
                    diffusiveDensity = dnet(diffusiveDensity, dm32, zhm32, xmm, O2_MASS);
                    /* Correction to specified mixing ratio at ground */
                    final double rl = FastMath.log(b28 * PDM[3][1] / b32);
                    final double hc32 = PDM[3][5] * PDL[1][7];
                    final double zc32 = PDM[3][4] * PDL[1][6];
                    diffusiveDensity *= ccor(alt, rl, hc32, zc32);
                }
                /* Correction for general departure from diffusive equilibrium above Zlb */
                final double hcc32  = PDM[3][7] * PDL[1][22];
                final double hcc232 = PDM[3][7] * PDL[0][22];
                final double zcc32  = PDM[3][6] * PDL[1][21];
                final double rc32   = PDM[3][3] * PDL[1][23] * (1. + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
                /* Net density corrected at Alt */
                setDensity(MOLECULAR_OXYGEN, diffusiveDensity * ccor2(alt, rc32, hcc32, zcc32, hcc232));
            }

            /**** Ar density ****/
            /* Density variation factor at Zlb */
            final double g40 = sw[21] * globe7(PD[5]);
            /* Diffusive density at Zlb */
            final double db40 = PDM[4][0] * FastMath.exp(g40) * PD[5][0];
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db40, tinf, tlb, AR_MASS, alpha[4], PTM[5], s);
            setDensity(ARGON, diffusiveDensity);
            if (sw[15] != 0 && alt <= altl[4]) {
                /* Turbopause */
                final double zh40 = PDM[4][2];
                /* Mixed density at Zlb */
                final double b40 = densu(zh40, db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
                /* Mixed density at Alt */
                final double dm40 = densu(alt, b40, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm40 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm40, zhm40, xmm, AR_MASS);
                /* Correction to specified mixing ratio at ground */
                final double rl = FastMath.log(b28 * PDM[4][1] / b40);
                final double hc40 = PDM[4][5] * PDL[1][9];
                final double zc40 = PDM[4][4] * PDL[1][8];
                /* Net density corrected at Alt */
                setDensity(ARGON, diffusiveDensity * ccor(alt, rl, hc40, zc40));
            }

            /**** H density ****/
            /* Density variation factor at Zlb */
            final double g1 = sw[21] * globe7(PD[6]);
            /* Diffusive density at Zlb */
            final double db01 = PDM[5][0] * FastMath.exp(g1) * PD[6][0];
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db01, tinf, tlb, H_MASS, alpha[6], PTM[5], s);
            setDensity(HYDROGEN, diffusiveDensity);
            if (sw[15] != 0 && alt <= altl[6]) {
                /* Turbopause */
                final double zh01 = PDM[5][2];
                /* Mixed density at Zlb */
                final double b01 = densu(zh01, db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
                /* Mixed density at Alt */
                final double dm01 = densu(alt, b01, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm01 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm01, zhm01, xmm, H_MASS);
                /* Correction to specified mixing ratio at ground */
                final double rl = FastMath.log(b28 * PDM[5][1] * FastMath.sqrt(PDL[1][17] * PDL[1][17]) / b01);
                final double hc01 = PDM[5][5] * PDL[1][11];
                final double zc01 = PDM[5][4] * PDL[1][10];
                diffusiveDensity *= ccor(alt, rl, hc01, zc01);
                /* Chemistry correction */
                final double hcc01 = PDM[5][7] * PDL[1][19];
                final double zcc01 = PDM[5][6] * PDL[1][18];
                final double rc01 = PDM[5][3] * PDL[1][20];
                /* Net density corrected at Alt */
                setDensity(HYDROGEN, diffusiveDensity * ccor(alt, rc01, hcc01, zcc01));
            }

            /**** N density ****/
            /* Density variation factor at Zlb */
            final double g14 = sw[21] * globe7(PD[7]);
            /* Diffusive density at Zlb */
            final double db14 = PDM[6][0] * FastMath.exp(g14) * PD[7][0];
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db14, tinf, tlb, N_MASS, alpha[7], PTM[5], s);
            setDensity(ATOMIC_NITROGEN, diffusiveDensity);
            if (sw[15] != 0 && alt <= altl[7]) {
                /* Turbopause */
                final double zh14 = PDM[6][2];
                /* Mixed density at Zlb */
                final double b14 = densu(zh14, db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
                /* Mixed density at Alt */
                final double dm14 = densu(alt, b14, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm14 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm14, zhm14, xmm, N_MASS);
                /* Correction to specified mixing ratio at ground */
                final double rl = FastMath.log(b28 * PDM[6][1] * PDL[0][2] / b14);
                final double hc14 = PDM[6][5] * PDL[0][1];
                final double zc14 = PDM[6][4] * PDL[0][0];
                diffusiveDensity *= ccor(alt, rl, hc14, zc14);
                /* Chemistry correction */
                final double hcc14 = PDM[6][7] * PDL[0][4];
                final double zcc14 = PDM[6][6] * PDL[0][3];
                final double rc14 = PDM[6][3] * PDL[0][5];
                /* Net density corrected at Alt */
                setDensity(ATOMIC_NITROGEN, diffusiveDensity * ccor(alt, rc14, hcc14, zcc14));
            }

            /**** Anomalous O density ****/
            final double g16h  = sw[21] * globe7(PD[8]);
            final double db16h = PDM[7][0] * FastMath.exp(g16h) * PD[8][0];
            final double tho   = PDM[7][9] * PDL[0][6];
            diffusiveDensity = densu(alt, db16h, tho, tho, O_MASS, alpha[8], PTM[5], s);
            final double zsht = PDM[7][5];
            final double zmho = PDM[7][4];
            final double zsho = scalh(zmho, O_MASS, tho);
            diffusiveDensity *= FastMath.exp(-zsht / zsho * (FastMath.exp((zmho - alt ) / zsht) - 1.));
            setDensity(ANOMALOUS_OXYGEN, diffusiveDensity);

            // Convert densities from cm-3 to m-3
            for (int i = 0; i < 9; i++) {
                setDensity(i, getDensity(i) * 1.0e+06);
            }

            /**** Total mass density ****/
            final double tmd = AMU * (HE_MASS * getDensity(HELIUM) +
                                      O_MASS  * getDensity(ATOMIC_OXYGEN) +
                                      N2_MASS * getDensity(MOLECULAR_NITROGEN) +
                                      O2_MASS * getDensity(MOLECULAR_OXYGEN) +
                                      AR_MASS * getDensity(ARGON) +
                                      H_MASS  * getDensity(HYDROGEN) +
                                      N_MASS  * getDensity(ATOMIC_NITROGEN));
            setDensity(TOTAL_MASS, tmd);

        }

        /** Calculate temperatures and densities not including anomalous oxygen.
         *  <p>NOTES ON INPUT VARIABLES:<br>
         *  Seconds, Local Time, and Longitude are used independently in the
         *  model and are not of equal importance for every situation.<br>
         *  For the most physically realistic calculation these three
         *  variables should be consistent (lst=sec/3600 + lon/15).<br>
         *  The Equation of Time departures from the above formula
         *  for apparent local time can be included if available but
         *  are of minor importance.<br><br>
         *
         *  f107 and f107A values used to generate the model correspond
         *  to the 10.7 cm radio flux at the actual distance of the Earth
         *  from the Sun rather than the radio flux at 1 AU. The following
         *  site provides both classes of values:<br>
         *  ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
         *
         *  f107, f107A, and ap effects are neither large nor well established below 80 km
         *  and these parameters should be set to 150., 150., and 4. respectively.
         *  </p>
         *  @param alt altitude (km)
         */
        void gtd7(final double alt) {

            // Calculates for thermosphere/mesosphere (above ZN2[0])
            final double altt = (alt > ZN2[0]) ? alt : ZN2[0];
            gts7(altt);
            if (alt >= ZN2[0]) {
                return;
            }

            // Calculates for lower mesosphere/upper stratosphere (between ZN2[0] and ZN3[0]):
            // Temperature at nodes and gradients at end nodes
            // Inverse temperature a linear function of spherical harmonics
            final double r = PMA[2][0] * PAVGM[2];
            meso_tgn2[0] = meso_tgn1[1];
            meso_tn2[0]  = meso_tn1[4];
            meso_tn2[1]  = PMA[0][0] * PAVGM[0] / (1.0 - sw[20] * glob7s(PMA[0]));
            meso_tn2[2]  = PMA[1][0] * PAVGM[1] / (1.0 - sw[20] * glob7s(PMA[1]));
            meso_tn2[3]  = PMA[2][0] * PAVGM[2] / (1.0 - sw[20] * sw[22] * glob7s(PMA[2]));
            meso_tgn2[1] = PMA[9][0] * PAVGM[8] * (1.0 + sw[20] * sw[22] * glob7s(PMA[9])) *
                           meso_tn2[3] * meso_tn2[3] / (r * r);
            meso_tn3[0]  = meso_tn2[3];

            // Calculates for lower stratosphere and troposphere (below ZN3[0])
            // Temperature at nodes and gradients at end nodes
            // Inverse temperature a linear function of spherical harmonics
            if (alt < ZN3[0]) {
                final double q = PMA[6][0] * PAVGM[6];
                meso_tgn3[0] = meso_tgn2[1];
                meso_tn3[1]  = PMA[3][0] * PAVGM[3] / (1.0 - sw[22] * glob7s(PMA[3]));
                meso_tn3[2]  = PMA[4][0] * PAVGM[4] / (1.0 - sw[22] * glob7s(PMA[4]));
                meso_tn3[3]  = PMA[5][0] * PAVGM[5] / (1.0 - sw[22] * glob7s(PMA[5]));
                meso_tn3[4]  = PMA[6][0] * PAVGM[6] / (1.0 - sw[22] * glob7s(PMA[6]));
                meso_tgn3[1] = PMA[7][0] * PAVGM[7] * (1.0 + sw[22] * glob7s(PMA[7])) *
                               meso_tn3[4] * meso_tn3[4] / (q * q);

            }

            // Linear transition to full mixing below ZN2[0]
            final double dmc = (alt > ZMIX) ? 1.0 - (ZN2[0] - alt) / (ZN2[0] - ZMIX) : 0.;
            final double dz28 = getDensity(MOLECULAR_NITROGEN);

            // N2 density
            final double dm28m = dm28 * 1.0e+06;
            double dmr = dz28 / dm28m - 1.0;
            double dst = densm(alt, dm28m, PDM[2][4]) * (1.0 + dmr * dmc);
            setDensity(MOLECULAR_NITROGEN, dst);

            // HE density
            dmr = getDensity(HELIUM) / (dz28 * PDM[0][1]) - 1.0;
            dst = getDensity(MOLECULAR_NITROGEN) * PDM[0][1] * (1.0 + dmr * dmc);
            setDensity(HELIUM, dst);

            // O density
            setDensity(ATOMIC_OXYGEN, 0.);
            setDensity(ANOMALOUS_OXYGEN, 0.);

            // O2 density
            dmr = getDensity(MOLECULAR_OXYGEN) / (dz28 * PDM[3][1]) - 1.0;
            dst = getDensity(MOLECULAR_NITROGEN) * PDM[3][1] * (1.0 + dmr * dmc);
            setDensity(MOLECULAR_OXYGEN, dst);

            // AR density
            dmr = getDensity(ARGON) / (dz28 * PDM[4][1]) - 1.0;
            dst = getDensity(MOLECULAR_NITROGEN) * PDM[4][1] * (1.0 + dmr * dmc);
            setDensity(ARGON, dst);

            // H density
            setDensity(HYDROGEN, 0.);

            // N density
            setDensity(ATOMIC_NITROGEN, 0.);

            // Total mass density
            final double tmd = AMU * (HE_MASS * getDensity(HELIUM) +
                                      O_MASS  * getDensity(ATOMIC_OXYGEN) +
                                      N2_MASS * getDensity(MOLECULAR_NITROGEN) +
                                      O2_MASS * getDensity(MOLECULAR_OXYGEN) +
                                      AR_MASS * getDensity(ARGON) +
                                      H_MASS  * getDensity(HYDROGEN) +
                                      N_MASS  * getDensity(ATOMIC_NITROGEN));
            setDensity(TOTAL_MASS, tmd);

            // Temperature at altitude
            setTemperature(ALTITUDE, densm(alt, 1.0, 0));

        }

        /** Calculate temperatures and densities including anomalous oxygen.
         *  <p></p>
         *  <p>NOTES ON INPUT VARIABLES:<br>
         *  Seconds, Local Time, and Longitude are used independently in the
         *  model and are not of equal importance for every situation.<br>
         *  For the most physically realistic calculation these three
         *  variables should be consistent (lst=sec/3600 + lon/15).<br>
         *  The Equation of Time departures from the above formula
         *  for apparent local time can be included if available but
         *  are of minor importance.<br>
         *  <br>
         *  f107 and f107A values used to generate the model correspond
         *  to the 10.7 cm radio flux at the actual distance of the Earth
         *  from the Sun rather than the radio flux at 1 AU. The following
         *  site provides both classes of values:<br>
         *  ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br>
         *  <br>
         *  f107, f107A, and ap effects are neither large nor well established below 80 km
         *  and these parameters should be set to 150., 150., and 4. respectively.
         *  </p>
         *  @param alt altitude (km)
         */
        void gtd7d(final double alt) {

            // Compute densities and temperatures
            gtd7(alt);

            // Update the total mass density with anomalous oxygen contribution
            final double dTot = getDensity(TOTAL_MASS) + AMU * O_MASS * getDensity(ANOMALOUS_OXYGEN);
            setDensity(TOTAL_MASS, dTot);

        }

        /** Set one density.
         * @param index one of the nine elements :
         * <ul>
         * <li>{@link #HELIUM}</li>
         * <li>{@link #ATOMIC_OXYGEN}</li>
         * <li>{@link #MOLECULAR_NITROGEN}</li>
         * <li>{@link #MOLECULAR_OXYGEN}</li>
         * <li>{@link #ARGON}</li>
         * <li>{@link #TOTAL_MASS}</li>
         * <li>{@link #HYDROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * </ul>
         * @param d the value of density to set
         */
        void setDensity(final int index, final double d) {
            densities[index] = d;
        }

        /** Set one temperature.
         * @param index one of the two elements :
         * <ul>
         * <li>{@link #EXOSPHERIC}</li>
         * <li>{@link #ALTITUDE}</li>
         * </ul>
         * @param t the value of temperature to set
         */
        void setTemperature(final int index, final double t) {
            temperatures[index] = t;
        }

        /** Get one of the stored densities.
         * @param index one of the nine elements :
         * <ul>
         * <li>{@link #HELIUM}</li>
         * <li>{@link #ATOMIC_OXYGEN}</li>
         * <li>{@link #MOLECULAR_NITROGEN}</li>
         * <li>{@link #MOLECULAR_OXYGEN}</li>
         * <li>{@link #ARGON}</li>
         * <li>{@link #TOTAL_MASS}</li>
         * <li>{@link #HYDROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * </ul>
         * @return the requested density
         */
        public double getDensity(final int index) {
            return densities[index];
        }

        /** Calculate G(L) function with upper thermosphere parameters.
         *  @param p array of parameters
         *  @return G(L) value
         */
        private double globe7(final double[] p) {

            final double[] t = new double[14];
            final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
            final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
            final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
            final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));

            // F10.7 effect
            final double df  = f107  - f107a;
            final double dfa = f107a - FLUX_REF;
            t[0] = p[19] * df * (1.0 + p[59] * dfa) + p[20] * df * df + p[21] * dfa + p[29] * dfa * dfa;

            final double f1 = 1.0 + (p[47] * dfa + p[19] * df + p[20] * df * df) * swc[1];
            final double f2 = 1.0 + (p[49] * dfa + p[19] * df + p[20] * df * df) * swc[1];

            // Time independent
            t[1] = (p[1]  * plg[0][2] + p[2] * plg[0][4] + p[22] * plg[0][6]) +
                   (p[14] * plg[0][2]) * dfa * swc[1] + p[26] * plg[0][1];

            // Symmetrical annual
            t[2] = p[18] * cd32;

            // Symmetrical semiannual
            t[3] = (p[15] + p[16] * plg[0][2]) * cd18;

            // Asymmetrical annual
            t[4] = f1 * (p[9] * plg[0][1] + p[10] * plg[0][3]) * cd14;

            // Asymmetrical semiannual
            t[5] = p[37] * plg[0][1] * cd39;

            // Diurnal
            if (sw[7] != 0) {
                final double t71 = (p[11] * plg[1][2]) * cd14 * swc[5];
                final double t72 = (p[12] * plg[1][2]) * cd14 * swc[5];
                t[6] = f2 * ((p[3] * plg[1][1] + p[4] * plg[1][3] + p[27] * plg[1][5] + t71) * ctloc +
                             (p[6] * plg[1][1] + p[7] * plg[1][3] + p[28] * plg[1][5] + t72) * stloc);
            }

            // Semidiurnal
            if (sw[8] != 0) {
                final double t81 = (p[23] * plg[2][3] + p[35] * plg[2][5]) * cd14 * swc[5];
                final double t82 = (p[33] * plg[2][3] + p[36] * plg[2][5]) * cd14 * swc[5];
                t[7] = f2 * ((p[5] * plg[2][2] + p[41] * plg[2][4] + t81) * c2tloc +
                             (p[8] * plg[2][2] + p[42] * plg[2][4] + t82) * s2tloc);
            }

            // Terdiurnal
            if (sw[14] != 0) {
                t[13] = f2 * ((p[39] * plg[3][3] + (p[93] * plg[3][4] + p[46] * plg[3][6]) * cd14 * swc[5]) * s3tloc +
                              (p[40] * plg[3][3] + (p[94] * plg[3][4] + p[48] * plg[3][6]) * cd14 * swc[5]) * c3tloc);
            }

            // magnetic activity based on daily ap
            if (sw[9] == -1) {
                if (p[51] != 0) {
                    final double exp1 = FastMath.exp(-10800.0 * FastMath.abs(p[51]) /
                                                     (1.0 + p[138] * (LAT_REF - FastMath.abs(lat))));
                    final double p24 = FastMath.max(p[24], 1.0e-4);
                    apt = sg0(FastMath.min(exp1, 0.99999), p24, p[25]);
                    t[8] = apt * (p[50] + p[96] * plg[0][2] + p[54] * plg[0][4] +
                                  (p[125] * plg[0][1] + p[126] * plg[0][3] + p[127] * plg[0][5]) * cd14 * swc[5] +
                                  (p[128] * plg[1][1] + p[129] * plg[1][3] + p[130] * plg[1][5]) * swc[7] *
                                  FastMath.cos(HOUR_TO_RAD * (hl - p[131])));
                }
            } else {
                final double apd = ap[0] - 4.0;
                final double p44 = (p[43] < 0.) ? 1.0E-5 : p[43];
                final double p45 = p[44];
                apdf = apd + (p45 - 1.0) * (apd + (FastMath.exp(-p44 * apd) - 1.0) / p44);
                if (sw[9] != 0) {
                    t[8] = apdf * (p[32] + p[45] * plg[0][2] + p[34] * plg[0][4] +
                                   (p[100] * plg[0][1] + p[101] * plg[0][3] + p[102] * plg[0][5]) * cd14 * swc[5] +
                                   (p[121] * plg[1][1] + p[122] * plg[1][3] + p[123] * plg[1][5]) * swc[7] *
                                   FastMath.cos(HOUR_TO_RAD * (hl - p[124])));
                }
            }

            if (sw[10] != 0) {
                final double lonr = DEG_TO_RAD * lon;
                // Longitudinal
                if (sw[11] != 0) {
                    t[10] = (1.0 + p[80] * dfa * swc[1]) *
                            ((p[64]  * plg[1][2] + p[65]  * plg[1][4] + p[66]  * plg[1][6] +
                              p[103] * plg[1][1] + p[104] * plg[1][3] + p[105] * plg[1][5] +
                             (p[109] * plg[1][1] + p[110] * plg[1][3] + p[111] * plg[1][5]) * swc[5] * cd14) *
                             FastMath.cos(lonr) +
                             (p[90]  * plg[1][2] + p[91]  * plg[1][4] + p[92]  * plg[1][6] +
                              p[106] * plg[1][1] + p[107] * plg[1][3] + p[108] * plg[1][5] +
                             (p[112] * plg[1][1] + p[113] * plg[1][3] + p[114] * plg[1][5]) * swc[5] * cd14) *
                             FastMath.sin(lonr));
                }

                // ut and mixed ut, longitude
                if (sw[12] != 0) {
                    t[11] = (1.0 + p[95]  * plg[0][1]) * (1.0 + p[81] * dfa * swc[1]) *
                            (1.0 + p[119] * plg[0][1] * swc[5] * cd14) *
                            (p[68] * plg[0][1] + p[69] * plg[0][3] + p[70] * plg[0][5]) *
                            FastMath.cos(SEC_TO_RAD * (sec - p[71]));
                    t[11] += swc[11] * (1.0 + p[137] * dfa * swc[1]) *
                            (p[76] * plg[2][3] + p[77] * plg[2][5] + p[78] * plg[2][7]) *
                            FastMath.cos(SEC_TO_RAD * (sec - p[79]) + 2.0 * lonr);
                }

                /* ut, longitude magnetic activity */
                if (sw[13] != 0) {
                    if (sw[9] == -1) {
                        if (p[51] != 0.) {
                            t[12] = apt * swc[11] * (1. + p[132] * plg[0][1]) *
                                    (p[52] * plg[1][2] + p[98] * plg[1][4] + p[67] * plg[1][6]) *
                                    FastMath.cos(DEG_TO_RAD * (lon - p[97])) +
                                    apt * swc[11] * swc[5] * cd14 *
                                    (p[133] * plg[1][1] + p[134] * plg[1][3] + p[135] * plg[1][5]) *
                                    FastMath.cos(DEG_TO_RAD * (lon - p[136])) +
                                    apt * swc[12] *
                                    (p[55] * plg[0][1] + p[56] * plg[0][3] + p[57] * plg[0][5]) *
                                    FastMath.cos(SEC_TO_RAD * (sec - p[58]));
                        }
                    } else {
                        t[12] = apdf * swc[11] * (1.0 + p[120] * plg[0][1]) *
                                ((p[60] * plg[1][2] + p[61] * plg[1][4] + p[62] * plg[1][6]) *
                                FastMath.cos(DEG_TO_RAD * (lon - p[63]))) +
                                apdf * swc[11] * swc[5] * cd14 *
                                (p[115] * plg[1][1] + p[116] * plg[1][3] + p[117] * plg[1][5]) *
                                FastMath.cos(DEG_TO_RAD * (lon - p[118])) +
                                apdf * swc[12] *
                                (p[83] * plg[0][1] + p[84] * plg[0][3] + p[85] * plg[0][5]) *
                                FastMath.cos(SEC_TO_RAD * (sec - p[75]));
                    }
                }
            }

            // Sum all effects (params not used: 82, 89, 99, 139-149)
            double tinf = p[30];
            for (int i = 0; i < 14; i++) {
                tinf += FastMath.abs(sw[i + 1]) * t[i];
            }

            // Return G(L)
            return tinf;

        }

        /** Calculate G(L) function with lower atmosphere parameters.
         *  @param p array of parameters
         *  @return G(L) value
         */
        private double glob7s(final double[] p) {

            final double[] t = new double[14];
            final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
            final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
            final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
            final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));

            // F10.7 effect
            t[0] = p[21] * (f107a - FLUX_REF);

            // Time independent
            t[1] = p[1]  * plg[0][2] + p[2]  * plg[0][4] + p[22] * plg[0][6] +
                   p[26] * plg[0][1] + p[14] * plg[0][3] + p[59] * plg[0][5];

            // Symmetrical annual
            t[2] = (p[18] + p[47] * plg[0][2] + p[29] * plg[0][4]) * cd32;

            // Symmetrical semiannual
            t[3] = (p[15] + p[16] * plg[0][2] + p[30] * plg[0][4]) * cd18;

            // Asymmetrical annual
            t[4] = (p[9] * plg[0][1] + p[10] * plg[0][3] + p[20] * plg[0][5]) * cd14;

            // Asymmetrical semiannual
            t[5] = (p[37] * plg[0][1]) * cd39;

            // Diurnal
            if (sw[7] != 0) {
                final double t71 = p[11] * plg[1][2] * cd14 * swc[5];
                final double t72 = p[12] * plg[1][2] * cd14 * swc[5];
                t[6] = (p[3] * plg[1][1] + p[4] * plg[1][3] + t71) * ctloc +
                       (p[6] * plg[1][1] + p[7] * plg[1][3] + t72) * stloc;
            }

            // Semidiurnal
            if (sw[8] != 0) {
                final double t81 = (p[23] * plg[2][3] + p[35] * plg[2][5]) * cd14 * swc[5];
                final double t82 = (p[33] * plg[2][3] + p[36] * plg[2][5]) * cd14 * swc[5];
                t[7] = (p[5] * plg[2][2] + p[41] * plg[2][4] + t81) * c2tloc +
                       (p[8] * plg[2][2] + p[42] * plg[2][4] + t82) * s2tloc;
            }

            // Terdiurnal
            if (sw[14] != 0) {
                t[13] = p[39] * plg[3][3] * s3tloc + p[40] * plg[3][3] * c3tloc;
            }

            // Magnetic activity
            if (sw[9] == 1) {
                t[8] = apdf * (p[32] + p[45] * plg[0][2] * swc[2]);
            } else if (sw[9] == -1) {
                t[8] = apt  * (p[50] + p[96] * plg[0][2] * swc[2]);
            }

            // Longitudinal
            if (!(sw[10] == 0 || sw[11] == 0)) {
                final double lonr = DEG_TO_RAD * lon;
                t[10] = (1.0 + plg[0][1] * (p[80] * swc[5] * FastMath.cos(DAY_TO_RAD * (doy - p[81])) +
                                            p[85] * swc[6] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[86]))) +
                               p[83] * swc[3] * FastMath.cos(DAY_TO_RAD * (doy - p[84])) +
                               p[87] * swc[4] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[88]))) *
                        ((p[64] * plg[1][2] + p[65] * plg[1][4] + p[66] * plg[1][6] +
                          p[74] * plg[1][1] + p[75] * plg[1][3] + p[76] * plg[1][5]) * FastMath.cos(lonr) +
                         (p[90] * plg[1][2] + p[91] * plg[1][4] + p[92] * plg[1][6] +
                          p[77] * plg[1][1] + p[78] * plg[1][3] + p[79] * plg[1][5]) * FastMath.sin(lonr));
            }

            // Sum all effects
            double gl = 0;
            for (int i = 0; i < 14; i++) {
                gl += FastMath.abs(sw[i + 1]) * t[i];
            }

            // Return G(L)
            return gl;
        }

        /** Implements sg0 function (Eq. A24a).
         * @param ex ex
         * @param p24 abs(p[24])
         * @param p25 p[25]
         * @return sg0
         */
        private double sg0(final double ex, final double p24, final double p25) {
            final double g01 = g0(ap[1], p24, p25);
            final double g02 = g0(ap[2], p24, p25);
            final double g03 = g0(ap[3], p24, p25);
            final double g04 = g0(ap[4], p24, p25);
            final double g05 = g0(ap[5], p24, p25);
            final double g06 = g0(ap[6], p24, p25);
            final double ex2 = ex * ex;
            final double ex3 = ex * ex2;
            final double ex4 = ex2 * ex2;
            final double ex8 = ex4 * ex4;
            final double ex12 = ex4 * ex8;
            final double g234 = g02 * ex + g03 * ex2 + g04 * ex3;
            final double g56  = g05 * ex4 + g06 * ex12;
            final double ex19 = ex3 * ex4 * ex12;
            final double omex = 1.0 - ex;
            final double sumex = 1.0 + (1.0 - ex19) / omex * FastMath.sqrt(ex);
            return (g01 + (g234 + g56 * (1.0 - ex8) / omex)) / sumex;
        }

        /** Implements go function (Eq. A24d).
         * @param apI 3 hrs ap
         * @param p24 abs(p[24])
         * @param p25 p[25]
         * @return go
         */
        private double g0(final double apI, final double p24, final double p25) {
            final double am4 = apI - 4.0;
            return am4 + (p25 - 1.0) * (am4 + (FastMath.exp(-p24 * am4) - 1.0) / p24);
        }

        /** Calculates chemistry/dissociation correction for MSIS models.
         * @param alt altitude
         * @param r target ratio
         * @param h1 transition scale length
         * @param zh altitude of 1/2 R
         * @return correction
         */
        private double ccor(final double alt, final double r, final double h1, final double zh) {
            final double e = (alt - zh) / h1;
            if (e > 70.) {
                return 1.;
            } else if (e < -70.) {
                return FastMath.exp(r);
            } else {
                return FastMath.exp(r / (1.0 + FastMath.exp(e)));
            }
        }


        /** Calculates O & O2 chemistry/dissociation correction for MSIS models.
         * @param alt altitude
         * @param r target ratio
         * @param h1 transition scale length
         * @param zh altitude of 1/2 R
         * @param h2 transition scale length
         * @return correction
         */
        private double ccor2(final double alt, final double r,
                             final double h1, final double zh, final double h2) {
            final double e1 = (alt - zh) / h1;
            final double e2 = (alt - zh) / h2;
            if ((e1 > 70.) || (e2 > 70.)) {
                return 1.;
            } else if ((e1 < -70.) && (e2 < -70.)) {
                return FastMath.exp(r);
            } else {
                final double ex1 = FastMath.exp(e1);
                final double ex2 = FastMath.exp(e2);
                return FastMath.exp(r / (1.0 + 0.5 * (ex1 + ex2)));
            }
        }

        /** Calculates scale height.
         * @param alt altitude
         * @param xm species molecular weight
         * @param temp temperature
         * @return scale height (km)
         */
        private double scalh(final double alt, final double xm, final double temp) {
            // Gravity at altitude
            final double denom = 1.0 + alt / rlat;
            final double galt = glat / (denom * denom);
            return R_GAS * temp / (galt * xm);
        }

        /** Calculates turbopause correction for MSIS models.
         * @param dd diffusive density
         * @param dm full mixed density
         * @param zhm transition scale length
         * @param xmm full mixed molecular weight
         * @param xm species molecular weight
         * @return combined density
         */
        private double dnet(final double dd, final double dm,
                            final double zhm, final double xmm, final double xm) {
            if (!(dm > 0 && dd > 0)) {
                double ddd = dd;
                if (dd == 0 && dm == 0) {
                    ddd = 1;
                }
                if (dm == 0) {
                    return ddd;
                }
                if (dd == 0) {
                    return dm;
                }
            }

            final double a  = zhm / (xmm - xm);
            final double ylog = a * FastMath.log(dm / dd);
            if (ylog < -10.) {
                return dd;
            } else if (ylog > 10.) {
                return dm;
            } else {
                return dd * FastMath.pow(1.0 + FastMath.exp(ylog), 1.0 / a);
            }
        }

        /** Integrate cubic spline function from xa[0] to x.
         * <p>ADAPTED FROM NUMERICAL RECIPES</p>
         * @param xa array of abscissas in ascending order
         * @param ya array of ordinates in ascending order by xa
         * @param y2a array of second derivatives in ascending order by xa
         * @param x abscissa end point
         * @return integral value
         */
        private double splini(final double[] xa, final double[] ya, final double[] y2a, final double x) {
            final int n = xa.length;
            double yi = 0;
            int klo = 0;
            int khi = 1;
            while (x > xa[klo] && khi < n) {
                double xx = x;
                if (khi < n - 1) {
                    xx = (x < xa[khi]) ? x : xa[khi];
                }
                final double h = xa[khi] - xa[klo];
                final double a = (xa[khi] - xx) / h;
                final double b = (xx - xa[klo]) / h;
                final double a2 = a * a;
                final double b2 = b * b;
                yi += ((1.0 - a2) * ya[klo] / 2.0 + b2 * ya[khi] / 2.0 +
                       ((-(1.0 + a2 * a2) / 4.0 + a2 / 2.0) * y2a[klo] +
                          (b2 * b2 / 4.0 - b2 / 2.0) * y2a[khi]) * h * h / 6.0) * h;
                klo++;
                khi++;
            }
            return yi;
        }

        /** Calculate cubic spline interpolated value.
         * <p>ADAPTED FROM NUMERICAL RECIPES</p>
         * @param xa array of abscissas in ascending order
         * @param ya array of ordinates in ascending order by xa
         * @param y2a array of second derivatives in ascending order by xa
         * @param x abscissa for interpolation
         * @return interpolated value
         */
        private double splint(final double[] xa, final double[] ya, final double[] y2a, final double x) {
            final int n = xa.length;
            int klo = 0;
            int khi = n - 1;
            while (khi - klo > 1) {
                final int k = (khi + klo) >>> 1;
                if (xa[k] > x) {
                    khi = k;
                } else {
                    klo = k;
                }
            }
            final double h = xa[khi] - xa[klo];
            final double a = (xa[khi] - x) / h;
            final double b = (x - xa[klo]) / h;
            return a * ya[klo] + b * ya[khi] +
                    ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * h * h / 6.0;
        }

        /** Calculate 2nd derivatives of cubic spline interpolation function.
         * <p>ADAPTED FROM NUMERICAL RECIPES</p>
         * @param x array of abscissas in ascending order
         * @param y array of ordinates in ascending order by x
         * @param yp1 derivative at x[0] (2nd derivatives null if > 1E30)
         * @param ypn derivative at x[n-1] (2nd derivatives null if > 1E30)
         * @return array of second derivatives
         */
        private double[] spline(final double[] x, final double[] y, final double yp1, final double ypn) {
            final int n = x.length;
            final double[] y2 = new double[n];
            final double[] u  = new double[n];

            if (yp1 < 1e+30) {
                y2[0] = -0.5;
                u[0]  = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
            }
            for (int i = 1; i < n - 1; i++) {
                final double sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
                final double p = sig * y2[i - 1] + 2.0;
                y2[i] = (sig - 1.0) / p;
                u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) /
                        (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
            }

            double qn = 0;
            double un = 0;
            if (ypn < 1e+30) {
                qn = 0.5;
                un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
            }

            y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
            for (int k = n - 2; k >= 0; k--) {
                y2[k] = y2[k] * y2[k + 1] + u[k];
            }

            return y2;
        }

        /** Calculate Temperature and Density Profiles for lower atmosphere.
         * @param alt altitude
         * @param d0 density
         * @param xm mixed density
         * @return temperature or density profile
         */
        private double densm(final double alt, final double d0, final double xm) {

            double densm = d0;

            // stratosphere/mesosphere temperature
            int mn = ZN2.length;
            double z = (alt > ZN2[mn - 1]) ? alt : ZN2[mn - 1];

            double z1 = ZN2[0];
            double z2 = ZN2[mn - 1];
            double t1 = meso_tn2[0];
            double t2 = meso_tn2[mn - 1];
            double zg  = zeta(z, z1);
            double zgdif = zeta(z2, z1);

            /* set up spline nodes */
            double[] xs = new double[mn];
            double[] ys = new double[mn];
            for (int k = 0; k < mn; k++) {
                xs[k] = zeta(ZN2[k], z1) / zgdif;
                ys[k] = 1.0 / meso_tn2[k];
            }
            final double qSM = (rlat + z2) / (rlat + z1);
            double yd1 = -meso_tgn2[0] / (t1 * t1) * zgdif;
            double yd2 = -meso_tgn2[1] / (t2 * t2) * zgdif * qSM * qSM;

            /* calculate spline coefficients */
            double[] y2out = spline(xs, ys, yd1, yd2);
            double x = zg / zgdif;
            double y = splint(xs, ys, y2out, x);

            /* temperature at altitude */
            double tz = 1.0 / y;

            if (xm != 0.0) {
                /* calculate stratosphere / mesospehere density */
                final double glb  = galt(z1);
                final double gamm = xm * glb * zgdif / R_GAS;

                /* Integrate temperature profile */
                final double yi = splini(xs, ys, y2out, x);
                final double expl = FastMath.min(50., gamm * yi);

                /* Density at altitude */
                densm *= (t1 / tz) * FastMath.exp(-expl);
            }

            if (alt > ZN3[0]) {
                return (xm == 0.0) ? tz : densm;
            }

            // troposhere/stratosphere temperature
            z = alt;
            mn = ZN3.length;
            z1 = ZN3[0];
            z2 = ZN3[mn - 1];
            t1 = meso_tn3[0];
            t2 = meso_tn3[mn - 1];
            zg = zeta(z, z1);
            zgdif = zeta(z2, z1);

            /* set up spline nodes */
            xs = new double[mn];
            ys = new double[mn];
            for (int k = 0; k < mn; k++) {
                xs[k] = zeta(ZN3[k], z1) / zgdif;
                ys[k] = 1.0 / meso_tn3[k];
            }
            final double qTS = (rlat + z2) / (rlat + z1);
            yd1 = -meso_tgn3[0] / (t1 * t1) * zgdif;
            yd2 = -meso_tgn3[1] / (t2 * t2) * zgdif * qTS * qTS;

            /* calculate spline coefficients */
            y2out = spline(xs, ys, yd1, yd2);
            x = zg / zgdif;
            y = splint(xs, ys, y2out, x);

            /* temperature at altitude */
            tz = 1.0 / y;

            if (xm != 0.0) {
                /* calculate tropospheric / stratosphere density */
                final double glb = galt(z1);
                final double gamm = xm * glb * zgdif / R_GAS;

                /* Integrate temperature profile */
                final double yi = splini(xs, ys, y2out, x);
                final double expl = FastMath.min(50., gamm * yi);

                /* Density at altitude */
                densm *= (t1 / tz) * FastMath.exp(-expl);
            }

            return (xm == 0.0) ? tz : densm;
        }

        /** Calculate temperature and density profiles according to new lower thermo polynomial.
         * @param alt altitude
         * @param dlb density at lower boundary
         * @param tinf exospheric temperature
         * @param tlb temperature at lower boundary
         * @param xm species molecular weight
         * @param alpha thermal diffusion coefficient
         * @param zlb altitude of the lower boundary
         * @param s2 slope
         * @return temperature or density profile
         */
        private double densu(final double alt, final double dlb, final double tinf,
                             final double tlb, final double xm, final double alpha,
                             final double zlb, final double s2) {
            /* joining altitudes of Bates and spline */
            double z = (alt > ZN1[0]) ? alt : ZN1[0];

            /* geopotential altitude difference from ZLB */
            final double zg2 = zeta(z, zlb);

            /* Bates temperature */
            final double tt = tinf - (tinf - tlb) * FastMath.exp(-s2 * zg2);
            final double ta = tt;
            double tz = tt;

            final int mn = ZN1.length;
            final double[] xs = new double[mn];
            final double[] ys = new double[mn];
            double x = 0.;
            double[] y2out =  new double[mn];
            double zgdif = 0.;
            if (alt < ZN1[0]) {
                /* calculate temperature below ZA
                 * temperature gradient at ZA from Bates profile */
                final double p = (rlat + zlb) / (rlat + ZN1[0]);
                final double dta = (tinf - ta) * s2 * p * p;
                meso_tgn1[0] = dta;
                meso_tn1[0] = ta;
                z = (alt > ZN1[mn - 1]) ? alt : ZN1[mn - 1];

                final double t1 = meso_tn1[0];
                final double t2 = meso_tn1[mn - 1];
                /* geopotental difference from z1 */
                final double zg = zeta(z, ZN1[0]);
                zgdif = zeta(ZN1[mn - 1], ZN1[0]);
                /* set up spline nodes */
                for (int k = 0; k < mn; k++) {
                    xs[k] = zeta(ZN1[k], ZN1[0]) / zgdif;
                    ys[k] = 1.0 / meso_tn1[k];
                }
                /* end node derivatives */
                final double q   = (rlat + ZN1[mn - 1]) / (rlat + ZN1[0]);
                final double yd1 = -meso_tgn1[0] / (t1 * t1) * zgdif;
                final double yd2 = -meso_tgn1[1] / (t2 * t2) * zgdif * q * q;
                /* calculate spline coefficients */
                y2out = spline(xs, ys, yd1, yd2);
                x = zg / zgdif;
                final double y = splint(xs, ys, y2out, x);
                /* temperature at altitude */
                tz = 1.0 / y;
            }

            if (xm == 0) {
                return tz;
            }

            /* calculate density above za */
            double glb   = galt(zlb);
            double gamma = xm * glb / (R_GAS * s2 * tinf);
            double expl  = (tt <= 0) ? 50. : FastMath.min(50., FastMath.exp(-s2 * gamma * zg2));
            double densu = dlb * FastMath.pow(tlb / tt, 1.0 + alpha + gamma) * expl;

            /* calculate density below za */
            if (alt < ZN1[0]) {
                glb   = galt(ZN1[0]);
                gamma = xm * glb * zgdif / R_GAS;
                /* integrate spline temperatures */
                expl  = (tz <= 0) ? 50.0 : FastMath.min(50., gamma * splini(xs, ys, y2out, x));
                /* correct density at altitude */
                densu *= FastMath.pow(meso_tn1[0] / tz, 1.0 + alpha) * FastMath.exp(-expl);
            }

            /* Return density at altitude */
            return densu;
        }

        /** Calculate gravity at altitude.
         * @param alt altitude (km)
         * @return gravity at altitude (cm/s2)
         */
        private double galt(final double alt) {
            final double r = 1.0 + alt / rlat;
            return glat / (r * r);
        }

        /** Calculate zeta function.
         * @param zz zz value
         * @param zl zl value
         * @return value of zeta function
         */
        private double zeta(final double zz, final double zl) {
            return (zz - zl) * (rlat + zl) / (rlat + zz);
        }

    }

    /**
     * This class is a placeholder for the computed densities and temperatures.
     * <p>
     * Densities are provided as an array d such as:
     * <ul>
     * <li>d[0] = He number density (1/m³)</li>
     * <li>d[1] = O number density (1/m³)</li>
     * <li>d[2] = N2 number density (1/m³)</li>
     * <li>d[3] = O2 number density (1/m³)</li>
     * <li>d[4] = Ar number density (1/m³)</li>
     * <li>d[5] = total mass density (kg/m³) (*)</li>
     * <li>d[6] = H number density (1/m³)</li>
     * <li>d[7] = N number density (1/m³)</li>
     * <li>d[8] = anomalous oxygen number density (1/m³)
     * </ul>
     * Total mass density, d[5], is NOT the same for methods gtd7 and gtd7d:
     * <ul>
     * <li>For gtd7: d[5] is the sum of the mass densities of the species
     * He, O, N2, O2, Ar, H and N but does NOT include anomalous oxygen.</li>
     * <li>For gtd7d: d[5] is the "effective total mass density for drag" and is the sum
     * of the mass densities of all species in this model, INCLUDING anomalous oxygen.</li>
     * </ul>
     * O, H, and N are set to zero below 72.5 km.
     * </p>
     * <p>
     * Temperatures are provided as an array t such as:
     * <ul>
     * <li>t[0] = exospheric temperature (K)</li>
     * <li>t[1] = temperature at altitude (K)</li>
     * </ul>
     * t[0] is set to global average for altitudes below 120 km.<br>
     * The 120 km gradient is left at global average value for altitudes below 72 km.
     * </p>
     * @param <T> type of the field elements
     * @since 9.0
     */
    public class FieldOutput<T extends RealFieldElement<T>> {

        /** Type of the field elements. */
        private final Field<T> field;

        /** Zero for the field. */
        private final T zero;

        /** Day of year (from 1 to 365 or 366). */
        private final int doy;

        /** Seconds in day (UT scale). */
        private final T sec;

        /** Geodetic latitude (°). */
        private final T lat;

        /** Geodetic longitude (°). */
        private final T lon;

        /** Local apparent solar time (hours). */
        private final T hl;

        /** 81 day average of F10.7 flux (centered on day). */
        private final double f107a;

        /** Daily F10.7 flux for previous day. */
        private final double f107;

        /** Array containing:
        *  <ul>
        *  <li>0: daily Ap</li>
        *  <li>1: 3 hr ap index for current time</li>
        *  <li>2: 3 hr ap index for 3 hrs before current time</li>
        *  <li>3: 3 hr ap index for 6 hrs before current time</li>
        *  <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
        *  <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
        *  <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
        *  </ul>. */
        private final double[] ap;

        /** Gravity at latitude (cm/s2). */
        private final T glat;

        /** Effective Earth radius at latitude (km). */
        private final T rlat;

        /** N2 mixed density at alt. */
        private T dm28;

        /** Legendre polynomials. */
        private final T[][] plg;

        /** Cosinus of local solar time. */
        private final T ctloc;
        /** Sinus of local solar time. */
        private final T stloc;
        /** Square of ctloc. */
        private final T c2tloc;
        /** Square of stloc. */
        private final T s2tloc;
        /** Cube of ctloc. */
        private final T c3tloc;
        /** Cube of stloc. */
        private final T s3tloc;

        /** Magnetic activity based on daily ap. */
        private double apdf;

        /** Magnetic activity based on daily ap. */
        private T apt;

        /** Temperature at nodes for ZN1 scale. */
        private final T[] meso_tn1;

        /** Temperature at nodes for ZN2 scale. */
        private final T[] meso_tn2;

        /** Temperature at nodes for ZN3 scale. */
        private final T[] meso_tn3;

        /** Temperature gradients at end nodes for ZN1 scale. */
        private final T[] meso_tgn1;

        /** Temperature gradients at end nodes for ZN2 scale. */
        private final T[] meso_tgn2;

        /** Temperature gradients at end nodes for ZN3 scale. */
        private final T[] meso_tgn3;

        /** Densities. */
        private final T[] densities;

        /** Temperatures. */
        private final T[] temperatures;

        /** Simple constructor.
         *  @param doy day of year (from 1 to 365 or 366)
         *  @param sec seconds in day (UT scale)
         *  @param lat geodetic latitude (°)
         *  @param lon geodetic longitude (°)
         *  @param hl local apparent solar time (hours)
         *  @param f107a 81 day average of F10.7 flux (centered on day)
         *  @param f107 daily F10.7 flux for previous day
         *  @param ap array containing:
         *  <ul>
         *  <li>0: daily Ap</li>
         *  <li>1: 3 hr ap index for current time</li>
         *  <li>2: 3 hr ap index for 3 hrs before current time</li>
         *  <li>3: 3 hr ap index for 6 hrs before current time</li>
         *  <li>4: 3 hr ap index for FOR 9 hrs before current time</li>
         *  <li>5: average of eight 3 hr ap indices from 12 to 33 hrs prior to current time</li>
         *  <li>6: average of eight 3 hr ap indices from 36 to 57 hrs prior to current time</li>
         *  </ul>
         */
        FieldOutput(final int doy, final T sec,
                    final T lat, final T lon, final T hl,
                    final double f107a, final double f107, final double[] ap) {

            this.field = sec.getField();
            this.zero = field.getZero();

            this.doy   = doy;
            this.sec   = sec;
            this.lat   = lat;
            this.lon   = lon;
            this.hl    = hl;
            this.f107a = f107a;
            this.f107  = f107;
            this.ap    = ap.clone();

            this.plg       = MathArrays.buildArray(field, 4, 8);

            this.meso_tn1  = MathArrays.buildArray(field, ZN1.length);
            this.meso_tn2  = MathArrays.buildArray(field, ZN2.length);
            this.meso_tn3  = MathArrays.buildArray(field, ZN3.length);
            this.meso_tgn1 = MathArrays.buildArray(field, 2);
            this.meso_tgn2 = MathArrays.buildArray(field, 2);
            this.meso_tgn3 = MathArrays.buildArray(field, 2);

            densities       = MathArrays.buildArray(field, 9);
            temperatures    = MathArrays.buildArray(field, 2);

            // Calculates latitude variable gravity and effective radius
            final T xlat = (sw[2] == 0) ? zero.add(LAT_REF) : lat;
            final T c2   = xlat.multiply(2 * DEG_TO_RAD).cos();
            glat = c2.multiply(-0.0026373).add(1).multiply(G_REF);
            rlat = glat.multiply(2).divide(c2.multiply(2.27e-9).add(3.085462e-6)).multiply(1.e-5);

            // Convert latitude into radians
            final T latr = lat.multiply(DEG_TO_RAD);

            // Calculate legendre polynomials
            final T c = latr.sin();
            final T s = latr.cos();

            plg[0][1] = c;
            plg[0][2] = c.multiply( 3.0).multiply(plg[0][1]).subtract(1.0).divide(2.0);
            plg[0][3] = c.multiply( 5.0).multiply(plg[0][2]).subtract(plg[0][1].multiply(2.0)).divide(3.0);
            plg[0][4] = c.multiply( 7.0).multiply(plg[0][3]).subtract(plg[0][2].multiply(3.0)).divide(4.0);
            plg[0][5] = c.multiply( 9.0).multiply(plg[0][4]).subtract(plg[0][3].multiply(4.0)).divide(5.0);
            plg[0][6] = c.multiply(11.0).multiply(plg[0][5]).subtract(plg[0][4].multiply(5.0)).divide(6.0);

            plg[1][1] = s;
            plg[1][2] = c.multiply( 3.0).multiply(plg[1][1]);
            plg[1][3] = c.multiply( 5.0).multiply(plg[1][2]).subtract(plg[1][1].multiply(3.0)).divide(2.0);
            plg[1][4] = c.multiply( 7.0).multiply(plg[1][3]).subtract(plg[1][2].multiply(4.0)).divide(3.0);
            plg[1][5] = c.multiply( 9.0).multiply(plg[1][4]).subtract(plg[1][3].multiply(5.0)).divide(4.0);
            plg[1][6] = c.multiply(11.0).multiply(plg[1][5]).subtract(plg[1][4].multiply(6.0)).divide(5.0);

            plg[2][2] = s.multiply( 3.0).multiply(plg[1][1]);
            plg[2][3] = c.multiply( 5.0).multiply(plg[2][2]);
            plg[2][4] = c.multiply( 7.0).multiply(plg[2][3]).subtract(plg[2][2].multiply(5.0)).divide(2.0);
            plg[2][5] = c.multiply( 9.0).multiply(plg[2][4]).subtract(plg[2][3].multiply(6.0)).divide(3.0);
            plg[2][6] = c.multiply(11.0).multiply(plg[2][5]).subtract(plg[2][4].multiply(7.0)).divide(4.0);
            plg[2][7] = c.multiply(13.0).multiply(plg[2][6]).subtract(plg[2][5].multiply(8.0)).divide(5.0);

            plg[3][3] = s.multiply( 5.0).multiply(plg[2][2]);
            plg[3][4] = c.multiply( 7.0).multiply(plg[3][3]);
            plg[3][5] = c.multiply( 9.0).multiply(plg[3][4]).subtract(plg[3][3].multiply(7.0)).divide(2.0);
            plg[3][6] = c.multiply(11.0).multiply(plg[3][5]).subtract(plg[3][4].multiply(8.0)).divide(3.0);

            // Calculate additional data
            if (!(sw[7] == 0 && sw[8] == 0 && sw[14] == 0)) {
                final T tloc = hl.multiply(HOUR_TO_RAD);
                final T tlx2 = tloc.add(tloc);
                final T tlx3 = tloc.add(tlx2);
                stloc  = tloc.sin();
                ctloc  = tloc.cos();
                s2tloc = tlx2.sin();
                c2tloc = tlx2.cos();
                s3tloc = tlx3.sin();
                c3tloc = tlx3.cos();
            } else {
                stloc  = zero;
                ctloc  = zero;
                s2tloc = zero;
                c2tloc = zero;
                s3tloc = zero;
                c3tloc = zero;
            }

        }

        /** Calculate temperatures and densities not including anomalous oxygen.
         *  <p>
         *  This method is the thermospheric portion of NRLMSISE-00 for alt > 72.5 km.
         *  </p>
         *  <p>NOTES ON INPUT VARIABLES:<br>
         *  Seconds, Local Time, and Longitude are used independently in the
         *  model and are not of equal importance for every situation.<br>
         *  For the most physically realistic calculation these three
         *  variables should be consistent (lst=sec/3600 + lon/15).<br>
         *  The Equation of Time departures from the above formula
         *  for apparent local time can be included if available but
         *  are of minor importance.<br><br>
         *
         *  f107 and f107A values used to generate the model correspond
         *  to the 10.7 cm radio flux at the actual distance of the Earth
         *  from the Sun rather than the radio flux at 1 AU. The following
         *  site provides both classes of values:<br>
         *  ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
         *
         *  f107, f107A, and ap effects are neither large nor well established below 80 km
         *  and these parameters should be set to 150., 150., and 4. respectively.
         *  </p>
         *  @param alt altitude (km)
         */
        void gts7(final T alt) {

            // Thermal diffusion coefficients for species
            final double[] alpha = {-0.38, 0.0, 0.0, 0.0, 0.17, 0.0, -0.38, 0.0, 0.0};
            // Altitude limits for net density computation for species
            final double[] altl  = {200.0, 300.0, 160.0, 250.0, 240.0, 450.0, 320.0, 450.0};
            // N2 mixed density
            final double xmm = PDM[2][4];

            /**** Exospheric temperature ****/
            T tinf = zero.add(PTM[0] * PT[0]);
            // Tinf variations not important below ZA or ZN[0]
            if (alt.getReal() > ZN1[0]) {
                tinf = tinf.multiply(globe7(PT).multiply(sw[16]).add(1));
            }
            setTemperature(EXOSPHERIC, tinf);

            // Gradient variations not important below ZN[4]
            T g0 = zero.add(PTM[3] * PS[0]);
            if (alt.getReal() > ZN1[4]) {
                g0 = g0.multiply(globe7(PS).multiply(sw[19]).add(1));
            }

            // Temperature at lower boundary
            T tlb = zero.add(PTM[1] * PD[3][0]);
            tlb = tlb.multiply(globe7(PD[3]).multiply(sw[17]).add(1));

            // Slope
            final T s = g0.divide(tinf.subtract(tlb));

            // Lower thermosphere temp variations not significant for density above 300 km
            meso_tn1[1]  = zero.add(PTM[6] * PTL[0][0]);
            meso_tn1[2]  = zero.add(PTM[2] * PTL[1][0]);
            meso_tn1[3]  = zero.add(PTM[7] * PTL[2][0]);
            meso_tn1[4]  = zero.add(PTM[4] * PTL[3][0]);
            meso_tgn1[1] = zero.add(PTM[8] * PMA[8][0]);
            if (alt.getReal() < 300.0) {
                final double r = PTM[4] * PTL[3][0];
                meso_tn1[1]  =  meso_tn1[1].divide(glob7s(PTL[0]).multiply(sw[18]         ).negate().add(1));
                meso_tn1[2]  =  meso_tn1[2].divide(glob7s(PTL[1]).multiply(sw[18]         ).negate().add(1));
                meso_tn1[3]  =  meso_tn1[3].divide(glob7s(PTL[2]).multiply(sw[18]         ).negate().add(1));
                meso_tn1[4]  =  meso_tn1[4].divide(glob7s(PTL[3]).multiply(sw[18] * sw[20]).negate().add(1));
                meso_tgn1[1] =  meso_tgn1[1].multiply(glob7s(PMA[8]).multiply(sw[18] * sw[20]).add(1));
                meso_tgn1[1] =  meso_tgn1[1].multiply(meso_tn1[4].multiply(meso_tn1[4]).divide(r * r));
            }

            /**** Temperature at altitude ****/
            setTemperature(ALTITUDE, densu(alt, zero.add(1.0), tinf, tlb, 0, 0, PTM[5], s));

            /**** N2 density ****/
            /*   Density variation factor at Zlb */
            final T g28 = globe7(PD[2]).multiply(sw[21]);
            /* Diffusive density at Zlb */
            final T db28 = g28.exp().multiply(PDM[2][0] * PD[2][0]);
            /* Diffusive density at Alt */
            T diffusiveDensity = densu(alt, db28, tinf, tlb, N2_MASS, alpha[2], PTM[5], s);
            setDensity(MOLECULAR_NITROGEN, diffusiveDensity);
            // Variation of turbopause height
            final T zhf = lat.multiply(DEG_TO_RAD).sin().
                            multiply(sw[5] * PDL[0][24] * FastMath.cos(DAY_TO_RAD * (doy - PT[13]))).
                            add(1).
                            multiply(PDL[1][24]);
            /* Turbopause */
            final T zh28  = zhf.multiply(PDM[2][2]);
            final double zhm28 = PDM[2][3] * PDL[1][5];
            /* Mixed density at Zlb */
            final T b28 = densu(zh28, db28, tinf, tlb, N2_MASS - xmm, alpha[2] - 1.0, PTM[5], s);
            if (sw[15] != 0 && alt.getReal() <= altl[2]) {
                /*  Mixed density at Alt */
                dm28 = densu(alt, b28, tinf, tlb, xmm, alpha[2], PTM[5], s);
                /*  Net density at Alt */
                setDensity(MOLECULAR_NITROGEN, dnet(diffusiveDensity, dm28, zhm28, xmm, N2_MASS));
            } else {
                dm28 = zero;
            }

            /**** He density ****/
            /*   Density variation factor at Zlb */
            final T g4 = globe7(PD[0]).multiply(sw[21]);
            /*  Diffusive density at Zlb */
            final T db04 = g4.exp().multiply(PDM[0][0] * PD[0][0]);
            /*  Diffusive density at Alt */
            diffusiveDensity = densu(alt, db04, tinf, tlb, HE_MASS, alpha[0], PTM[5], s);
            setDensity(HELIUM, diffusiveDensity);
            if (sw[15] != 0 && alt.getReal() < altl[0]) {
                /*  Turbopause */
                final double zh04 = PDM[0][2];
                /*  Mixed density at Zlb */
                final T b04 = densu(zero.add(zh04), db04, tinf, tlb, HE_MASS - xmm, alpha[0] - 1., PTM[5], s);
                /*  Mixed density at Alt */
                final T dm04 = densu(alt, b04, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm04 = zhm28;
                /*  Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm04, zhm04, xmm, HE_MASS);
                /*  Correction to specified mixing ratio at ground */
                final T rl = b28.multiply(PDM[0][1]).divide(b04).log();
                final double zc04 = PDM[0][4] * PDL[1][0];
                final double hc04 = PDM[0][5] * PDL[1][1];
                /*  Net density corrected at Alt */
                setDensity(HELIUM, diffusiveDensity.multiply(ccor(alt, rl, hc04, zc04)));
            }

            /**** O density ****/
            /* Density variation factor at Zlb */
            final T g16 = globe7(PD[1]).multiply(sw[21]);
            /* Diffusive density at Zlb */
            final T db16 = g16.exp().multiply(PDM[1][0] * PD[1][0]);
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db16, tinf, tlb, O_MASS, alpha[1], PTM[5], s);
            setDensity(ATOMIC_OXYGEN, diffusiveDensity);
            if (sw[15] != 0 && alt.getReal() < altl[1]) {
                /* Turbopause */
                final double zh16 = PDM[1][2];
                /* Mixed density at Zlb */
                final T b16 = densu(zero.add(zh16), db16, tinf, tlb, O_MASS - xmm, alpha[1] - 1.0, PTM[5], s);
                /* Mixed density at Alt */
                final T dm16 = densu(alt, b16, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm16 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm16, zhm16, xmm, O_MASS);
                final double rl = PDM[1][1] * PDL[1][16] * (1.0 + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
                final double hc16 = PDM[1][5] * PDL[1][3];
                final double zc16 = PDM[1][4] * PDL[1][2];
                final double hc216 = PDM[1][5] * PDL[1][4];
                diffusiveDensity = diffusiveDensity.multiply(ccor2(alt, rl, hc16, zc16, hc216));
                /* Chemistry correction */
                final double hcc16 = PDM[1][7] * PDL[1][13];
                final double zcc16 = PDM[1][6] * PDL[1][12];
                final double rc16  = PDM[1][3] * PDL[1][14];
                /* Net density corrected at Alt */
                setDensity(ATOMIC_OXYGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc16), hcc16, zcc16)));
            }

            /**** O2 density ****/
            /* Density variation factor at Zlb */
            final T g32 = globe7(PD[4]).multiply(sw[21]);
            /* Diffusive density at Zlb */
            final T db32 = g32.exp().multiply(PDM[3][0] * PD[4][0]);
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db32, tinf, tlb, O2_MASS, alpha[3], PTM[5], s);
            setDensity(MOLECULAR_OXYGEN, diffusiveDensity);
            if (sw[15] != 0) {
                if (alt.getReal() <= altl[3]) {
                    /* Turbopause */
                    final double zh32 = PDM[3][2];
                    /* Mixed density at Zlb */
                    final T b32 = densu(zero.add(zh32), db32, tinf, tlb, O2_MASS - xmm, alpha[3] - 1., PTM[5], s);
                    /* Mixed density at Alt */
                    final T dm32 = densu(alt, b32, tinf, tlb, xmm, 0., PTM[5], s);
                    final double zhm32 = zhm28;
                    /* Net density at Alt */
                    diffusiveDensity = dnet(diffusiveDensity, dm32, zhm32, xmm, O2_MASS);
                    /* Correction to specified mixing ratio at ground */
                    final T rl = b28.multiply(PDM[3][1]).divide(b32).log();
                    final double hc32 = PDM[3][5] * PDL[1][7];
                    final double zc32 = PDM[3][4] * PDL[1][6];
                    diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc32, zc32));
                }
                /* Correction for general departure from diffusive equilibrium above Zlb */
                final double hcc32  = PDM[3][7] * PDL[1][22];
                final double hcc232 = PDM[3][7] * PDL[0][22];
                final double zcc32  = PDM[3][6] * PDL[1][21];
                final double rc32   = PDM[3][3] * PDL[1][23] * (1. + sw[1] * PDL[0][23] * (f107a - FLUX_REF));
                /* Net density corrected at Alt */
                setDensity(MOLECULAR_OXYGEN, diffusiveDensity.multiply(ccor2(alt, rc32, hcc32, zcc32, hcc232)));
            }

            /**** Ar density ****/
            /* Density variation factor at Zlb */
            final T g40 = globe7(PD[5]).multiply(sw[21]);
            /* Diffusive density at Zlb */
            final T db40 = g40.exp().multiply(PDM[4][0] * PD[5][0]);
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db40, tinf, tlb, AR_MASS, alpha[4], PTM[5], s);
            setDensity(ARGON, diffusiveDensity);
            if (sw[15] != 0 && alt.getReal() <= altl[4]) {
                /* Turbopause */
                final double zh40 = PDM[4][2];
                /* Mixed density at Zlb */
                final T b40 = densu(zero.add(zh40), db40, tinf, tlb, AR_MASS - xmm, alpha[4] - 1., PTM[5], s);
                /* Mixed density at Alt */
                final T dm40 = densu(alt, b40, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm40 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm40, zhm40, xmm, AR_MASS);
                /* Correction to specified mixing ratio at ground */
                final T rl = b28.multiply(PDM[4][1]).divide(b40).log();
                final double hc40 = PDM[4][5] * PDL[1][9];
                final double zc40 = PDM[4][4] * PDL[1][8];
                /* Net density corrected at Alt */
                setDensity(ARGON, diffusiveDensity.multiply(ccor(alt, rl, hc40, zc40)));
            }

            /**** H density ****/
            /* Density variation factor at Zlb */
            final T g1 = globe7(PD[6]).multiply(sw[21]);
            /* Diffusive density at Zlb */
            final T db01 = g1.exp().multiply(PDM[5][0] * PD[6][0]);
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db01, tinf, tlb, H_MASS, alpha[6], PTM[5], s);
            setDensity(HYDROGEN, diffusiveDensity);
            if (sw[15] != 0 && alt.getReal() <= altl[6]) {
                /* Turbopause */
                final double zh01 = PDM[5][2];
                /* Mixed density at Zlb */
                final T b01 = densu(zero.add(zh01), db01, tinf, tlb, H_MASS - xmm, alpha[6] - 1., PTM[5], s);
                /* Mixed density at Alt */
                final T dm01 = densu(alt, b01, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm01 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm01, zhm01, xmm, H_MASS);
                /* Correction to specified mixing ratio at ground */
                final T rl = b28.multiply(PDM[5][1] * FastMath.sqrt(PDL[1][17] * PDL[1][17])).divide(b01).log();
                final double hc01 = PDM[5][5] * PDL[1][11];
                final double zc01 = PDM[5][4] * PDL[1][10];
                diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc01, zc01));
                /* Chemistry correction */
                final double hcc01 = PDM[5][7] * PDL[1][19];
                final double zcc01 = PDM[5][6] * PDL[1][18];
                final double rc01 = PDM[5][3] * PDL[1][20];
                /* Net density corrected at Alt */
                setDensity(HYDROGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc01), hcc01, zcc01)));
            }

            /**** N density ****/
            /* Density variation factor at Zlb */
            final T g14 = globe7(PD[7]).multiply(sw[21]);
            /* Diffusive density at Zlb */
            final T db14 = g14.exp().multiply(PDM[6][0] * PD[7][0]);
            /* Diffusive density at Alt */
            diffusiveDensity = densu(alt, db14, tinf, tlb, N_MASS, alpha[7], PTM[5], s);
            setDensity(ATOMIC_NITROGEN, diffusiveDensity);
            if (sw[15] != 0 && alt.getReal() <= altl[7]) {
                /* Turbopause */
                final double zh14 = PDM[6][2];
                /* Mixed density at Zlb */
                final T b14 = densu(zero.add(zh14), db14, tinf, tlb, N_MASS - xmm, alpha[7] - 1., PTM[5], s);
                /* Mixed density at Alt */
                final T dm14 = densu(alt, b14, tinf, tlb, xmm, 0., PTM[5], s);
                final double zhm14 = zhm28;
                /* Net density at Alt */
                diffusiveDensity = dnet(diffusiveDensity, dm14, zhm14, xmm, N_MASS);
                /* Correction to specified mixing ratio at ground */
                final T rl = b28.multiply(PDM[6][1] * PDL[0][2]).divide(b14).log();
                final double hc14 = PDM[6][5] * PDL[0][1];
                final double zc14 = PDM[6][4] * PDL[0][0];
                diffusiveDensity = diffusiveDensity.multiply(ccor(alt, rl, hc14, zc14));
                /* Chemistry correction */
                final double hcc14 = PDM[6][7] * PDL[0][4];
                final double zcc14 = PDM[6][6] * PDL[0][3];
                final double rc14 = PDM[6][3] * PDL[0][5];
                /* Net density corrected at Alt */
                setDensity(ATOMIC_NITROGEN, diffusiveDensity.multiply(ccor(alt, zero.add(rc14), hcc14, zcc14)));
            }

            /**** Anomalous O density ****/
            final T g16h = globe7(PD[8]).multiply(sw[21]);
            final T db16h = g16h.exp().multiply(PDM[7][0] * PD[8][0]);
            final double tho   = PDM[7][9] * PDL[0][6];
            diffusiveDensity = densu(alt, db16h, zero.add(tho), zero.add(tho), O_MASS, alpha[8], PTM[5], s);
            final double zsht = PDM[7][5];
            final double zmho = PDM[7][4];
            final T zsho = scalh(zmho, O_MASS, tho);
            diffusiveDensity = diffusiveDensity.multiply(alt.negate().add(zmho).divide(zsht).exp().subtract(1).multiply(-zsht).divide(zsho).exp());
            setDensity(ANOMALOUS_OXYGEN, diffusiveDensity);

            // Convert densities from cm-3 to m-3
            for (int i = 0; i < 9; i++) {
                setDensity(i, getDensity(i).multiply(1.0e+06));
            }

            /**** Total mass density ****/
            final T tmd =     getDensity(HELIUM)            .multiply(HE_MASS).
                          add(getDensity(ATOMIC_OXYGEN)     .multiply( O_MASS)).
                          add(getDensity(MOLECULAR_NITROGEN).multiply(N2_MASS)).
                          add(getDensity(MOLECULAR_OXYGEN)  .multiply(O2_MASS)).
                          add(getDensity(ARGON)             .multiply(AR_MASS)).
                          add(getDensity(HYDROGEN)          .multiply( H_MASS)).
                          add(getDensity(ATOMIC_NITROGEN)   .multiply( N_MASS)).
                          multiply(AMU);
            setDensity(TOTAL_MASS, tmd);

        }

        /** Calculate temperatures and densities not including anomalous oxygen.
         *  <p>NOTES ON INPUT VARIABLES:<br>
         *  Seconds, Local Time, and Longitude are used independently in the
         *  model and are not of equal importance for every situation.<br>
         *  For the most physically realistic calculation these three
         *  variables should be consistent (lst=sec/3600 + lon/15).<br>
         *  The Equation of Time departures from the above formula
         *  for apparent local time can be included if available but
         *  are of minor importance.<br><br>
         *
         *  f107 and f107A values used to generate the model correspond
         *  to the 10.7 cm radio flux at the actual distance of the Earth
         *  from the Sun rather than the radio flux at 1 AU. The following
         *  site provides both classes of values:<br>
         *  ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br><br>
         *
         *  f107, f107A, and ap effects are neither large nor well established below 80 km
         *  and these parameters should be set to 150., 150., and 4. respectively.
         *  </p>
         *  @param alt altitude (km)
         */
        void gtd7(final T alt) {

            // Calculates for thermosphere/mesosphere (above ZN2[0])
            final T altt = (alt.getReal() > ZN2[0]) ? alt : zero.add(ZN2[0]);
            gts7(altt);
            if (alt.getReal() >= ZN2[0]) {
                return;
            }

            // Calculates for lower mesosphere/upper stratosphere (between ZN2[0] and ZN3[0]):
            // Temperature at nodes and gradients at end nodes
            // Inverse temperature a linear function of spherical harmonics
            final double r = PMA[2][0] * PAVGM[2];
            meso_tgn2[0] = meso_tgn1[1];
            meso_tn2[0]  = meso_tn1[4];
            meso_tn2[1]  = glob7s(PMA[0]).multiply(sw[20]         ).negate().add(1).reciprocal().multiply(PMA[0][0] * PAVGM[0]);
            meso_tn2[2]  = glob7s(PMA[1]).multiply(sw[20]         ).negate().add(1).reciprocal().multiply(PMA[1][0] * PAVGM[1]);
            meso_tn2[3]  = glob7s(PMA[2]).multiply(sw[20] * sw[22]).negate().add(1).reciprocal().multiply(PMA[2][0] * PAVGM[2]);
            meso_tgn2[1] = glob7s(PMA[9]).multiply(sw[20] * sw[22]).add(1).multiply(PMA[9][0] * PAVGM[8]).
                           multiply(meso_tn2[3]).multiply(meso_tn2[3]).divide(r * r);
            meso_tn3[0]  = meso_tn2[3];

            // Calculates for lower stratosphere and troposphere (below ZN3[0])
            // Temperature at nodes and gradients at end nodes
            // Inverse temperature a linear function of spherical harmonics
            if (alt.getReal() < ZN3[0]) {
                final double q = PMA[6][0] * PAVGM[6];
                meso_tgn3[0] = meso_tgn2[1];
                meso_tn3[1]  = glob7s(PMA[3]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[3][0] * PAVGM[3]);
                meso_tn3[2]  = glob7s(PMA[4]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[4][0] * PAVGM[4]);
                meso_tn3[3]  = glob7s(PMA[5]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[5][0] * PAVGM[5]);
                meso_tn3[4]  = glob7s(PMA[6]).multiply(sw[22]).negate().add(1).reciprocal().multiply(PMA[6][0] * PAVGM[6]);
                meso_tgn3[1] = glob7s(PMA[7]).multiply(sw[22])         .add(1).multiply(PMA[7][0] * PAVGM[7]).
                               multiply(meso_tn3[4]).multiply(meso_tn3[4]).divide(q * q);

            }

            // Linear transition to full mixing below ZN2[0]
            final T dmc = (alt.getReal() > ZMIX) ?
                           alt.subtract(ZN2[0]).divide(ZN2[0] - ZMIX).add(1) :
                           zero;
            final T dz28 = getDensity(MOLECULAR_NITROGEN);

            // N2 density
            final T dm28m = dm28.multiply(1.0e+06);
            T dmr = dz28.divide(dm28m).subtract(1);
            T dst = densm(alt, dm28m, PDM[2][4]).multiply(dmr.multiply(dmc).add(1));
            setDensity(MOLECULAR_NITROGEN, dst);

            // HE density
            dmr = getDensity(HELIUM).divide(dz28.multiply(PDM[0][1])).subtract(1);
            dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[0][1]).multiply(dmr.multiply(dmc).add(1));
            setDensity(HELIUM, dst);

            // O density
            setDensity(ATOMIC_OXYGEN, zero);
            setDensity(ANOMALOUS_OXYGEN, zero);

            // O2 density
            dmr = getDensity(MOLECULAR_OXYGEN).divide(dz28.multiply(PDM[3][1])).subtract(1);
            dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[3][1]).multiply(dmr.multiply(dmc).add(1));
            setDensity(MOLECULAR_OXYGEN, dst);

            // AR density
            dmr = getDensity(ARGON).divide(dz28.multiply(PDM[4][1])).subtract(1);
            dst = getDensity(MOLECULAR_NITROGEN).multiply(PDM[4][1]).multiply(dmr.multiply(dmc).add(1));
            setDensity(ARGON, dst);

            // H density
            setDensity(HYDROGEN, zero);

            // N density
            setDensity(ATOMIC_NITROGEN, zero);

            // Total mass density
            final T tmd =       getDensity(HELIUM)            .multiply(HE_MASS).
                            add(getDensity(ATOMIC_OXYGEN)     .multiply( O_MASS)).
                            add(getDensity(MOLECULAR_NITROGEN).multiply(N2_MASS)).
                            add(getDensity(MOLECULAR_OXYGEN)  .multiply(O2_MASS)).
                            add(getDensity(ARGON)             .multiply(AR_MASS)).
                            add(getDensity(HYDROGEN)          .multiply( H_MASS)).
                            add(getDensity(ATOMIC_NITROGEN)   .multiply( N_MASS)).
                            multiply(AMU);
            setDensity(TOTAL_MASS, tmd);

            // Temperature at altitude
            setTemperature(ALTITUDE, densm(alt, field.getOne(), 0));

        }

        /** Calculate temperatures and densities including anomalous oxygen.
         *  <p></p>
         *  <p>NOTES ON INPUT VARIABLES:<br>
         *  Seconds, Local Time, and Longitude are used independently in the
         *  model and are not of equal importance for every situation.<br>
         *  For the most physically realistic calculation these three
         *  variables should be consistent (lst=sec/3600 + lon/15).<br>
         *  The Equation of Time departures from the above formula
         *  for apparent local time can be included if available but
         *  are of minor importance.<br>
         *  <br>
         *  f107 and f107A values used to generate the model correspond
         *  to the 10.7 cm radio flux at the actual distance of the Earth
         *  from the Sun rather than the radio flux at 1 AU. The following
         *  site provides both classes of values:<br>
         *  ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/<br>
         *  <br>
         *  f107, f107A, and ap effects are neither large nor well established below 80 km
         *  and these parameters should be set to 150., 150., and 4. respectively.
         *  </p>
         *  @param alt altitude (km)
         */
        void gtd7d(final T alt) {

            // Compute densities and temperatures
            gtd7(alt);

            // Update the total mass density with anomalous oxygen contribution
            final T dTot = getDensity(TOTAL_MASS).add(getDensity(ANOMALOUS_OXYGEN).multiply( AMU * O_MASS));
            setDensity(TOTAL_MASS, dTot);

        }

        /** Set one density.
         * @param index one of the nine elements :
         * <ul>
         * <li>{@link #HELIUM}</li>
         * <li>{@link #ATOMIC_OXYGEN}</li>
         * <li>{@link #MOLECULAR_NITROGEN}</li>
         * <li>{@link #MOLECULAR_OXYGEN}</li>
         * <li>{@link #ARGON}</li>
         * <li>{@link #TOTAL_MASS}</li>
         * <li>{@link #HYDROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * </ul>
         * @param d the value of density to set
         */
        void setDensity(final int index, final T d) {
            densities[index] = d;
        }

        /** Set one temperature.
         * @param index one of the two elements :
         * <ul>
         * <li>{@link #EXOSPHERIC}</li>
         * <li>{@link #ALTITUDE}</li>
         * </ul>
         * @param t the value of temperature to set
         */
        void setTemperature(final int index, final T t) {
            temperatures[index] = t;
        }

        /** Get one of the stored densities.
         * @param index one of the nine elements :
         * <ul>
         * <li>{@link #HELIUM}</li>
         * <li>{@link #ATOMIC_OXYGEN}</li>
         * <li>{@link #MOLECULAR_NITROGEN}</li>
         * <li>{@link #MOLECULAR_OXYGEN}</li>
         * <li>{@link #ARGON}</li>
         * <li>{@link #TOTAL_MASS}</li>
         * <li>{@link #HYDROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * <li>{@link #ATOMIC_NITROGEN}</li>
         * </ul>
         * @return the requested density
         */
        public T getDensity(final int index) {
            return densities[index];
        }

        /** Calculate G(L) function with upper thermosphere parameters.
         *  @param p array of parameters
         *  @return G(L) value
         */
        private T globe7(final double[] p) {

            final T[] t = MathArrays.buildArray(field, 14);
            final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
            final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
            final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
            final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));

            // F10.7 effect
            final double df  = f107  - f107a;
            final double dfa = f107a - FLUX_REF;
            t[0] = zero.add(p[19] * df * (1.0 + p[59] * dfa) +
                            p[20] * df * df +
                            p[21] * dfa +
                            p[29] * dfa * dfa);

            final double f1 = 1.0 + (p[47] * dfa + p[19] * df + p[20] * df * df) * swc[1];
            final double f2 = 1.0 + (p[49] * dfa + p[19] * df + p[20] * df * df) * swc[1];

            // Time independent
            t[1] =     plg[0][2].multiply(p[ 1]).
                   add(plg[0][4].multiply(p[ 2])).
                   add(plg[0][6].multiply(p[22])).
                   add(plg[0][2].multiply(p[14] * dfa * swc[1])).
                   add(plg[0][1].multiply(p[26]));

            // Symmetrical annual
            t[2] = zero.add(p[18] * cd32);

            // Symmetrical semiannual
            t[3] = plg[0][2].multiply(p[16]).add(p[15]).multiply(cd18);

            // Asymmetrical annual
            t[4] = plg[0][1].multiply(p[9]).add(plg[0][3].multiply(p[10])).multiply(f1 * cd14);

            // Asymmetrical semiannual
            t[5] = plg[0][1].multiply(p[37] * cd39);

            // Diurnal
            if (sw[7] != 0) {
                final T t71 = plg[1][2].multiply(p[11] * cd14 * swc[5]);
                final T t72 = plg[1][2].multiply(p[12] * cd14 * swc[5]);
                t[6] =      plg[1][1].multiply(p[3]).add(plg[1][3].multiply(p[4])).add(plg[1][5].multiply(p[27])).add(t71).multiply(ctloc).
                        add(plg[1][1].multiply(p[6]).add(plg[1][3].multiply(p[7])).add(plg[1][5].multiply(p[28])).add(t72).multiply(stloc)).
                        multiply(f2);
            }

            // Semidiurnal
            if (sw[8] != 0) {
                final T t81 = plg[2][3].multiply(p[23]).add(plg[2][5].multiply(p[35])).multiply(cd14 * swc[5]);
                final T t82 = plg[2][3].multiply(p[33]).add(plg[2][5].multiply(p[36])).multiply(cd14 * swc[5]);
                t[7] =     plg[2][2].multiply(p[5]).add(plg[2][4].multiply(p[41])).add(t81).multiply(c2tloc).
                       add(plg[2][2].multiply(p[8]).add(plg[2][4].multiply(p[42])).add(t82).multiply(s2tloc)).
                       multiply(f2);
            }

            // Terdiurnal
            if (sw[14] != 0) {
                t[13] =     plg[3][3].multiply(p[39]).add(plg[3][4].multiply(p[93]).add(plg[3][6].multiply(p[46])).multiply(cd14 * swc[5])).multiply(s3tloc).
                        add(plg[3][3].multiply(p[40]).add(plg[3][4].multiply(p[94]).add(plg[3][6].multiply(p[48])).multiply(cd14 * swc[5])).multiply(c3tloc)).
                        multiply(f2);
            }

            // magnetic activity based on daily ap
            if (sw[9] == -1) {
                if (p[51] != 0) {
                    final T exp1 = lat.abs().negate().add(LAT_REF).multiply(p[138]).add(1).
                                    reciprocal().multiply(-10800.0 * FastMath.abs(p[51])).
                                    exp();
                    final double p24 = FastMath.max(p[24], 1.0e-4);
                    apt = sg0(min(0.99999, exp1), p24, p[25]);
                    t[8] =      plg[0][2].multiply(p[96]).add(plg[0][4].multiply(p[54])).add(p[50]).
                           add((plg[0][1].multiply(p[125]).add(plg[0][3].multiply(p[126])).add(plg[0][5].multiply(p[127]))).multiply(cd14 * swc[5])).
                           add((plg[1][1].multiply(p[128]).add(plg[1][3].multiply(p[129])).add(plg[1][5].multiply(p[130]))).multiply(swc[7]).multiply(hl.subtract(p[131]).multiply(HOUR_TO_RAD).cos())).
                           multiply(apt);
                }
            } else {
                final double apd = ap[0] - 4.0;
                final double p44 = (p[43] < 0.) ? 1.0E-5 : p[43];
                final double p45 = p[44];
                apdf = apd + (p45 - 1.0) * (apd + (FastMath.exp(-p44 * apd) - 1.0) / p44);
                if (sw[9] != 0) {
                    t[8] =      plg[0][2].multiply(p[45]).add(plg[0][4].multiply(p[34])).add(p[32]).
                           add((plg[0][1].multiply(p[100]).add(plg[0][3].multiply(p[101])).add(plg[0][5].multiply(p[102]))).multiply(cd14 * swc[5])).
                           add((plg[1][1].multiply(p[121]).add(plg[1][3].multiply(p[122])).add(plg[1][5].multiply(p[123]))).multiply(swc[7]).multiply(hl.subtract(p[124]).multiply(HOUR_TO_RAD).cos())).
                           multiply(apdf);
                }
            }

            if (sw[10] != 0) {
                final T lonr = lon.multiply(DEG_TO_RAD);
                // Longitudinal
                if (sw[11] != 0) {
                    t[10] =         plg[1][2].multiply(p[ 64]) .add(plg[1][4].multiply(p[ 65])).add(plg[1][6].multiply(p[ 66])).
                                add(plg[1][1].multiply(p[103])).add(plg[1][3].multiply(p[104])).add(plg[1][5].multiply(p[105])).
                                add((plg[1][1].multiply(p[109])).add(plg[1][3].multiply(p[110])).add(plg[1][5].multiply(p[111])).multiply(swc[5] * cd14)).
                                multiply(lonr.cos()).
                            add(    plg[1][2].multiply(p[ 90]) .add(plg[1][4].multiply(p[ 91])).add(plg[1][6].multiply(p[ 92])).
                                add(plg[1][1].multiply(p[106])).add(plg[1][3].multiply(p[107])).add(plg[1][5].multiply(p[108])).
                                add((plg[1][1].multiply(p[112])).add(plg[1][3].multiply(p[113])).add(plg[1][5].multiply(p[114])).multiply(swc[5] * cd14)).
                                multiply(lonr.sin())).
                            multiply(1.0 + p[80] * dfa * swc[1]);
                }

                // ut and mixed ut, longitude
                if (sw[12] != 0) {
                    t[11] =          plg[0][1].multiply(p[95]).add(1).multiply(1.0 + p[81] * dfa * swc[1]).
                            multiply(plg[0][1].multiply(p[119] * swc[5] * cd14).add(1)).
                            multiply(plg[0][1].multiply(p[68]).add(plg[0][3].multiply(p[69])).add(plg[0][5].multiply(p[70]))).
                            multiply(sec.subtract(p[71]).multiply(SEC_TO_RAD).cos());
                    t[11] = t[11].
                            add(plg[2][3].multiply(p[76]).add(plg[2][5].multiply(p[77])).add(plg[2][7].multiply(p[78])).
                                multiply(swc[11] * (1.0 + p[137] * dfa * swc[1])).
                                multiply(sec.subtract(p[79]).multiply(SEC_TO_RAD).add(lonr.multiply(2)).cos()));
                }

                /* ut, longitude magnetic activity */
                if (sw[13] != 0) {
                    if (sw[9] == -1) {
                        if (p[51] != 0.) {
                            t[12] = apt.multiply(swc[11]).multiply(plg[0][1].multiply(p[132]).add(1)).
                                    multiply(plg[1][2].multiply(p[52]).add(plg[1][4].multiply(p[98])).add(plg[1][6].multiply(p[67]))).
                                    multiply(lon.subtract(p[97]).multiply(DEG_TO_RAD).cos()).
                                    add(apt.multiply(swc[11] * swc[5] * cd14).
                                        multiply(plg[1][1].multiply(p[133]).add(plg[1][3].multiply(p[134])).add(plg[1][5].multiply(p[135]))).
                                        multiply(lon.subtract(p[136]).multiply(DEG_TO_RAD).cos())).
                                    add(apt.multiply(swc[12]).
                                        multiply(plg[0][1].multiply(p[55]).add(plg[0][3].multiply(p[56])).add(plg[0][5].multiply(p[57]))).
                                        multiply(sec.subtract(p[58]).multiply(SEC_TO_RAD).cos()));
                        }
                    } else {
                        t[12] = plg[0][1].multiply(p[120]).add(1).multiply(apdf * swc[11]).
                                multiply(plg[1][2].multiply(p[60]).add(plg[1][4].multiply(p[61])).add(plg[1][6].multiply(p[62]))).
                                multiply(lon.subtract(p[63]).multiply(DEG_TO_RAD).cos()).
                                add(plg[1][1].multiply(p[115]).add(plg[1][3].multiply(p[116])).add(plg[1][5].multiply(p[117])).
                                    multiply(apdf * swc[11] * swc[5] * cd14).
                                    multiply(lon.subtract(p[118]).multiply(DEG_TO_RAD).cos())).
                                add(plg[0][1].multiply(p[83]).add(plg[0][3].multiply(p[84])).add(plg[0][5].multiply(p[85])).
                                    multiply(apdf * swc[12]).
                                    multiply(sec.subtract(p[75]).multiply(SEC_TO_RAD).cos()));
                    }
                }
            }

            // Sum all effects (params not used: 82, 89, 99, 139-149)
            T tinf = zero.add(p[30]);
            for (int i = 0; i < 14; i++) {
                tinf = tinf.add(t[i].multiply(FastMath.abs(sw[i + 1])));
            }

            // Return G(L)
            return tinf;

        }

        /** Calculate G(L) function with lower atmosphere parameters.
         *  @param p array of parameters
         *  @return G(L) value
         */
        private T glob7s(final double[] p) {

            final T[] t = MathArrays.buildArray(field, 14);
            final double cd32 = FastMath.cos(DAY_TO_RAD * (doy - p[31]));
            final double cd18 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[17]));
            final double cd14 = FastMath.cos(DAY_TO_RAD * (doy - p[13]));
            final double cd39 = FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[38]));

            // F10.7 effect
            t[0] = zero.add(p[21] * (f107a - FLUX_REF));

            // Time independent
            t[1] =     plg[0][2].multiply(p[1]).
                   add(plg[0][4].multiply(p[2])).
                   add(plg[0][6].multiply(p[22])).
                   add(plg[0][1].multiply(p[26])).
                   add(plg[0][3].multiply(p[14])).
                   add(plg[0][5].multiply(p[59]));

            // Symmetrical annual
            t[2] = plg[0][2].multiply(p[47]).add(plg[0][4].multiply(p[29])).add(p[18]).multiply(cd32);

            // Symmetrical semiannual
            t[3] = plg[0][2].multiply(p[16]).add(plg[0][4].multiply(p[30])).add(p[15]).multiply(cd18);

            // Asymmetrical annual
            t[4] = plg[0][1].multiply(p[9]).add(plg[0][3].multiply(p[10])).add(plg[0][5].multiply(p[20])).multiply(cd14);

            // Asymmetrical semiannual
            t[5] = plg[0][1].multiply(p[37]).multiply(cd39);

            // Diurnal
            if (sw[7] != 0) {
                final T t71 = plg[1][2].multiply(p[11]).multiply(cd14 * swc[5]);
                final T t72 = plg[1][2].multiply(p[12]).multiply(cd14 * swc[5]);
                t[6] =     plg[1][1].multiply(p[3]).add(plg[1][3].multiply(p[4])).add(t71).multiply(ctloc).
                       add(plg[1][1].multiply(p[6]).add(plg[1][3].multiply(p[7])).add(t72).multiply(stloc));
            }

            // Semidiurnal
            if (sw[8] != 0) {
                final T t81 = plg[2][3].multiply(p[23]).add(plg[2][5].multiply(p[35])).multiply(cd14 * swc[5]);
                final T t82 = plg[2][3].multiply(p[33]).add(plg[2][5].multiply(p[36])).multiply(cd14 * swc[5]);
                t[7] =     plg[2][2].multiply(p[5]).add(plg[2][4].multiply(p[41])).add(t81).multiply(c2tloc).
                       add(plg[2][2].multiply(p[8]).add(plg[2][4].multiply(p[42])).add(t82).multiply(s2tloc));
            }

            // Terdiurnal
            if (sw[14] != 0) {
                t[13] = plg[3][3].multiply(p[39]).multiply(s3tloc).add(plg[3][3].multiply(p[40]).multiply(c3tloc));
            }

            // Magnetic activity
            if (sw[9] == 1) {
                t[8] = plg[0][2].multiply(p[45] * swc[2]).add(p[32]).multiply(apdf);
            } else if (sw[9] == -1) {
                t[8] = plg[0][2].multiply(p[96] * swc[2]).add(p[50]).multiply(apt);
            }

            // Longitudinal
            if (!(sw[10] == 0 || sw[11] == 0)) {
                final T lonr = lon.multiply(DEG_TO_RAD);
                t[10] = plg[0][1].multiply(p[80] * swc[5] * FastMath.cos(DAY_TO_RAD * (doy - p[81])) +
                                           p[85] * swc[6] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[86]))).
                       add(1.0 +
                           p[83] * swc[3] * FastMath.cos(DAY_TO_RAD * (doy - p[84])) +
                           p[87] * swc[4] * FastMath.cos(2.0 * DAY_TO_RAD * (doy - p[88]))).
                       multiply(    plg[1][2].multiply(p[64]).
                                add(plg[1][4].multiply(p[65])).
                                add(plg[1][6].multiply(p[66])).
                                add(plg[1][1].multiply(p[74])).
                                add(plg[1][3].multiply(p[75])).
                                add(plg[1][5].multiply(p[76])).multiply(lonr.cos()).
                          add(      plg[1][2].multiply(p[90]).
                                add(plg[1][4].multiply(p[91])).
                                add(plg[1][6].multiply(p[92])).
                                add(plg[1][1].multiply(p[77])).
                                add(plg[1][3].multiply(p[78])).
                                add(plg[1][5].multiply(p[79])).multiply(lonr.sin())));
            }

            // Sum all effects
            T gl = zero;
            for (int i = 0; i < 14; i++) {
                gl = gl.add(t[i].multiply(FastMath.abs(sw[i + 1])));
            }

            // Return G(L)
            return gl;
        }

        /** Implements sg0 function (Eq. A24a).
         * @param ex ex
         * @param p24 abs(p[24])
         * @param p25 p[25]
         * @return sg0
         */
        private T sg0(final T ex, final double p24, final double p25) {
            final double g01 = g0(ap[1], p24, p25);
            final double g02 = g0(ap[2], p24, p25);
            final double g03 = g0(ap[3], p24, p25);
            final double g04 = g0(ap[4], p24, p25);
            final double g05 = g0(ap[5], p24, p25);
            final double g06 = g0(ap[6], p24, p25);
            final T ex2      = ex.multiply(ex);
            final T ex3      = ex.multiply(ex2);
            final T ex4      = ex2.multiply(ex2);
            final T ex8      = ex4.multiply(ex4);
            final T ex12     = ex4.multiply(ex8);
            final T g234     = ex.multiply(g02).add(ex2.multiply(g03)).add(ex3.multiply(g04));
            final T g56      = ex4.multiply(g05).add(ex12.multiply(g06));
            final T ex19     = ex3.multiply(ex4).multiply(ex12);
            final T omex     = ex.negate().add(1);
            final T sumex    = ex19.negate().add(1).divide(omex).multiply(ex.sqrt()).add(1);
            return ex8.negate().add(1).multiply(g56).divide(omex).add(g234).add(g01).divide(sumex);
        }

        /** Implements go function (Eq. A24d).
         * @param apI 3 hrs ap
         * @param p24 abs(p[24])
         * @param p25 p[25]
         * @return go
         */
        private double g0(final double apI, final double p24, final double p25) {
            final double am4 = apI - 4.0;
            return am4 + (p25 - 1.0) * (am4 + (FastMath.exp(-p24 * am4) - 1.0) / p24);
        }

        /** Calculates chemistry/dissociation correction for MSIS models.
         * @param alt altitude
         * @param r target ratio
         * @param h1 transition scale length
         * @param zh altitude of 1/2 R
         * @return correction
         */
        private T ccor(final T alt, final T r, final double h1, final double zh) {
            final T e = alt.subtract(zh).divide(h1);
            if (e.getReal() > 70.) {
                return field.getOne();
            } else if (e.getReal() < -70.) {
                return r.exp();
            } else {
                return r.divide(e.exp().add(1)).exp();
            }
        }


        /** Calculates O & O2 chemistry/dissociation correction for MSIS models.
         * @param alt altitude
         * @param r target ratio
         * @param h1 transition scale length
         * @param zh altitude of 1/2 R
         * @param h2 transition scale length
         * @return correction
         */
        private T ccor2(final T alt, final double r, final double h1, final double zh, final double h2) {
            final T e1 = alt.subtract(zh).divide(h1);
            final T e2 = alt.subtract(zh).divide(h2);
            if ((e1.getReal() > 70.) || (e2.getReal() > 70.)) {
                return field.getOne();
            } else if ((e1.getReal() < -70.) && (e2.getReal() < -70.)) {
                return zero.add(FastMath.exp(r));
            } else {
                final T ex1 = e1.exp();
                final T ex2 = e2.exp();
                return ex1.add(ex2).multiply(0.5).add(1).reciprocal().multiply(r).exp();
            }
        }

        /** Calculates scale height.
         * @param alt altitude
         * @param xm species molecular weight
         * @param temp temperature
         * @return scale height (km)
         */
        private T scalh(final double alt, final double xm, final double temp) {
            // Gravity at altitude
            final T denom = rlat.reciprocal().multiply(alt).add(1);
            final T galt = glat.divide(denom.multiply(denom));
            return galt.reciprocal().multiply(R_GAS * temp / xm);
        }

        /** Calculates turbopause correction for MSIS models.
         * @param dd diffusive density
         * @param dm full mixed density
         * @param zhm transition scale length
         * @param xmm full mixed molecular weight
         * @param xm species molecular weight
         * @return combined density
         */
        private T dnet(final T dd, final T dm, final double zhm, final double xmm, final double xm) {
            if (!(dm.getReal() > 0 && dd.getReal() > 0)) {
                T ddd = dd;
                if (dd.getReal() == 0 && dm.getReal() == 0) {
                    ddd = field.getOne();
                }
                if (dm.getReal() == 0) {
                    return ddd;
                }
                if (dd.getReal() == 0) {
                    return dm;
                }
            }

            final double a  = zhm / (xmm - xm);
            final T ylog = dm.divide(dd).log().multiply(a);
            if (ylog.getReal() < -10.) {
                return dd;
            } else if (ylog.getReal() > 10.) {
                return dm;
            } else {
                return ylog.exp().add(1).pow(1.0 / a).multiply(dd);
            }
        }

        /** Integrate cubic spline function from xa[0] to x.
         * <p>ADAPTED FROM NUMERICAL RECIPES</p>
         * @param xa array of abscissas in ascending order
         * @param ya array of ordinates in ascending order by xa
         * @param y2a array of second derivatives in ascending order by xa
         * @param x abscissa end point
         * @return integral value
         */
        private T splini(final T[] xa, final T[] ya, final T[] y2a, final T x) {
            final int n = xa.length;
            T yi = zero;
            int klo = 0;
            int khi = 1;
            while (x.getReal() > xa[klo].getReal() && khi < n) {
                T xx = x;
                if (khi < n - 1) {
                    xx = (x.getReal() < xa[khi].getReal()) ? x : xa[khi];
                }
                final T h = xa[khi].subtract(xa[klo]);
                final T a = xa[khi].subtract(xx).divide(h);
                final T b = xx.subtract(xa[klo]).divide(h);
                final T a2 = a.multiply(a);
                final T b2 = b.multiply(b);

                final T z =
                           a2.divide(2).subtract(a2.multiply(a2).add(1).divide(4)).multiply(y2a[klo]).
                           add(b2.multiply(b2).divide(4).subtract(b2.divide(2)).multiply(y2a[khi]));
                yi = yi.add(    a2.negate().add(1).multiply(ya[klo]).divide(2).
                            add(b2.multiply(ya[khi]).divide(2)).
                            add(z.multiply(h).multiply(h).divide(6)).
                            multiply(h));
                klo++;
                khi++;
            }
            return yi;
        }

        /** Calculate cubic spline interpolated value.
         * <p>ADAPTED FROM NUMERICAL RECIPES</p>
         * @param xa array of abscissas in ascending order
         * @param ya array of ordinates in ascending order by xa
         * @param y2a array of second derivatives in ascending order by xa
         * @param x abscissa for interpolation
         * @return interpolated value
         */
        private T splint(final T[] xa, final T[] ya, final T[] y2a, final T x) {
            final int n = xa.length;
            int klo = 0;
            int khi = n - 1;
            while (khi - klo > 1) {
                final int k = (khi + klo) >>> 1;
                if (xa[k].getReal() > x.getReal()) {
                    khi = k;
                } else {
                    klo = k;
                }
            }
            final T h = xa[khi].subtract(xa[klo]);
            final T a = xa[khi].subtract(x).divide(h);
            final T b = x.subtract(xa[klo]).divide(h);
            return a.multiply(ya[klo]).add(b.multiply(ya[khi])).
                   add((    a.multiply(a).multiply(a).subtract(a).multiply(y2a[klo]).
                        add(b.multiply(b).multiply(b).subtract(b).multiply(y2a[khi]))
                       ).multiply(h).multiply(h).divide(6));
        }

        /** Calculate 2nd derivatives of cubic spline interpolation function.
         * <p>ADAPTED FROM NUMERICAL RECIPES</p>
         * @param x array of abscissas in ascending order
         * @param y array of ordinates in ascending order by x
         * @param yp1 derivative at x[0] (2nd derivatives null if > 1E30)
         * @param ypn derivative at x[n-1] (2nd derivatives null if > 1E30)
         * @return array of second derivatives
         */
        private T[] spline(final T[] x, final T[] y, final T yp1, final T ypn) {
            final int n = x.length;
            final T[] y2 = MathArrays.buildArray(field, n);
            final T[] u  = MathArrays.buildArray(field, n);

            if (yp1.getReal() < 1e+30) {
                y2[0] = zero.add(-0.5);
                final T dx = x[1].subtract(x[0]);
                final T dy = y[1].subtract(y[0]);
                u[0]  = dx.reciprocal().multiply(3.0).multiply(dy.divide(dx).subtract(yp1));
            }
            for (int i = 1; i < n - 1; i++) {
                final T dx0m = x[i].subtract(x[i - 1]);
                final T dy0m = y[i].subtract(y[i - 1]);
                final T dxpm = x[i + 1].subtract(x[i - 1]);
                final T dxp0 = x[i + 1].subtract(x[i]);
                final T dyp0 = y[i + 1].subtract(y[i]);
                final T sig = dx0m.divide(dxpm);
                final T p = sig.multiply(y2[i - 1]).add(2.0);
                y2[i] = sig.subtract(1.0).divide(p);
                u[i] = dyp0.divide(dxp0).subtract(dy0m.divide(dx0m)).multiply(6).divide(dxpm).subtract(sig.multiply(u[i - 1])).divide(p);
            }

            double qn = 0;
            T un = zero;
            if (ypn.getReal() < 1e+30) {
                final T dx12 = x[n - 1].subtract(x[n - 2]);
                final T dy12 = y[n - 1].subtract(y[n - 2]);
                qn = 0.5;
                un = dx12.reciprocal().multiply(3.0).multiply(ypn.subtract(dy12.divide(dx12)));
            }

            y2[n - 1] = un.subtract(u[n - 2].multiply(qn)).divide(y2[n - 2].multiply(qn).add(1.0));
            for (int k = n - 2; k >= 0; k--) {
                y2[k] = y2[k].multiply(y2[k + 1]).add(u[k]);
            }

            return y2;

        }

        /** Calculate Temperature and Density Profiles for lower atmosphere.
         * @param alt altitude
         * @param d0 density
         * @param xm mixed density
         * @return temperature or density profile
         */
        private T densm(final T alt, final T d0, final double xm) {

            T densm = d0;

            // stratosphere/mesosphere temperature
            int mn = ZN2.length;
            T z = (alt.getReal() > ZN2[mn - 1]) ? alt : zero.add(ZN2[mn - 1]);

            double z1 = ZN2[0];
            double z2 = ZN2[mn - 1];
            T t1 = meso_tn2[0];
            T t2 = meso_tn2[mn - 1];
            T zg  = zeta(z, z1);
            T zgdif = zeta(zero.add(z2), z1);

            /* set up spline nodes */
            T[] xs = MathArrays.buildArray(field, mn);
            T[] ys = MathArrays.buildArray(field, mn);
            for (int k = 0; k < mn; k++) {
                xs[k] = zeta(zero.add(ZN2[k]), z1).divide(zgdif);
                ys[k] = meso_tn2[k].reciprocal();
            }
            final T qSM = rlat.add(z2).divide(rlat.add(z1));
            T yd1 = meso_tgn2[0].negate().divide(t1.multiply(t1)).multiply(zgdif);
            T yd2 = meso_tgn2[1].negate().divide(t2.multiply(t2)).multiply(zgdif).multiply(qSM).multiply(qSM);

            /* calculate spline coefficients */
            T[] y2out = spline(xs, ys, yd1, yd2);
            T x = zg.divide(zgdif);
            T y = splint(xs, ys, y2out, x);

            /* temperature at altitude */
            T tz = y.reciprocal();

            if (xm != 0.0) {
                /* calculate stratosphere / mesospehere density */
                final T glb  = galt(zero.add(z1));
                final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);

                /* Integrate temperature profile */
                final T yi = splini(xs, ys, y2out, x);
                final T expl = min(50., gamm.multiply(yi));

                /* Density at altitude */
                densm = densm.multiply(t1.divide(tz).multiply(expl.negate().exp()));
            }

            if (alt.getReal() > ZN3[0]) {
                return (xm == 0.0) ? tz : densm;
            }

            // troposhere/stratosphere temperature
            z = alt;
            mn = ZN3.length;
            z1 = ZN3[0];
            z2 = ZN3[mn - 1];
            t1 = meso_tn3[0];
            t2 = meso_tn3[mn - 1];
            zg = zeta(z, z1);
            zgdif = zeta(zero.add(z2), z1);

            /* set up spline nodes */
            xs = MathArrays.buildArray(field, mn);
            ys = MathArrays.buildArray(field, mn);
            for (int k = 0; k < mn; k++) {
                xs[k] = zeta(zero.add(ZN3[k]), z1).divide(zgdif);
                ys[k] = meso_tn3[k].reciprocal();
            }
            final T qTS = rlat.add(z2) .divide(rlat.add(z1));
            yd1 = meso_tgn3[0].negate().divide(t1.multiply(t1)).multiply(zgdif);
            yd2 = meso_tgn3[1].negate().divide(t2.multiply(t2)).multiply(zgdif).multiply(qTS).multiply(qTS);

            /* calculate spline coefficients */
            y2out = spline(xs, ys, yd1, yd2);
            x = zg.divide(zgdif);
            y = splint(xs, ys, y2out, x);

            /* temperature at altitude */
            tz = y.reciprocal();

            if (xm != 0.0) {
                /* calculate tropospheric / stratosphere density */
                final T glb = galt(zero.add(z1));
                final T gamm = glb.multiply(zgdif).multiply(xm / R_GAS);

                /* Integrate temperature profile */
                final T yi = splini(xs, ys, y2out, x);
                final T expl = min(50., gamm.multiply(yi));

                /* Density at altitude */
                densm = densm.multiply(t1.divide(tz).multiply(expl.negate().exp()));
            }

            return (xm == 0.0) ? tz : densm;
        }

        /** Calculate temperature and density profiles according to new lower thermo polynomial.
         * @param alt altitude
         * @param dlb density at lower boundary
         * @param tinf exospheric temperature
         * @param tlb temperature at lower boundary
         * @param xm species molecular weight
         * @param alpha thermal diffusion coefficient
         * @param zlb altitude of the lower boundary
         * @param s2 slope
         * @return temperature or density profile
         */
        private T densu(final T alt, final T dlb, final T tinf,
                        final T tlb, final double xm,  final double alpha,
                        final double zlb, final T s2) {
            /* joining altitudes of Bates and spline */
            T z = (alt.getReal() > ZN1[0]) ? alt : zero.add(ZN1[0]);

            /* geopotential altitude difference from ZLB */
            final T zg2 = zeta(z, zlb);

            /* Bates temperature */
            final T tt = tinf.subtract(tinf.subtract(tlb).multiply(s2.negate().multiply(zg2).exp()));
            final T ta = tt;
            T tz = tt;

            final int mn = ZN1.length;
            final T[] xs = MathArrays.buildArray(field, mn);
            final T[] ys = MathArrays.buildArray(field, mn);
            T x = zero;
            T[] y2out =  MathArrays.buildArray(field, mn);
            T zgdif = zero;
            if (alt.getReal() < ZN1[0]) {
                /* calculate temperature below ZA
                 * temperature gradient at ZA from Bates profile */
                final T p = rlat.add(zlb).divide(rlat.add(ZN1[0]));
                final T dta = tinf.subtract(ta).multiply(s2).multiply(p.multiply(p));
                meso_tgn1[0] = dta;
                meso_tn1[0] = ta;
                final T tzn1mn1 = zero.add(ZN1[mn - 1]);
                z = (alt.getReal() > ZN1[mn - 1]) ? alt : tzn1mn1;

                final T t1 = meso_tn1[0];
                final T t2 = meso_tn1[mn - 1];
                /* geopotental difference from z1 */
                final T zg = zeta(z, ZN1[0]);
                zgdif = zeta(tzn1mn1, ZN1[0]);
                /* set up spline nodes */
                for (int k = 0; k < mn; k++) {
                    xs[k] = zeta(zero.add(ZN1[k]), ZN1[0]).divide(zgdif);
                    ys[k] =  meso_tn1[k].reciprocal();
                }
                /* end node derivatives */
                final T q   = rlat.add(ZN1[mn - 1]).divide(rlat.add(ZN1[0]));
                final T yd1 = meso_tgn1[0].negate().divide(t1.multiply(t1)).multiply(zgdif);
                final T yd2 = meso_tgn1[1].negate().divide(t2.multiply(t2)).multiply(zgdif).multiply(q.multiply(q));
                /* calculate spline coefficients */
                y2out = spline(xs, ys, yd1, yd2);
                x = zg.divide(zgdif);
                final T y = splint(xs, ys, y2out, x);
                /* temperature at altitude */
                tz = y.reciprocal();
            }

            if (xm == 0) {
                return tz;
            }

            /* calculate density above za */
            T glb   = galt(zero.add(zlb));
            T gamma = glb.divide(s2.multiply(tinf)).multiply(xm / R_GAS);
            T expl = tt.getReal() <= 0 ?
                     zero.add(50) :
                     min(50.0, s2.negate().multiply(gamma).multiply(zg2).exp());
            T densu = dlb.multiply(tlb.divide(tt).pow(gamma.add(alpha + 1))).multiply(expl);

            /* calculate density below za */
            if (alt.getReal() < ZN1[0]) {
                glb   = galt(zero.add(ZN1[0]));
                gamma = glb.multiply(zgdif).multiply(xm / R_GAS);
                /* integrate spline temperatures */
                expl = tz.getReal() <= 0 ?
                       zero.add(50.0) :
                       min(50.0, gamma.multiply(splini(xs, ys, y2out, x)));
                /* correct density at altitude */
                densu = densu.multiply(meso_tn1[0].divide(tz).pow(alpha + 1).multiply(expl.negate().exp()));
            }

            /* Return density at altitude */
            return densu;
        }

        /** Compute min of two values, one double and one field element.
         * @param d double value
         * @param f field element
         * @return min value
         */
        private T min(final double d, final T f) {
            return (f.getReal() > d) ? zero.add(d) : f;
        }

        /** Calculate gravity at altitude.
         * @param alt altitude (km)
         * @return gravity at altitude (cm/s2)
         */
        private T galt(final T alt) {
            final T r = alt.divide(rlat).add(1);
            return glat.divide(r.multiply(r));
        }

        /** Calculate zeta function.
         * @param zz zz value
         * @param zl zl value
         * @return value of zeta function
         */
        private T zeta(final T zz, final double zl) {
            return zz.subtract(zl).multiply(rlat.add(zl)).divide(rlat.add(zz));
        }

    }

}