JB2008.java

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

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.MathArrays;
  25. import org.hipparchus.util.MathUtils;
  26. import org.orekit.bodies.BodyShape;
  27. import org.orekit.bodies.FieldGeodeticPoint;
  28. import org.orekit.bodies.GeodeticPoint;
  29. import org.orekit.errors.OrekitException;
  30. import org.orekit.errors.OrekitMessages;
  31. import org.orekit.frames.Frame;
  32. import org.orekit.time.AbsoluteDate;
  33. import org.orekit.time.FieldAbsoluteDate;
  34. import org.orekit.utils.Constants;
  35. import org.orekit.utils.PVCoordinatesProvider;

  37. /** This is the realization of the Jacchia-Bowman 2008 atmospheric model.
  38.  * <p>
  39.  * It is described in the paper:<br>
  40.  * <a href="http://sol.spacenvironment.net/~JB2008/pubs/AIAA_2008-6438_JB2008_Model.pdf">A
  41.  * New Empirical Thermospheric Density Model JB2008 Using New Solar Indices</a><br>
  42.  * <i>Bruce R. Bowman & al.</i><br>
  43.  * AIAA 2008-6438<br>
  44.  * </p>
  45.  * <p>
  46.  * Two computation methods are proposed to the user:
  47.  * <ul>
  48.  * <li>one OREKIT independent and compliant with initial FORTRAN routine entry values:
  49.  *     {@link #getDensity(double, double, double, double, double, double, double, double, double, double, double, double, double, double, double)}. </li>
  50.  * <li>one compliant with OREKIT Atmosphere interface, necessary to the
  51.  *     {@link org.orekit.forces.drag.DragForce drag force model} computation.</li>
  52.  * </ul>
  53.  * </p>
  54.  * <p>
  55.  * This model provides dense output for all altitudes and positions. Output data are :
  56.  * <ul>
  57.  * <li>Exospheric Temperature above Input Position (deg K)</li>
  58.  * <li>Temperature at Input Position (deg K)</li>
  59.  * <li>Total Mass-Density at Input Position (kg/m³)</li>
  60.  * </ul>
  61.  * </p>
  62.  * <p>
  63.  * The model needs geographical and time information to compute general values,
  64.  * but also needs space weather data : mean and daily solar flux, retrieved through
  65.  * different indices, and planetary geomagnetic indices.<br>
  66.  * More information on these indices can be found on the <a
  67.  * href="http://sol.spacenvironment.net/~JB2008/indices.html">
  68.  * official JB2008 website.</a>
  69.  * </p>
  70.  *
  71.  * @author Bruce R Bowman (HQ AFSPC, Space Analysis Division), 2008: FORTRAN routine
  72.  * @author Pascal Parraud (java translation)
  73.  */
  74. public class JB2008 implements Atmosphere {

  76.     /** Serializable UID. */
  77.     private static final long serialVersionUID = -4201270765122160831L;

  79.     /** Minimum altitude (m) for JB2008 use. */
  80.     private static final double ALT_MIN = 90000.;

  82.     /** Earth radius (km). */
  83.     private static final double EARTH_RADIUS = 6356.766;

  85.     /** Natural logarithm of 10.0. */
  86.     private static final double LOG10  = FastMath.log(10.);

  88.     /** Avogadro's number in mks units (molecules/kmol). */
  89.     private static final double AVOGAD = 6.02257e26;

  91.     /** Universal gas-constant in mks units (joules/K/kmol). */
  92.     private static final double RSTAR  = 8.31432;

  94.     /** The alpha are the thermal diffusion coefficients in Equation (6). */
  95.     private static final double[] ALPHA = {
  96.         0, 0, 0, 0, -0.38
  97.     };

  99.     /** Molecular weights in order: N2, O2, O, Ar, He and H. */
  100.     private static final double[] AMW = {
  101.         28.0134, 31.9988, 15.9994, 39.9480, 4.0026, 1.00797
  102.     };

  104.     /** The FRAC are the assumed sea-level volume fractions in order: N2, O2, Ar, and He. */
  105.     private static final double[] FRAC = {
  106.         0.78110, 0.20955, 9.3400e-3, 1.2890e-5
  107.     };

  109.     /** Value used to establish height step sizes in the regime 90km to 105km. */
  110.     private static final double R1 = 0.010;

  112.     /** Value used to establish height step sizes in the regime 105km to 500km. */
  113.     private static final double R2 = 0.025;

  115.     /** Value used to establish height step sizes in the regime above 500km. */
  116.     private static final double R3 = 0.075;

  118.     /** Weights for the Newton-Cotes five-points quadrature formula. */
  119.     private static final double[] WT = {
  120.         0.311111111111111, 1.422222222222222, 0.533333333333333, 1.422222222222222, 0.311111111111111
  121.     };

  123.     /** Coefficients for high altitude density correction. */
  124.     private static final double[] CHT = {
  125.         0.22, -0.20e-02, 0.115e-02, -0.211e-05
  126.     };

  128.     /** FZ global model values (1997-2006 fit).  */
  129.     private static final double[] FZM = {
  130.         0.2689e+00, -0.1176e-01, 0.2782e-01, -0.2782e-01, 0.3470e-03
  131.     };

  133.     /** GT global model values (1997-2006 fit). */
  134.     private static final double[] GTM = {
  135.         -0.3633e+00, 0.8506e-01,  0.2401e+00, -0.1897e+00, -0.2554e+00,
  136.         -0.1790e-01, 0.5650e-03, -0.6407e-03, -0.3418e-02, -0.1252e-02
  137.     };

  139.     /** Mbar polynomial coeffcients. */
  140.     private static final double[] CMB = {
  141.         28.15204, -8.5586e-2, +1.2840e-4, -1.0056e-5, -1.0210e-5, +1.5044e-6, +9.9826e-8
  142.     };

  144.     /** DTC relative data. */
  145.     private static final double[] BDTC = {
  146.         -0.457512297e+01, -0.512114909e+01, -0.693003609e+02,
  147.         0.203716701e+03,  0.703316291e+03, -0.194349234e+04,
  148.         0.110651308e+04, -0.174378996e+03,  0.188594601e+04,
  149.         -0.709371517e+04,  0.922454523e+04, -0.384508073e+04,
  150.         -0.645841789e+01,  0.409703319e+02, -0.482006560e+03,
  151.         0.181870931e+04, -0.237389204e+04,  0.996703815e+03,
  152.         0.361416936e+02
  153.     };

  155.     /** DTC relative data. */
  156.     private static final double[] CDTC = {
  157.         -0.155986211e+02, -0.512114909e+01, -0.693003609e+02,
  158.         0.203716701e+03,  0.703316291e+03, -0.194349234e+04,
  159.         0.110651308e+04, -0.220835117e+03,  0.143256989e+04,
  160.         -0.318481844e+04,  0.328981513e+04, -0.135332119e+04,
  161.         0.199956489e+02, -0.127093998e+02,  0.212825156e+02,
  162.         -0.275555432e+01,  0.110234982e+02,  0.148881951e+03,
  163.         -0.751640284e+03,  0.637876542e+03,  0.127093998e+02,
  164.         -0.212825156e+02,  0.275555432e+01
  165.     };

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

  170.     /** External data container. */
  171.     private JB2008InputParameters inputParams;

  173.     /** Earth body shape. */
  174.     private BodyShape earth;

  176.     /** Constructor with space environment information for internal computation.
  177.      * @param parameters the solar and magnetic activity data
  178.      * @param sun the sun position
  179.      * @param earth the earth body shape
  180.      */
  181.     public JB2008(final JB2008InputParameters parameters,
  182.                   final PVCoordinatesProvider sun, final BodyShape earth) {
  183.         this.earth = earth;
  184.         this.sun = sun;
  185.         this.inputParams = parameters;
  186.     }

  188.     /** {@inheritDoc} */
  189.     public Frame getFrame() {
  190.         return earth.getBodyFrame();
  191.     }

  193.     /** Get the local density with initial entries.
  194.      * @param dateMJD date and time, in modified julian days and fraction
  195.      * @param sunRA Right Ascension of Sun (radians)
  196.      * @param sunDecli Declination of Sun (radians)
  197.      * @param satLon Right Ascension of position (radians)
  198.      * @param satLat Geocentric latitude of position (radians)
  199.      * @param satAlt Height of position (m)
  200.      * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
  201.      *        (Tabular time 1.0 day earlier)
  202.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
  203.      *        (Tabular time 1.0 day earlier)
  204.      * @param s10 EUV index (26-34 nm) scaled to F10<br>
  205.      *        (Tabular time 1 day earlier)
  206.      * @param s10B UV 81-day averaged centered index
  207.      *        (Tabular time 1 day earlier)
  208.      * @param xm10 MG2 index scaled to F10<br>
  209.      *        (Tabular time 2.0 days earlier)
  210.      * @param xm10B MG2 81-day ave. centered index<br>
  211.      *        (Tabular time 2.0 days earlier)
  212.      * @param y10 Solar X-Ray & Lya index scaled to F10<br>
  213.      *        (Tabular time 5.0 days earlier)
  214.      * @param y10B Solar X-Ray & Lya 81-day ave. centered index<br>
  215.      *        (Tabular time 5.0 days earlier)
  216.      * @param dstdtc Temperature change computed from Dst index
  217.      * @return total mass-Density at input position (kg/m³)
  218.      * @exception OrekitException if altitude is below 90 km
  219.      */
  220.     public double getDensity(final double dateMJD, final double sunRA, final double sunDecli,
  221.                              final double satLon, final double satLat, final double satAlt,
  222.                              final double f10, final double f10B, final double s10,
  223.                              final double s10B, final double xm10, final double xm10B,
  224.                              final double y10, final double y10B, final double dstdtc)
  225.         throws OrekitException {

  227.         if (satAlt < ALT_MIN) {
  228.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, satAlt, ALT_MIN);
  229.         }

  231.         final double altKm = satAlt / 1000.0;

  233.         // Equation (14)
  234.         final double fn  = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
  235.         final double fsb = f10B * fn + s10B * (1. - fn);
  236.         final double tsubc = 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) +
  237.                              0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);

  239.         // Equation (15)
  240.         final double eta   = 0.5 * FastMath.abs(satLat - sunDecli);
  241.         final double theta = 0.5 * FastMath.abs(satLat + sunDecli);

  243.         // Equation (16)
  244.         final double h   = satLon - sunRA;
  245.         final double tau = h - 0.64577182 + 0.10471976 * FastMath.sin(h + 0.75049158);
  246.         double solarTime = FastMath.toDegrees(h + FastMath.PI) / 15.0;
  247.         while (solarTime >= 24) {
  248.             solarTime -= 24.;
  249.         }
  250.         while (solarTime < 0) {
  251.             solarTime += 24.;
  252.         }

  254.         // Equation (17)
  255.         final double cosEta  = FastMath.pow(FastMath.cos(eta), 2.5);
  256.         final double sinTeta = FastMath.pow(FastMath.sin(theta), 2.5);
  257.         final double cosTau  = FastMath.abs(FastMath.cos(0.5 * tau));
  258.         final double df = sinTeta + (cosEta - sinTeta) * cosTau * cosTau * cosTau;
  259.         final double tsubl = tsubc * (1. + 0.31 * df);

  261.         // Compute correction to dTc for local solar time and lat correction
  262.         final double dtclst = dTc(f10, solarTime, satLat, altKm);

  264.         // Compute the local exospheric temperature.
  265.         // Add geomagnetic storm effect from input dTc value
  266.         final double[] temp = new double[2];
  267.         temp[0] = tsubl + dstdtc;
  268.         final double tinf = temp[0] + dtclst;

  270.         // Equation (9)
  271.         final double tsubx = 444.3807 + 0.02385 * tinf - 392.8292 * FastMath.exp(-0.0021357 * tinf);

  273.         // Equation (11)
  274.         final double gsubx = 0.054285714 * (tsubx - 183.);

  276.         // The TC array will be an argument in the call to localTemp,
  277.         // which evaluates Eq. (10) or (13)
  278.         final double[] tc = new double[4];
  279.         tc[0] = tsubx;
  280.         tc[1] = gsubx;
  281.         // A of Equation (13)
  282.         tc[2] = (tinf - tsubx) * 2. / FastMath.PI;
  283.         tc[3] = gsubx / tc[2];

  285.         // Equation (5)
  286.         final double z1 = 90.;
  287.         final double z2 = FastMath.min(altKm, 105.0);
  288.         double al = FastMath.log(z2 / z1);
  289.         int n = 1 + (int) FastMath.floor(al / R1);
  290.         double zr = FastMath.exp(al / n);
  291.         final double mb1 = mBar(z1);
  292.         final double tloc1 = localTemp(z1, tc);
  293.         double zend  = z1;
  294.         double sub2  = 0.;
  295.         double ain   = mb1 * gravity(z1) / tloc1;
  296.         double mb2   = 0;
  297.         double tloc2 = 0;
  298.         double z     = 0;
  299.         double gravl = 0;

  301.         for (int i = 0; i < n; ++i) {
  302.             z = zend;
  303.             zend = zr * z;
  304.             final double dz = 0.25 * (zend - z);
  305.             double sum1 = WT[0] * ain;
  306.             for (int j = 1; j < 5; ++j) {
  307.                 z += dz;
  308.                 mb2   = mBar(z);
  309.                 tloc2 = localTemp(z, tc);
  310.                 gravl = gravity(z);
  311.                 ain   = mb2 * gravl / tloc2;
  312.                 sum1 += WT[j] * ain;
  313.             }
  314.             sub2 += dz * sum1;
  315.         }

  317.         double rho = 3.46e-6 * mb2 * tloc1 / FastMath.exp(sub2 / RSTAR) / (mb1 * tloc2);

  319.         // Equation (2)
  320.         final double anm = AVOGAD * rho;
  321.         final double an  = anm / mb2;

  323.         // Equation (3)
  324.         double fact2  = anm / 28.960;
  325.         final double[] aln = new double[6];
  326.         aln[0] = FastMath.log(FRAC[0] * fact2);
  327.         aln[3] = FastMath.log(FRAC[2] * fact2);
  328.         aln[4] = FastMath.log(FRAC[3] * fact2);
  329.         // Equation (4)
  330.         aln[1] = FastMath.log(fact2 * (1. + FRAC[1]) - an);
  331.         aln[2] = FastMath.log(2. * (an - fact2));

  333.         if (altKm <= 105.0) {
  334.             temp[1] = tloc2;
  335.             // Put in negligible hydrogen for use in DO-LOOP 13
  336.             aln[5] = aln[4] - 25.0;
  337.         } else {
  338.             // Equation (6)
  339.             al   = FastMath.log(FastMath.min(altKm, 500.0) / z);
  340.             n    = 1 + (int) FastMath.floor(al / R2);
  341.             zr   = FastMath.exp(al / n);
  342.             sub2 = 0.;
  343.             ain  = gravl / tloc2;

  345.             double tloc3 = 0;
  346.             for (int i = 0; i < n; ++i) {
  347.                 z = zend;
  348.                 zend = zr * z;
  349.                 final double dz = 0.25 * (zend - z);
  350.                 double sum1 = WT[0] * ain;
  351.                 for (int j = 1; j < 5; ++j) {
  352.                     z += dz;
  353.                     tloc3 = localTemp(z, tc);
  354.                     gravl = gravity(z);
  355.                     ain   = gravl / tloc3;
  356.                     sum1 += WT[j] * ain;
  357.                 }
  358.                 sub2 += dz * sum1;
  359.             }

  361.             al = FastMath.log(FastMath.max(altKm, 500.0) / z);
  362.             final double r = (altKm > 500.0) ? R3 : R2;
  363.             n = 1 + (int) FastMath.floor(al / r);
  364.             zr = FastMath.exp(al / n);
  365.             double sum3 = 0.;
  366.             double tloc4 = 0;
  367.             for (int i = 0; i < n; ++i) {
  368.                 z = zend;
  369.                 zend = zr * z;
  370.                 final double dz = 0.25 * (zend - z);
  371.                 double sum1 = WT[0] * ain;
  372.                 for (int j = 1; j < 5; ++j) {
  373.                     z += dz;
  374.                     tloc4 = localTemp(z, tc);
  375.                     gravl = gravity(z);
  376.                     ain   = gravl / tloc4;
  377.                     sum1 += WT[j] * ain;
  378.                 }
  379.                 sum3 += dz * sum1;
  380.             }
  381.             final double altr;
  382.             final double hSign;
  383.             if (altKm <= 500.) {
  384.                 temp[1] = tloc3;
  385.                 altr = FastMath.log(tloc3 / tloc2);
  386.                 fact2 = sub2 / RSTAR;
  387.                 hSign = 1.0;
  388.             } else {
  389.                 temp[1] = tloc4;
  390.                 altr = FastMath.log(tloc4 / tloc2);
  391.                 fact2 = (sub2 + sum3) / RSTAR;
  392.                 hSign = -1.0;
  393.             }
  394.             for (int i = 0; i < 5; ++i) {
  395.                 aln[i] = aln[i] - (1.0 + ALPHA[i]) * altr - fact2 * AMW[i];
  396.             }

  398.             // Equation (7)
  399.             final double al10t5 = FastMath.log10(tinf);
  400.             final double alnh5 = (5.5 * al10t5 - 39.40) * al10t5 + 73.13;
  401.             aln[5] = LOG10 * (alnh5 + 6.) + hSign * (FastMath.log(tloc4 / tloc3) + sum3 * AMW[5] / RSTAR);
  402.         }

  404.         // Equation (24) - J70 Seasonal-Latitudinal Variation
  405.         final double capPhi = ((dateMJD - 36204.0) / 365.2422) % 1;
  406.         final int signum = (satLat >= 0.) ? 1 : -1;
  407.         final double sinLat = FastMath.sin(satLat);
  408.         final double hm90  = altKm - 90.;
  409.         final double dlrsl = 0.02 * hm90 * FastMath.exp(-0.045 * hm90) *
  410.                              signum * FastMath.sin(MathUtils.TWO_PI * capPhi + 1.72) * sinLat * sinLat;

  412.         // Equation (23) - Computes the semiannual variation
  413.         double dlrsa = 0;
  414.         if (z < 2000.0) {
  415.             // Use new semiannual model dLog(rho)
  416.             dlrsa = semian08(dayOfYear(dateMJD), altKm, f10B, s10B, xm10B);
  417.         }

  419.         // Sum the delta-log-rhos and apply to the number densities.
  420.         // In CIRA72 the following equation contains an actual sum,
  421.         // namely DLR = LOG10 * (DLRGM + DLRSA + DLRSL)
  422.         // However, for Jacchia 70, there is no DLRGM or DLRSA.
  423.         final double dlr = LOG10 * (dlrsl + dlrsa);
  424.         for (int i = 0; i < 6; ++i) {
  425.             aln[i] += dlr;
  426.         }

  428.         // Compute mass-density and mean-molecular-weight and
  429.         // convert number density logs from natural to common.
  430.         double sumnm = 0.0;
  431.         for (int i = 0; i < 6; ++i) {
  432.             sumnm += FastMath.exp(aln[i]) * AMW[i];
  433.         }
  434.         rho = sumnm / AVOGAD;

  436.         // Compute the high altitude exospheric density correction factor
  437.         double fex = 1.0;
  438.         if ((altKm >= 1000.0) && (altKm < 1500.0)) {
  439.             final double zeta = (altKm - 1000.) * 0.002;
  440.             final double f15c = CHT[0] + CHT[1] * f10B + (CHT[2] + CHT[3] * f10B) * 1500.0;
  441.             final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
  442.             final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
  443.             final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
  444.             fex += zeta * zeta * (fex2 + fex3 * zeta);
  445.         } else if (altKm >= 1500.0) {
  446.             fex = CHT[0] + CHT[1] * f10B + CHT[2] * altKm + CHT[3] * f10B * altKm;
  447.         }

  449.         // Apply the exospheric density correction factor.
  450.         rho *= fex;

  452.         return rho;
  453.     }

  455.     /** Get the local density with initial entries.
  456.      * @param dateMJD date and time, in modified julian days and fraction
  457.      * @param sunRA Right Ascension of Sun (radians)
  458.      * @param sunDecli Declination of Sun (radians)
  459.      * @param satLon Right Ascension of position (radians)
  460.      * @param satLat Geocentric latitude of position (radians)
  461.      * @param satAlt Height of position (m)
  462.      * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
  463.      *        (Tabular time 1.0 day earlier)
  464.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
  465.      *        (Tabular time 1.0 day earlier)
  466.      * @param s10 EUV index (26-34 nm) scaled to F10<br>
  467.      *        (Tabular time 1 day earlier)
  468.      * @param s10B UV 81-day averaged centered index
  469.      *        (Tabular time 1 day earlier)
  470.      * @param xm10 MG2 index scaled to F10<br>
  471.      *        (Tabular time 2.0 days earlier)
  472.      * @param xm10B MG2 81-day ave. centered index<br>
  473.      *        (Tabular time 2.0 days earlier)
  474.      * @param y10 Solar X-Ray & Lya index scaled to F10<br>
  475.      *        (Tabular time 5.0 days earlier)
  476.      * @param y10B Solar X-Ray & Lya 81-day ave. centered index<br>
  477.      *        (Tabular time 5.0 days earlier)
  478.      * @param dstdtc Temperature change computed from Dst index
  479.      * @param <T> type fo the field elements
  480.      * @return total mass-Density at input position (kg/m³)
  481.      * @exception OrekitException if altitude is below 90 km
  482.      */
  483.     public <T extends RealFieldElement<T>> T getDensity(final T dateMJD, final T sunRA, final T sunDecli,
  484.                                                         final T satLon, final T satLat, final T satAlt,
  485.                                                         final double f10, final double f10B, final double s10,
  486.                                                         final double s10B, final double xm10, final double xm10B,
  487.                                                         final double y10, final double y10B, final double dstdtc)
  488.         throws OrekitException {

  490.         if (satAlt.getReal() < ALT_MIN) {
  491.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
  492.                                       satAlt.getReal(), ALT_MIN);
  493.         }

  495.         final Field<T> field  = satAlt.getField();
  496.         final T altKm = satAlt.divide(1000.0);

  498.         // Equation (14)
  499.         final double fn  = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
  500.         final double fsb = f10B * fn + s10B * (1. - fn);
  501.         final double tsubc = 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) +
  502.                              0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);

  504.         // Equation (15)
  505.         final T eta   = satLat.subtract(sunDecli).abs().multiply(0.5);
  506.         final T theta = satLat.add(sunDecli).abs().multiply(0.5);

  508.         // Equation (16)
  509.         final T h   = satLon.subtract(sunRA);
  510.         final T tau = h.subtract(0.64577182).add(h.add(0.75049158).sin().multiply(0.10471976));
  511.         T solarTime = h.add(FastMath.PI).multiply(12.0 / FastMath.PI);
  512.         while (solarTime.getReal() >= 24) {
  513.             solarTime = solarTime.subtract(24);
  514.         }
  515.         while (solarTime.getReal() < 0) {
  516.             solarTime = solarTime.add(24);
  517.         }

  519.         // Equation (17)
  520.         final T cos     = eta.cos();
  521.         final T cosEta  = cos.multiply(cos).multiply(cos.sqrt());
  522.         final T sin     = theta.sin();
  523.         final T sinTeta = sin.multiply(sin).multiply(sin.sqrt());
  524.         final T cosTau  = tau.multiply(0.5).cos().abs();
  525.         final T df      = sinTeta.add(cosEta.subtract(sinTeta).multiply(cosTau).multiply(cosTau).multiply(cosTau));
  526.         final T tsubl   = df.multiply(0.31).add(1).multiply(tsubc);

  528.         // Compute correction to dTc for local solar time and lat correction
  529.         final T dtclst = dTc(f10, solarTime, satLat, altKm);

  531.         // Compute the local exospheric temperature.
  532.         // Add geomagnetic storm effect from input dTc value
  533.         final T[] temp = MathArrays.buildArray(field, 2);
  534.         temp[0] = tsubl.add(dstdtc);
  535.         final T tinf = temp[0].add(dtclst);

  537.         // Equation (9)
  538.         final T tsubx = tinf.multiply(0.02385).add(444.3807).subtract(tinf.multiply(-0.0021357).exp().multiply(392.8292));

  540.         // Equation (11)
  541.         final T gsubx = tsubx.subtract(183.).multiply(0.054285714);

  543.         // The TC array will be an argument in the call to localTemp,
  544.         // which evaluates Eq. (10) or (13)
  545.         final T[] tc = MathArrays.buildArray(field, 4);
  546.         tc[0] = tsubx;
  547.         tc[1] = gsubx;
  548.         // A of Equation (13)
  549.         tc[2] = tinf.subtract(tsubx).multiply(2. / FastMath.PI);
  550.         tc[3] = gsubx.divide(tc[2]);

  552.         // Equation (5)
  553.         final T z1 = field.getZero().add(90.);
  554.         final T z2 = min(105.0, altKm);
  555.         T al = z2.divide(z1).log();
  556.         int n = 1 + (int) FastMath.floor(al.getReal() / R1);
  557.         T zr = al.divide(n).exp();
  558.         final T mb1 = mBar(z1);
  559.         final T tloc1 = localTemp(z1, tc);
  560.         T zend  = z1;
  561.         T sub2  = field.getZero();
  562.         T ain   = mb1.multiply(gravity(z1)).divide(tloc1);
  563.         T mb2   = field.getZero();
  564.         T tloc2 = field.getZero();
  565.         T z     = field.getZero();
  566.         T gravl = field.getZero();
  567.         for (int i = 0; i < n; ++i) {
  568.             z = zend;
  569.             zend = zr.multiply(z);
  570.             final T dz = zend.subtract(z).multiply(0.25);
  571.             T sum1 = ain.multiply(WT[0]);
  572.             for (int j = 1; j < 5; ++j) {
  573.                 z = z.add(dz);
  574.                 mb2   = mBar(z);
  575.                 tloc2 = localTemp(z, tc);
  576.                 gravl = gravity(z);
  577.                 ain   = mb2.multiply(gravl).divide(tloc2);
  578.                 sum1  = sum1.add(ain.multiply(WT[j]));
  579.             }
  580.             sub2 = sub2.add(dz.multiply(sum1));
  581.         }

  583.         T rho = mb2.multiply(3.46e-6).multiply(tloc1).divide(sub2.divide(RSTAR).exp().multiply(mb1.multiply(tloc2)));

  585.         // Equation (2)
  586.         final T anm = rho.multiply(AVOGAD);
  587.         final T an  = anm.divide(mb2);

  589.         // Equation (3)
  590.         T fact2  = anm.divide(28.960);
  591.         final T[] aln = MathArrays.buildArray(field, 6);
  592.         aln[0] = fact2.multiply(FRAC[0]).log();
  593.         aln[3] = fact2.multiply(FRAC[2]).log();
  594.         aln[4] = fact2.multiply(FRAC[3]).log();
  595.         // Equation (4)
  596.         aln[1] = fact2.multiply(1. + FRAC[1]).subtract(an).log();
  597.         aln[2] = an.subtract(fact2).multiply(2).log();

  599.         if (altKm.getReal() <= 105.0) {
  600.             temp[1] = tloc2;
  601.             // Put in negligible hydrogen for use in DO-LOOP 13
  602.             aln[5] = aln[4].subtract(25.0);
  603.         } else {
  604.             // Equation (6)
  605.             al   = min(500.0, altKm).divide(z).log();
  606.             n    = 1 + (int) FastMath.floor(al.getReal() / R2);
  607.             zr   = al.divide(n).exp();
  608.             sub2 = field.getZero();
  609.             ain  = gravl.divide(tloc2);

  611.             T tloc3 = field.getZero();
  612.             for (int i = 0; i < n; ++i) {
  613.                 z = zend;
  614.                 zend = zr.multiply(z);
  615.                 final T dz = zend.subtract(z).multiply(0.25);
  616.                 T sum1 = ain.multiply(WT[0]);
  617.                 for (int j = 1; j < 5; ++j) {
  618.                     z = z.add(dz);
  619.                     tloc3 = localTemp(z, tc);
  620.                     gravl = gravity(z);
  621.                     ain   = gravl.divide(tloc3);
  622.                     sum1  = sum1.add(ain.multiply(WT[j]));
  623.                 }
  624.                 sub2 = sub2.add(dz.multiply(sum1));
  625.             }

  627.             al = max(500.0, altKm).divide(z).log();
  628.             final double r = (altKm.getReal() > 500.0) ? R3 : R2;
  629.             n = 1 + (int) FastMath.floor(al.getReal() / r);
  630.             zr = al.divide(n).exp();
  631.             T sum3 = field.getZero();
  632.             T tloc4 = field.getZero();
  633.             for (int i = 0; i < n; ++i) {
  634.                 z = zend;
  635.                 zend = zr.multiply(z);
  636.                 final T dz = zend.subtract(z).multiply(0.25);
  637.                 T sum1 = ain.multiply(WT[0]);
  638.                 for (int j = 1; j < 5; ++j) {
  639.                     z = z.add(dz);
  640.                     tloc4 = localTemp(z, tc);
  641.                     gravl = gravity(z);
  642.                     ain   = gravl.divide(tloc4);
  643.                     sum1  = sum1.add(ain.multiply(WT[j]));
  644.                 }
  645.                 sum3 = sum3.add(dz.multiply(sum1));
  646.             }
  647.             final T altr;
  648.             final double hSign;
  649.             if (altKm.getReal() <= 500.) {
  650.                 temp[1] = tloc3;
  651.                 altr = tloc3.divide(tloc2).log();
  652.                 fact2 = sub2.divide(RSTAR);
  653.                 hSign = 1.0;
  654.             } else {
  655.                 temp[1] = tloc4;
  656.                 altr = tloc4.divide(tloc2).log();
  657.                 fact2 = sub2.add(sum3).divide(RSTAR);
  658.                 hSign = -1.0;
  659.             }
  660.             for (int i = 0; i < 5; ++i) {
  661.                 aln[i] = aln[i].subtract(altr.multiply(1.0 + ALPHA[i])).subtract(fact2.multiply(AMW[i]));
  662.             }

  664.             // Equation (7)
  665.             final T al10t5 = tinf.log10();
  666.             final T alnh5 = al10t5.multiply(5.5).subtract(39.40).multiply(al10t5).add(73.13);
  667.             aln[5] = alnh5.add(6.).multiply(LOG10).
  668.                      add(tloc4.divide(tloc3).log().add(sum3.multiply(AMW[5] / RSTAR)).multiply(hSign));
  669.         }

  671.         // Equation (24) - J70 Seasonal-Latitudinal Variation
  672.         T capPhi = dateMJD.subtract(36204.0).divide(365.2422);
  673.         capPhi = capPhi.subtract(FastMath.floor(capPhi.getReal()));
  674.         final int signum = (satLat.getReal() >= 0.) ? 1 : -1;
  675.         final T sinLat = satLat.sin();
  676.         final T hm90  = altKm.subtract(90.);
  677.         final T dlrsl = hm90.multiply(0.02).multiply(hm90.multiply(-0.045).exp()).
  678.                         multiply(capPhi.multiply(MathUtils.TWO_PI).add(1.72).sin()).
  679.                         multiply(signum).multiply(sinLat).multiply(sinLat);

  681.         // Equation (23) - Computes the semiannual variation
  682.         T dlrsa = field.getZero();
  683.         if (z.getReal() < 2000.0) {
  684.             // Use new semiannual model dLog(rho)
  685.             dlrsa = semian08(dayOfYear(dateMJD), altKm, f10B, s10B, xm10B);
  686.         }

  688.         // Sum the delta-log-rhos and apply to the number densities.
  689.         // In CIRA72 the following equation contains an actual sum,
  690.         // namely DLR = LOG10 * (DLRGM + DLRSA + DLRSL)
  691.         // However, for Jacchia 70, there is no DLRGM or DLRSA.
  692.         final T dlr = dlrsl.add(dlrsa).multiply(LOG10);
  693.         for (int i = 0; i < 6; ++i) {
  694.             aln[i] = aln[i].add(dlr);
  695.         }

  697.         // Compute mass-density and mean-molecular-weight and
  698.         // convert number density logs from natural to common.
  699.         T sumnm = field.getZero();
  700.         for (int i = 0; i < 6; ++i) {
  701.             sumnm = sumnm.add(aln[i].exp().multiply(AMW[i]));
  702.         }
  703.         rho = sumnm.divide(AVOGAD);

  705.         // Compute the high altitude exospheric density correction factor
  706.         T fex = field.getOne();
  707.         if ((altKm.getReal() >= 1000.0) && (altKm.getReal() < 1500.0)) {
  708.             final T zeta = altKm.subtract(1000.).multiply(0.002);
  709.             final double f15c     = CHT[0] + CHT[1] * f10B + (CHT[2] + CHT[3] * f10B) * 1500.0;
  710.             final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
  711.             final double fex2     = 3.0 * f15c - f15cZeta - 3.0;
  712.             final double fex3     = f15cZeta - 2.0 * f15c + 2.0;
  713.             fex = fex.add(zeta.multiply(zeta).multiply(zeta.multiply(fex3).add(fex2)));
  714.         } else if (altKm.getReal() >= 1500.0) {
  715.             fex = altKm.multiply(CHT[3] * f10B).add(altKm.multiply(CHT[2])).add(CHT[0] + CHT[1] * f10B);
  716.         }

  718.         // Apply the exospheric density correction factor.
  719.         rho = rho.multiply(fex);

  721.         return rho;
  722.     }

  724.     /** Compute daily temperature correction for Jacchia-Bowman model.
  725.      * @param f10 solar flux index
  726.      * @param solarTime local solar time (hours in [0, 24[)
  727.      * @param satLat sat lat (radians)
  728.      * @param satAlt height (km)
  729.      * @return dTc correction
  730.      */
  731.     private static double dTc(final double f10, final double solarTime,
  732.                               final double satLat, final double satAlt) {
  733.         double dTc = 0.;
  734.         final double st = solarTime / 24.0;
  735.         final double cs = FastMath.cos(satLat);
  736.         final double fs = (f10 - 100.0) / 100.0;

  738.         // Calculates dTc according to height
  739.         if ((satAlt >= 120) && (satAlt <= 200)) {
  740.             final double dtc200 = poly2CDTC(fs, st, cs);
  741.             final double dtc200dz = poly1CDTC(fs, st, cs);
  742.             final double cc = 3.0 * dtc200 - dtc200dz;
  743.             final double dd = dtc200 - cc;
  744.             final double zp = (satAlt - 120.0) / 80.0;
  745.             dTc = zp * zp * (cc + dd * zp);
  746.         } else if ((satAlt > 200.0) && (satAlt <= 240.0)) {
  747.             final double h = (satAlt - 200.0) / 50.0;
  748.             dTc = poly1CDTC(fs, st, cs) * h + poly2CDTC(fs, st, cs);
  749.         } else if ((satAlt > 240.0) && (satAlt <= 300.0)) {
  750.             final double h = 0.8;
  751.             final double bb = poly1CDTC(fs, st, cs);
  752.             final double aa = bb * h + poly2CDTC(fs, st, cs);
  753.             final double p2BDT = poly2BDTC(st);
  754.             final double dtc300 = poly1BDTC(fs, st, cs, 3 * p2BDT);
  755.             final double dtc300dz = cs * p2BDT;
  756.             final double cc = 3.0 * dtc300 - dtc300dz - 3.0 * aa - 2.0 * bb;
  757.             final double dd = dtc300 - aa - bb - cc;
  758.             final double zp = (satAlt - 240.0) / 60.0;
  759.             dTc = aa + zp * (bb + zp * (cc + zp * dd));
  760.         } else if ((satAlt > 300.0) && (satAlt <= 600.0)) {
  761.             final double h = satAlt / 100.0;
  762.             dTc = poly1BDTC(fs, st, cs, h * poly2BDTC(st));
  763.         } else if ((satAlt > 600.0) && (satAlt <= 800.0)) {
  764.             final double poly2 = poly2BDTC(st);
  765.             final double aa = poly1BDTC(fs, st, cs, 6 * poly2);
  766.             final double bb = cs * poly2;
  767.             final double cc = -(3.0 * aa + 4.0 * bb) / 4.0;
  768.             final double dd = (aa + bb) / 4.0;
  769.             final double zp = (satAlt - 600.0) / 100.0;
  770.             dTc = aa + zp * (bb + zp * (cc + zp * dd));
  771.         }

  773.         return dTc;
  774.     }

  776.     /** Compute daily temperature correction for Jacchia-Bowman model.
  777.      * @param f10 solar flux index
  778.      * @param solarTime local solar time (hours in [0, 24[)
  779.      * @param satLat sat lat (radians)
  780.      * @param satAlt height (km)
  781.      * @param <T> type of the filed elements
  782.      * @return dTc correction
  783.      */
  784.     private static <T extends RealFieldElement<T>> T dTc(final double f10, final T solarTime,
  785.                                                          final T satLat, final T satAlt) {
  786.         T dTc = solarTime.getField().getZero();
  787.         final T      st = solarTime.divide(24.0);
  788.         final T      cs = satLat.cos();
  789.         final double fs = (f10 - 100.0) / 100.0;

  791.         // Calculates dTc according to height
  792.         if ((satAlt.getReal() >= 120) && (satAlt.getReal() <= 200)) {
  793.             final T dtc200   = poly2CDTC(fs, st, cs);
  794.             final T dtc200dz = poly1CDTC(fs, st, cs);
  795.             final T cc       = dtc200.multiply(3).subtract(dtc200dz);
  796.             final T dd       = dtc200.subtract(cc);
  797.             final T zp       = satAlt.subtract(120.0).divide(80.0);
  798.             dTc = zp.multiply(zp).multiply(cc.add(dd.multiply(zp)));
  799.         } else if ((satAlt.getReal() > 200.0) && (satAlt.getReal() <= 240.0)) {
  800.             final T h = satAlt.subtract(200.0).divide(50.0);
  801.             dTc = poly1CDTC(fs, st, cs).multiply(h).add(poly2CDTC(fs, st, cs));
  802.         } else if ((satAlt.getReal() > 240.0) && (satAlt.getReal() <= 300.0)) {
  803.             final T h = solarTime.getField().getZero().add(0.8);
  804.             final T bb = poly1CDTC(fs, st, cs);
  805.             final T aa = bb.multiply(h).add(poly2CDTC(fs, st, cs));
  806.             final T p2BDT = poly2BDTC(st);
  807.             final T dtc300 = poly1BDTC(fs, st, cs, p2BDT.multiply(3));
  808.             final T dtc300dz = cs.multiply(p2BDT);
  809.             final T cc = dtc300.multiply(3).subtract(dtc300dz).subtract(aa.multiply(3)).subtract(bb.multiply(2));
  810.             final T dd = dtc300.subtract(aa).subtract(bb).subtract(cc);
  811.             final T zp = satAlt.subtract(240.0).divide(60.0);
  812.             dTc = aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
  813.         } else if ((satAlt.getReal() > 300.0) && (satAlt.getReal() <= 600.0)) {
  814.             final T h = satAlt.divide(100.0);
  815.             dTc = poly1BDTC(fs, st, cs, h.multiply(poly2BDTC(st)));
  816.         } else if ((satAlt.getReal() > 600.0) && (satAlt.getReal() <= 800.0)) {
  817.             final T poly2 = poly2BDTC(st);
  818.             final T aa = poly1BDTC(fs, st, cs, poly2.multiply(6));
  819.             final T bb = cs.multiply(poly2);
  820.             final T cc = aa.multiply(3).add(bb.multiply(4)).divide(-4.0);
  821.             final T dd = aa.add(bb).divide(4.0);
  822.             final T zp = satAlt.subtract(600.0).divide(100.0);
  823.             dTc = aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
  824.         }

  826.         return dTc;
  827.     }

  829.     /** Calculates first polynomial with CDTC array.
  830.      * @param fs scaled flux f10
  831.      * @param st local solar time in [0, 1[
  832.      * @param cs cosine of satLat
  833.      * @return the value of the polynomial
  834.      */
  835.     private static double poly1CDTC(final double fs, final double st, final double cs) {
  836.         return CDTC[0] +
  837.                 fs * (CDTC[1] + st * (CDTC[2] + st * (CDTC[3] + st * (CDTC[4] + st * (CDTC[5] + st * CDTC[6]))))) +
  838.                 cs * st * (CDTC[7] + st * (CDTC[8] + st * (CDTC[9] + st * (CDTC[10] + st * CDTC[11])))) +
  839.                 cs * (CDTC[12] + fs * (CDTC[13] + st * (CDTC[14] + st * CDTC[15])));
  840.     }

  842.     /** Calculates first polynomial with CDTC array.
  843.      * @param fs scaled flux f10
  844.      * @param st local solar time in [0, 1[
  845.      * @param cs cosine of satLat
  846.      * @param <T> type fo the field elements
  847.      * @return the value of the polynomial
  848.      */
  849.     private static <T extends RealFieldElement<T>>  T poly1CDTC(final double fs, final T st, final T cs) {
  850.         return    st.multiply(CDTC[6]).
  851.               add(CDTC[5]).multiply(st).
  852.               add(CDTC[4]).multiply(st).
  853.               add(CDTC[3]).multiply(st).
  854.               add(CDTC[2]).multiply(st).
  855.               add(CDTC[1]).multiply(fs).
  856.               add(st.multiply(CDTC[11]).
  857.                   add(CDTC[10]).multiply(st).
  858.                   add(CDTC[ 9]).multiply(st).
  859.                   add(CDTC[ 8]).multiply(st).
  860.                   add(CDTC[7]).multiply(st).multiply(cs)).
  861.               add(st.multiply(CDTC[15]).
  862.                   add(CDTC[14]).multiply(st).
  863.                   add(CDTC[13]).multiply(fs).
  864.                   add(CDTC[12]).multiply(cs)).
  865.               add(CDTC[0]);
  866.     }

  868.     /** Calculates second polynomial with CDTC array.
  869.      * @param fs scaled flux f10
  870.      * @param st local solar time in [0, 1[
  871.      * @param cs cosine of satLat
  872.      * @return the value of the polynomial
  873.      */
  874.     private static double poly2CDTC(final double fs, final double st, final double cs) {
  875.         return CDTC[16] + st * cs * (CDTC[17] + st * (CDTC[18] + st * CDTC[19])) +
  876.                           fs * cs * (CDTC[20] + st * (CDTC[21] + st * CDTC[22]));
  877.     }

  879.     /** Calculates second polynomial with CDTC array.
  880.      * @param fs scaled flux f10
  881.      * @param st local solar time in [0, 1[
  882.      * @param cs cosine of satLat
  883.      * @param <T> type fo the field elements
  884.      * @return the value of the polynomial
  885.      */
  886.     private static <T extends RealFieldElement<T>>  T poly2CDTC(final double fs, final T st, final T cs) {
  887.         return         st.multiply(CDTC[19]).
  888.                    add(CDTC[18]).multiply(st).
  889.                    add(CDTC[17]).multiply(cs).multiply(st).
  890.                add(    st.multiply(CDTC[22]).
  891.                    add(CDTC[21]).multiply(st).
  892.                    add(CDTC[20]).multiply(cs).multiply(fs)).
  893.                add(CDTC[16]);
  894.     }

  896.     /** Calculates first polynomial with BDTC array.
  897.      * @param fs scaled flux f10
  898.      * @param st local solar time in [0, 1[
  899.      * @param cs cosine of satLat
  900.      * @param hp scaled height * poly2BDTC
  901.      * @return the value of the polynomial
  902.      */
  903.     private static double poly1BDTC(final double fs, final double st, final double cs, final double hp) {
  904.         return BDTC[0] +
  905.                 fs * (BDTC[1] + st * (BDTC[2] + st * (BDTC[3] + st * (BDTC[4] + st * (BDTC[5] + st * BDTC[6]))))) +
  906.                 cs * (st * (BDTC[7] + st * (BDTC[8] + st * (BDTC[9] + st * (BDTC[10] + st * BDTC[11])))) + hp + BDTC[18]);
  907.     }

  909.     /** Calculates first polynomial with BDTC array.
  910.      * @param fs scaled flux f10
  911.      * @param st local solar time in [0, 1[
  912.      * @param cs cosine of satLat
  913.      * @param hp scaled height * poly2BDTC
  914.      * @param <T> type fo the field elements
  915.      * @return the value of the polynomial
  916.      */
  917.     private static <T extends RealFieldElement<T>>  T poly1BDTC(final double fs, final T st, final T cs, final T hp) {
  918.         return     st.multiply(BDTC[6]).
  919.                add(BDTC[5]).multiply(st).
  920.                add(BDTC[4]).multiply(st).
  921.                add(BDTC[3]).multiply(st).
  922.                add(BDTC[2]).multiply(st).
  923.                add(BDTC[1]).multiply(fs).
  924.                add(    st.multiply(BDTC[11]).
  925.                    add(BDTC[10]).multiply(st).
  926.                    add(BDTC[ 9]).multiply(st).
  927.                    add(BDTC[ 8]).multiply(st).
  928.                    add(BDTC[ 7]).multiply(st).
  929.                    add(hp).add(BDTC[18]).multiply(cs)).
  930.                add(BDTC[0]);
  931.     }

  933.     /** Calculates second polynomial with BDTC array.
  934.      * @param st local solar time in [0, 1[
  935.      * @return the value of the polynomial
  936.      */
  937.     private static double poly2BDTC(final double st) {
  938.         return BDTC[12] + st * (BDTC[13] + st * (BDTC[14] + st * (BDTC[15] + st * (BDTC[16] + st * BDTC[17]))));
  939.     }

  941.     /** Calculates second polynomial with BDTC array.
  942.      * @param st local solar time in [0, 1[
  943.      * @param <T> type fo the field elements
  944.      * @return the value of the polynomial
  945.      */
  946.     private static <T extends RealFieldElement<T>>  T poly2BDTC(final T st) {
  947.         return     st.multiply(BDTC[17]).
  948.                add(BDTC[16]).multiply(st).
  949.                add(BDTC[15]).multiply(st).
  950.                add(BDTC[14]).multiply(st).
  951.                add(BDTC[13]).multiply(st).
  952.                add(BDTC[12]);
  953.     }

  955.     /** Evaluates mean molecualr mass - Equation (1).
  956.      * @param z altitude (km)
  957.      * @return mean molecular mass
  958.      */
  959.     private static double mBar(final double z) {
  960.         final double dz = z - 100.;
  961.         double amb = CMB[6];
  962.         for (int i = 5; i >= 0; --i) {
  963.             amb = dz * amb + CMB[i];
  964.         }
  965.         return amb;
  966.     }

  968.     /** Evaluates mean molecualr mass - Equation (1).
  969.      * @param z altitude (km)
  970.      * @return mean molecular mass
  971.      * @param <T> type fo the field elements
  972.      */
  973.     private static <T extends RealFieldElement<T>>  T mBar(final T z) {
  974.         final T dz = z.subtract(100.);
  975.         T amb = z.getField().getZero().add(CMB[6]);
  976.         for (int i = 5; i >= 0; --i) {
  977.             amb = dz.multiply(amb).add(CMB[i]);
  978.         }
  979.         return amb;
  980.     }

  982.     /** Evaluates the local temperature, Eq. (10) or (13) depending on altitude.
  983.      * @param z altitude
  984.      * @param tc tc array
  985.      * @return temperature profile
  986.      */
  987.     private static double localTemp(final double z, final double[] tc) {
  988.         final double dz = z - 125.;
  989.         if (dz <= 0.) {
  990.             return ((-9.8204695e-6 * dz - 7.3039742e-4) * dz * dz + 1.0) * dz * tc[1] + tc[0];
  991.         } else {
  992.             return tc[0] + tc[2] * FastMath.atan(tc[3] * dz * (1 + 4.5e-6 * FastMath.pow(dz, 2.5)));
  993.         }
  994.     }

  996.     /** Evaluates the local temperature, Eq. (10) or (13) depending on altitude.
  997.      * @param z altitude
  998.      * @param tc tc array
  999.      * @return temperature profile
  1000.      * @param <T> type fo the field elements
  1001.      */
  1002.     private static <T extends RealFieldElement<T>>  T localTemp(final T z, final T[] tc) {
  1003.         final T dz = z.subtract(125.);
  1004.         if (dz.getReal() <= 0.) {
  1005.             return dz.multiply(-9.8204695e-6).subtract(7.3039742e-4).multiply(dz).multiply(dz).add(1.0).multiply(dz).multiply(tc[1]).add(tc[0]);
  1006.         } else {
  1007.             return dz.multiply(dz).multiply(dz.sqrt()).multiply(4.5e-6).add(1).multiply(dz).multiply(tc[3]).atan().multiply(tc[2]).add(tc[0]);
  1008.         }
  1009.     }

  1011.     /** Evaluates the gravity at the altitude - Equation (8).
  1012.      * @param z altitude (km)
  1013.      * @return the gravity (m/s2)
  1014.      */
  1015.     private static double gravity(final double z) {
  1016.         final double tmp = 1.0 + z / EARTH_RADIUS;
  1017.         return Constants.G0_STANDARD_GRAVITY / (tmp * tmp);
  1018.     }

  1020.     /** Evaluates the gravity at the altitude - Equation (8).
  1021.      * @param z altitude (km)
  1022.      * @return the gravity (m/s2)
  1023.      * @param <T> type fo the field elements
  1024.      */
  1025.     private static <T extends RealFieldElement<T>>  T gravity(final T z) {
  1026.         final T tmp = z.divide(EARTH_RADIUS).add(1);
  1027.         return tmp.multiply(tmp).reciprocal().multiply(Constants.G0_STANDARD_GRAVITY);
  1028.     }

  1030.     /** Compute semi-annual variation (delta log(rho)).
  1031.      * @param doy day of year
  1032.      * @param alt height (km)
  1033.      * @param f10B average 81-day centered f10
  1034.      * @param s10B average 81-day centered s10
  1035.      * @param xm10B average 81-day centered xn10
  1036.      * @return semi-annual variation
  1037.      */
  1038.     private static double semian08(final double doy, final double alt,
  1039.                                    final double f10B, final double s10B, final double xm10B) {

  1041.         final double htz = alt / 1000.0;

  1043.         // COMPUTE NEW 81-DAY CENTERED SOLAR INDEX FOR FZ
  1044.         double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;

  1046.         // SEMIANNUAL AMPLITUDE
  1047.         final double fzz = FZM[0] + fsmb * (FZM[1] + htz * (FZM[2] + FZM[3] * htz + FZM[4] * fsmb));

  1049.         // COMPUTE DAILY 81-DAY CENTERED SOLAR INDEX FOR GT
  1050.         fsmb  = f10B - 0.75 * s10B - 0.37 * xm10B;

  1052.         // SEMIANNUAL PHASE FUNCTION
  1053.         final double tau   = MathUtils.TWO_PI * (doy - 1.0) / 365;
  1054.         final double taux2 = tau + tau;
  1055.         final double sin1P = FastMath.sin(tau);
  1056.         final double cos1P = FastMath.cos(tau);
  1057.         final double sin2P = FastMath.sin(taux2);
  1058.         final double cos2P = FastMath.cos(taux2);
  1059.         final double gtz = GTM[0] + GTM[1] * sin1P + GTM[2] * cos1P + GTM[3] * sin2P + GTM[4] * cos2P +
  1060.                            fsmb * (GTM[5] + GTM[6] * sin1P + GTM[7] * cos1P + GTM[8] * sin2P + GTM[9] * cos2P);

  1062.         return FastMath.max(1.0e-6, fzz) * gtz;

  1064.     }

  1066.     /** Compute semi-annual variation (delta log(rho)).
  1067.      * @param doy day of year
  1068.      * @param alt height (km)
  1069.      * @param f10B average 81-day centered f10
  1070.      * @param s10B average 81-day centered s10
  1071.      * @param xm10B average 81-day centered xn10
  1072.      * @return semi-annual variation
  1073.      * @param <T> type fo the field elements
  1074.      */
  1075.     private static <T extends RealFieldElement<T>>  T semian08(final T doy, final T alt,
  1076.                                                                final double f10B, final double s10B, final double xm10B) {

  1078.         final T htz = alt.divide(1000.0);

  1080.         // COMPUTE NEW 81-DAY CENTERED SOLAR INDEX FOR FZ
  1081.         double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;

  1083.         // SEMIANNUAL AMPLITUDE
  1084.         final T fzz = htz.multiply(FZM[3]).add(FZM[2] + FZM[4] * fsmb).multiply(htz).add(FZM[1]).multiply(fsmb).add(FZM[0]);

  1086.         // COMPUTE DAILY 81-DAY CENTERED SOLAR INDEX FOR GT
  1087.         fsmb  = f10B - 0.75 * s10B - 0.37 * xm10B;

  1089.         // SEMIANNUAL PHASE FUNCTION
  1090.         final T tau   = doy.subtract(1).divide(365).multiply(MathUtils.TWO_PI);
  1091.         final T taux2 = tau.add(tau);
  1092.         final T sin1P = tau.sin();
  1093.         final T cos1P = tau.cos();
  1094.         final T sin2P = taux2.sin();
  1095.         final T cos2P = taux2.cos();
  1096.         final T gtz =           cos2P.multiply(GTM[9]).
  1097.                             add(sin2P.multiply(GTM[8])).
  1098.                             add(cos1P.multiply(GTM[7])).
  1099.                             add(sin1P.multiply(GTM[6])).
  1100.                             add(GTM[5]).multiply(fsmb).
  1101.                         add(    cos2P.multiply(GTM[4]).
  1102.                             add(sin2P.multiply(GTM[3])).
  1103.                             add(cos1P.multiply(GTM[2])).
  1104.                             add(sin1P.multiply(GTM[1])).
  1105.                             add(GTM[0]));

  1107.         return fzz.getReal() > 1.0e-6 ? gtz.multiply(fzz) : gtz.multiply(1.0e-6);

  1109.     }

  1111.     /** Compute day of year.
  1112.      * @param dateMJD Modified Julian date
  1113.      * @return the number days in year
  1114.      */
  1115.     private static double dayOfYear(final double dateMJD) {
  1116.         final double d1950 = dateMJD - 33281;

  1118.         int iyday = (int) d1950;
  1119.         final double frac = d1950 - iyday;
  1120.         iyday = iyday + 364;

  1122.         int itemp = iyday / 1461;

  1124.         iyday = iyday - itemp * 1461;
  1125.         itemp = iyday / 365;
  1126.         if (itemp >= 3) {
  1127.             itemp = 3;
  1128.         }
  1129.         iyday = iyday - 365 * itemp + 1;
  1130.         return iyday + frac;
  1131.     }

  1133.     /** Compute day of year.
  1134.      * @param dateMJD Modified Julian date
  1135.      * @param <T> type of the field elements
  1136.      * @return the number days in year
  1137.      */
  1138.     private static <T extends RealFieldElement<T>> T dayOfYear(final T dateMJD) {
  1139.         final T d1950 = dateMJD.subtract(33281);

  1141.         int iyday = (int) d1950.getReal();
  1142.         final T frac = d1950.subtract(iyday);
  1143.         iyday = iyday + 364;

  1145.         int itemp = iyday / 1461;

  1147.         iyday = iyday - itemp * 1461;
  1148.         itemp = iyday / 365;
  1149.         if (itemp >= 3) {
  1150.             itemp = 3;
  1151.         }
  1152.         iyday = iyday - 365 * itemp + 1;
  1153.         return frac.add(iyday);
  1154.     }

  1156.     // OUTPUT:

  1158.     /** Compute min of two values, one double and one field element.
  1159.      * @param d double value
  1160.      * @param f field element
  1161.      * @param <T> type fo the field elements
  1162.      * @return min value
  1163.      */
  1164.     private <T extends RealFieldElement<T>> T min(final double d, final T f) {
  1165.         return (f.getReal() > d) ? f.getField().getZero().add(d) : f;
  1166.     }

  1168.     /** Compute max of two values, one double and one field element.
  1169.      * @param d double value
  1170.      * @param f field element
  1171.      * @param <T> type fo the field elements
  1172.      * @return max value
  1173.      */
  1174.     private <T extends RealFieldElement<T>> T max(final double d, final T f) {
  1175.         return (f.getReal() <= d) ? f.getField().getZero().add(d) : f;
  1176.     }

  1178.     /** Get the local density.
  1179.      * @param date current date
  1180.      * @param position current position in frame
  1181.      * @param frame the frame in which is defined the position
  1182.      * @return local density (kg/m³)
  1183.      * @exception OrekitException if date is out of range of solar activity
  1184.      */
  1185.     public double getDensity(final AbsoluteDate date, final Vector3D position,
  1186.                              final Frame frame)
  1187.         throws OrekitException {
  1188.         // check if data are available :
  1189.         if (date.compareTo(inputParams.getMaxDate()) > 0 ||
  1190.             date.compareTo(inputParams.getMinDate()) < 0) {
  1191.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  1192.                                       date, inputParams.getMinDate(), inputParams.getMaxDate());
  1193.         }

  1195.         // compute MJD date
  1196.         final double dateMJD = date.durationFrom(AbsoluteDate.MODIFIED_JULIAN_EPOCH) / Constants.JULIAN_DAY;

  1198.         // compute geodetic position
  1199.         final GeodeticPoint inBody = earth.transform(position, frame, date);

  1201.         // compute sun position
  1202.         final Frame ecef = earth.getBodyFrame();
  1203.         final Vector3D sunPos = sun.getPVCoordinates(date, ecef).getPosition();
  1204.         final GeodeticPoint sunInBody = earth.transform(sunPos, ecef, date);

  1206.         return getDensity(dateMJD,
  1207.                           sunInBody.getLongitude(), sunInBody.getLatitude(),
  1208.                           inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude(),
  1209.                           inputParams.getF10(date), inputParams.getF10B(date),
  1210.                           inputParams.getS10(date), inputParams.getS10B(date),
  1211.                           inputParams.getXM10(date), inputParams.getXM10B(date),
  1212.                           inputParams.getY10(date), inputParams.getY10B(date),
  1213.                           inputParams.getDSTDTC(date));

  1215.     }

  1217.     /** {@inheritDoc} */
  1218.     @Override
  1219.     public <T extends RealFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
  1220.                                                         final FieldVector3D<T> position,
  1221.                                                         final Frame frame)
  1222.         throws OrekitException {

  1224.         // check if data are available :
  1225.         final AbsoluteDate dateD = date.toAbsoluteDate();
  1226.         if ((dateD.compareTo(inputParams.getMaxDate()) > 0) ||
  1227.                         (dateD.compareTo(inputParams.getMinDate()) < 0)) {
  1228.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  1229.                                       dateD, inputParams.getMinDate(), inputParams.getMaxDate());
  1230.         }

  1232.         // compute MJD date
  1233.         final T dateMJD = date.durationFrom(AbsoluteDate.MODIFIED_JULIAN_EPOCH).divide(Constants.JULIAN_DAY);

  1235.         // compute geodetic position (km and °)
  1236.         final FieldGeodeticPoint<T> inBody = earth.transform(position, frame, date);

  1238.         // compute sun position
  1239.         final Frame ecef = earth.getBodyFrame();
  1240.         final FieldVector3D<T> sunPos = new FieldVector3D<>(date.getField(),
  1241.                         sun.getPVCoordinates(dateD, ecef).getPosition());
  1242.         final FieldGeodeticPoint<T> sunInBody = earth.transform(sunPos, ecef, date);

  1244.         return getDensity(dateMJD,
  1245.                           sunInBody.getLongitude(), sunInBody.getLatitude(),
  1246.                           inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude(),
  1247.                           inputParams.getF10(dateD), inputParams.getF10B(dateD),
  1248.                           inputParams.getS10(dateD), inputParams.getS10B(dateD),
  1249.                           inputParams.getXM10(dateD), inputParams.getXM10B(dateD),
  1250.                           inputParams.getY10(dateD), inputParams.getY10B(dateD),
  1251.                           inputParams.getDSTDTC(dateD));

  1253.     }

  1255. }
Created with JaCoCo 0.7.9.201702052155