JB2008.java

  1. /* Copyright 2002-2025 CS GROUP
  2.  * Licensed to CS GROUP (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.models.earth.atmosphere;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.CalculusFieldElement;
  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.FieldSinCos;
  25. import org.hipparchus.util.MathArrays;
  26. import org.hipparchus.util.MathUtils;
  27. import org.hipparchus.util.SinCos;
  28. import org.orekit.annotation.DefaultDataContext;
  29. import org.orekit.bodies.BodyShape;
  30. import org.orekit.bodies.FieldGeodeticPoint;
  31. import org.orekit.bodies.GeodeticPoint;
  32. import org.orekit.data.DataContext;
  33. import org.orekit.errors.OrekitException;
  34. import org.orekit.errors.OrekitMessages;
  35. import org.orekit.frames.Frame;
  36. import org.orekit.time.AbsoluteDate;
  37. import org.orekit.time.DateTimeComponents;
  38. import org.orekit.time.FieldAbsoluteDate;
  39. import org.orekit.time.TimeScale;
  40. import org.orekit.utils.Constants;
  41. import org.orekit.utils.ExtendedPositionProvider;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  168.     /** External data container. */
  169.     private JB2008InputParameters inputParams;

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

  174.     /** UTC time scale. */
  175.     private final TimeScale utc;

  177.     /** Constructor with space environment information for internal computation.
  178.      *
  179.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  180.      *
  181.      * @param parameters the solar and magnetic activity data
  182.      * @param sun the sun position
  183.      * @param earth the earth body shape
  184.      * @see #JB2008(JB2008InputParameters, ExtendedPositionProvider, BodyShape, TimeScale)
  185.      */
  186.     @DefaultDataContext
  187.     public JB2008(final JB2008InputParameters parameters,
  188.                   final ExtendedPositionProvider sun, final BodyShape earth) {
  189.         this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
  190.     }

  192.     /**
  193.      * Constructor with space environment information for internal computation.
  194.      *
  195.      * @param parameters the solar and magnetic activity data
  196.      * @param sun        the sun position
  197.      * @param earth      the earth body shape
  198.      * @param utc        UTC time scale. Used to computed the day fraction.
  199.      * @since 10.1
  200.      */
  201.     public JB2008(final JB2008InputParameters parameters,
  202.                   final ExtendedPositionProvider sun,
  203.                   final BodyShape earth,
  204.                   final TimeScale utc) {
  205.         super(sun);
  206.         this.earth = earth;
  207.         this.inputParams = parameters;
  208.         this.utc = utc;
  209.     }

  211.     /** {@inheritDoc} */
  212.     public Frame getFrame() {
  213.         return earth.getBodyFrame();
  214.     }

  216.     /** Get the local density with initial entries.
  217.      * @param dateMJD date and time, in modified julian days and fraction
  218.      * @param sunRA Right Ascension of Sun (radians)
  219.      * @param sunDecli Declination of Sun (radians)
  220.      * @param satLon Right Ascension of position (radians)
  221.      * @param satLat Geocentric latitude of position (radians)
  222.      * @param satAlt Height of position (m)
  223.      * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
  224.      *        (Tabular time 1.0 day earlier)
  225.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
  226.      *        (Tabular time 1.0 day earlier)
  227.      * @param s10 EUV index (26-34 nm) scaled to F10<br>
  228.      *        (Tabular time 1 day earlier)
  229.      * @param s10B UV 81-day averaged centered index
  230.      *        (Tabular time 1 day earlier)
  231.      * @param xm10 MG2 index scaled to F10<br>
  232.      *        (Tabular time 2.0 days earlier)
  233.      * @param xm10B MG2 81-day ave. centered index<br>
  234.      *        (Tabular time 2.0 days earlier)
  235.      * @param y10 Solar X-Ray &amp; Lya index scaled to F10<br>
  236.      *        (Tabular time 5.0 days earlier)
  237.      * @param y10B Solar X-Ray &amp; Lya 81-day ave. centered index<br>
  238.      *        (Tabular time 5.0 days earlier)
  239.      * @param dstdtc Temperature change computed from Dst index
  240.      * @return total mass-Density at input position (kg/m³)
  241.      */
  242.     public double getDensity(final double dateMJD, final double sunRA, final double sunDecli,
  243.                              final double satLon, final double satLat, final double satAlt,
  244.                              final double f10, final double f10B, final double s10,
  245.                              final double s10B, final double xm10, final double xm10B,
  246.                              final double y10, final double y10B, final double dstdtc) {

  248.         if (satAlt < ALT_MIN) {
  249.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, satAlt, ALT_MIN);
  250.         }

  252.         final double altKm = satAlt / 1000.0;

  254.         // Equation (14)
  255.         final double fn  = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
  256.         final double fsb = f10B * fn + s10B * (1. - fn);
  257.         final double tsubc = 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) +
  258.                              0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);

  260.         // Equation (15)
  261.         final double eta   = 0.5 * FastMath.abs(satLat - sunDecli);
  262.         final double theta = 0.5 * FastMath.abs(satLat + sunDecli);

  264.         // Equation (16)
  265.         final double h   = satLon - sunRA;
  266.         final double tau = h - 0.64577182 + 0.10471976 * FastMath.sin(h + 0.75049158);
  267.         double solarTime = FastMath.toDegrees(h + FastMath.PI) / 15.0;
  268.         while (solarTime >= 24) {
  269.             solarTime -= 24.;
  270.         }
  271.         while (solarTime < 0) {
  272.             solarTime += 24.;
  273.         }

  275.         // Equation (17)
  276.         final double cosEta  = FastMath.pow(FastMath.cos(eta), 2.5);
  277.         final double sinTeta = FastMath.pow(FastMath.sin(theta), 2.5);
  278.         final double cosTau  = FastMath.abs(FastMath.cos(0.5 * tau));
  279.         final double df = sinTeta + (cosEta - sinTeta) * cosTau * cosTau * cosTau;
  280.         final double tsubl = tsubc * (1. + 0.31 * df);

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

  285.         // Compute the local exospheric temperature.
  286.         // Add geomagnetic storm effect from input dTc value
  287.         final double[] temp = new double[2];
  288.         temp[0] = tsubl + dstdtc;
  289.         final double tinf = temp[0] + dtclst;

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

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

  297.         // The TC array will be an argument in the call to localTemp,
  298.         // which evaluates Eq. (10) or (13)
  299.         final double[] tc = new double[4];
  300.         tc[0] = tsubx;
  301.         tc[1] = gsubx;
  302.         // A of Equation (13)
  303.         tc[2] = (tinf - tsubx) * 2. / FastMath.PI;
  304.         tc[3] = gsubx / tc[2];

  306.         // Equation (5)
  307.         final double z1 = 90.;
  308.         final double z2 = FastMath.min(altKm, 105.0);
  309.         double al = FastMath.log(z2 / z1);
  310.         int n = 1 + (int) FastMath.floor(al / R1);
  311.         double zr = FastMath.exp(al / n);
  312.         final double mb1 = mBar(z1);
  313.         final double tloc1 = localTemp(z1, tc);
  314.         double zend  = z1;
  315.         double sub2  = 0.;
  316.         double ain   = mb1 * gravity(z1) / tloc1;
  317.         double mb2   = 0;
  318.         double tloc2 = 0;
  319.         double z     = 0;
  320.         double gravl = 0;

  322.         for (int i = 0; i < n; ++i) {
  323.             z = zend;
  324.             zend = zr * z;
  325.             final double dz = 0.25 * (zend - z);
  326.             double sum1 = WT[0] * ain;
  327.             for (int j = 1; j < 5; ++j) {
  328.                 z += dz;
  329.                 mb2   = mBar(z);
  330.                 tloc2 = localTemp(z, tc);
  331.                 gravl = gravity(z);
  332.                 ain   = mb2 * gravl / tloc2;
  333.                 sum1 += WT[j] * ain;
  334.             }
  335.             sub2 += dz * sum1;
  336.         }

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

  340.         // Equation (2)
  341.         final double anm = AVOGAD * rho;
  342.         final double an  = anm / mb2;

  344.         // Equation (3)
  345.         double fact2  = anm / 28.960;
  346.         final double[] aln = new double[6];
  347.         aln[0] = FastMath.log(FRAC[0] * fact2);
  348.         aln[3] = FastMath.log(FRAC[2] * fact2);
  349.         aln[4] = FastMath.log(FRAC[3] * fact2);
  350.         // Equation (4)
  351.         aln[1] = FastMath.log(fact2 * (1. + FRAC[1]) - an);
  352.         aln[2] = FastMath.log(2. * (an - fact2));

  354.         if (altKm <= 105.0) {
  355.             temp[1] = tloc2;
  356.             // Put in negligible hydrogen for use in DO-LOOP 13
  357.             aln[5] = aln[4] - 25.0;
  358.         } else {
  359.             // Equation (6)
  360.             al   = FastMath.log(FastMath.min(altKm, 500.0) / z);
  361.             n    = 1 + (int) FastMath.floor(al / R2);
  362.             zr   = FastMath.exp(al / n);
  363.             sub2 = 0.;
  364.             ain  = gravl / tloc2;

  366.             double tloc3 = 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.                     tloc3 = localTemp(z, tc);
  375.                     gravl = gravity(z);
  376.                     ain   = gravl / tloc3;
  377.                     sum1 += WT[j] * ain;
  378.                 }
  379.                 sub2 += dz * sum1;
  380.             }

  382.             al = FastMath.log(FastMath.max(altKm, 500.0) / z);
  383.             final double r = (altKm > 500.0) ? R3 : R2;
  384.             n = 1 + (int) FastMath.floor(al / r);
  385.             zr = FastMath.exp(al / n);
  386.             double sum3 = 0.;
  387.             double tloc4 = 0;
  388.             for (int i = 0; i < n; ++i) {
  389.                 z = zend;
  390.                 zend = zr * z;
  391.                 final double dz = 0.25 * (zend - z);
  392.                 double sum1 = WT[0] * ain;
  393.                 for (int j = 1; j < 5; ++j) {
  394.                     z += dz;
  395.                     tloc4 = localTemp(z, tc);
  396.                     gravl = gravity(z);
  397.                     ain   = gravl / tloc4;
  398.                     sum1 += WT[j] * ain;
  399.                 }
  400.                 sum3 += dz * sum1;
  401.             }
  402.             final double altr;
  403.             final double hSign;
  404.             if (altKm <= 500.) {
  405.                 temp[1] = tloc3;
  406.                 altr = FastMath.log(tloc3 / tloc2);
  407.                 fact2 = sub2 / RSTAR;
  408.                 hSign = 1.0;
  409.             } else {
  410.                 temp[1] = tloc4;
  411.                 altr = FastMath.log(tloc4 / tloc2);
  412.                 fact2 = (sub2 + sum3) / RSTAR;
  413.                 hSign = -1.0;
  414.             }
  415.             for (int i = 0; i < 5; ++i) {
  416.                 aln[i] = aln[i] - (1.0 + ALPHA[i]) * altr - fact2 * AMW[i];
  417.             }

  419.             // Equation (7)
  420.             final double al10t5 = FastMath.log10(tinf);
  421.             final double alnh5 = (5.5 * al10t5 - 39.40) * al10t5 + 73.13;
  422.             aln[5] = LOG10 * (alnh5 + 6.) + hSign * (FastMath.log(tloc4 / tloc3) + sum3 * AMW[5] / RSTAR);
  423.         }

  425.         // Equation (24) - J70 Seasonal-Latitudinal Variation
  426.         final double capPhi = ((dateMJD - 36204.0) / 365.2422) % 1;
  427.         final int signum = (satLat >= 0.) ? 1 : -1;
  428.         final double sinLat = FastMath.sin(satLat);
  429.         final double hm90  = altKm - 90.;
  430.         final double dlrsl = 0.02 * hm90 * FastMath.exp(-0.045 * hm90) *
  431.                              signum * FastMath.sin(MathUtils.TWO_PI * capPhi + 1.72) * sinLat * sinLat;

  433.         // Equation (23) - Computes the semiannual variation
  434.         double dlrsa = 0;
  435.         if (z < 2000.0) {
  436.             // Use new semiannual model dLog(rho)
  437.             dlrsa = semian08(dayOfYear(dateMJD), altKm, f10B, s10B, xm10B);
  438.         }

  440.         // Sum the delta-log-rhos and apply to the number densities.
  441.         // In CIRA72 the following equation contains an actual sum,
  442.         // namely DLR = LOG10 * (DLRGM + DLRSA + DLRSL)
  443.         // However, for Jacchia 70, there is no DLRGM or DLRSA.
  444.         final double dlr = LOG10 * (dlrsl + dlrsa);
  445.         for (int i = 0; i < 6; ++i) {
  446.             aln[i] += dlr;
  447.         }

  449.         // Compute mass-density and mean-molecular-weight and
  450.         // convert number density logs from natural to common.
  451.         double sumnm = 0.0;
  452.         for (int i = 0; i < 6; ++i) {
  453.             sumnm += FastMath.exp(aln[i]) * AMW[i];
  454.         }
  455.         rho = sumnm / AVOGAD;

  457.         // Compute the high altitude exospheric density correction factor
  458.         double fex = 1.0;
  459.         if (altKm >= 1000.0 && altKm < 1500.0) {
  460.             final double zeta = (altKm - 1000.) * 0.002;
  461.             final double f15c = CHT[0] + CHT[1] * f10B + (CHT[2] + CHT[3] * f10B) * 1500.0;
  462.             final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
  463.             final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
  464.             final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
  465.             fex += zeta * zeta * (fex2 + fex3 * zeta);
  466.         } else if (altKm >= 1500.0) {
  467.             fex = CHT[0] + CHT[1] * f10B + CHT[2] * altKm + CHT[3] * f10B * altKm;
  468.         }

  470.         // Apply the exospheric density correction factor.
  471.         rho *= fex;

  473.         return rho;
  474.     }

  476.     /** Get the local density with initial entries.
  477.      * @param dateMJD date and time, in modified julian days and fraction
  478.      * @param sunRA Right Ascension of Sun (radians)
  479.      * @param sunDecli Declination of Sun (radians)
  480.      * @param satLon Right Ascension of position (radians)
  481.      * @param satLat Geocentric latitude of position (radians)
  482.      * @param satAlt Height of position (m)
  483.      * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
  484.      *        (Tabular time 1.0 day earlier)
  485.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
  486.      *        (Tabular time 1.0 day earlier)
  487.      * @param s10 EUV index (26-34 nm) scaled to F10<br>
  488.      *        (Tabular time 1 day earlier)
  489.      * @param s10B UV 81-day averaged centered index
  490.      *        (Tabular time 1 day earlier)
  491.      * @param xm10 MG2 index scaled to F10<br>
  492.      *        (Tabular time 2.0 days earlier)
  493.      * @param xm10B MG2 81-day ave. centered index<br>
  494.      *        (Tabular time 2.0 days earlier)
  495.      * @param y10 Solar X-Ray &amp; Lya index scaled to F10<br>
  496.      *        (Tabular time 5.0 days earlier)
  497.      * @param y10B Solar X-Ray &amp; Lya 81-day ave. centered index<br>
  498.      *        (Tabular time 5.0 days earlier)
  499.      * @param dstdtc Temperature change computed from Dst index
  500.      * @param <T> type of the field elements
  501.      * @return total mass-Density at input position (kg/m³)
  502.      */
  503.     public <T extends CalculusFieldElement<T>> T getDensity(final T dateMJD, final T sunRA, final T sunDecli,
  504.                                                         final T satLon, final T satLat, final T satAlt,
  505.                                                         final double f10, final double f10B, final double s10,
  506.                                                         final double s10B, final double xm10, final double xm10B,
  507.                                                         final double y10, final double y10B, final double dstdtc) {

  509.         if (satAlt.getReal() < ALT_MIN) {
  510.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
  511.                                       satAlt.getReal(), ALT_MIN);
  512.         }

  514.         final Field<T> field  = satAlt.getField();
  515.         final T pi    = field.getOne().getPi();
  516.         final T altKm = satAlt.divide(1000.0);

  518.         // Equation (14)
  519.         final double fn  = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
  520.         final double fsb = f10B * fn + s10B * (1. - fn);
  521.         final double tsubc = 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) +
  522.                              0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);

  524.         // Equation (15)
  525.         final T eta   = satLat.subtract(sunDecli).abs().multiply(0.5);
  526.         final T theta = satLat.add(sunDecli).abs().multiply(0.5);

  528.         // Equation (16)
  529.         final T h   = satLon.subtract(sunRA);
  530.         final T tau = h.subtract(0.64577182).add(h.add(0.75049158).sin().multiply(0.10471976));
  531.         T solarTime = FastMath.toDegrees(h.add(pi)).divide(15.0);
  532.         while (solarTime.getReal() >= 24) {
  533.             solarTime = solarTime.subtract(24);
  534.         }
  535.         while (solarTime.getReal() < 0) {
  536.             solarTime = solarTime.add(24);
  537.         }

  539.         // Equation (17)
  540.         final T cos     = eta.cos();
  541.         final T cosEta  = cos.square().multiply(cos.sqrt());
  542.         final T sin     = theta.sin();
  543.         final T sinTeta = sin.square().multiply(sin.sqrt());
  544.         final T cosTau  = tau.multiply(0.5).cos().abs();
  545.         final T df      = sinTeta.add(cosEta.subtract(sinTeta).multiply(cosTau).multiply(cosTau).multiply(cosTau));
  546.         final T tsubl   = df.multiply(0.31).add(1).multiply(tsubc);

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

  551.         // Compute the local exospheric temperature.
  552.         // Add geomagnetic storm effect from input dTc value
  553.         final T[] temp = MathArrays.buildArray(field, 2);
  554.         temp[0] = tsubl.add(dstdtc);
  555.         final T tinf = temp[0].add(dtclst);

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

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

  563.         // The TC array will be an argument in the call to localTemp,
  564.         // which evaluates Eq. (10) or (13)
  565.         final T[] tc = MathArrays.buildArray(field, 4);
  566.         tc[0] = tsubx;
  567.         tc[1] = gsubx;
  568.         // A of Equation (13)
  569.         tc[2] = tinf.subtract(tsubx).multiply(pi.reciprocal().multiply(2.0));
  570.         tc[3] = gsubx.divide(tc[2]);

  572.         // Equation (5)
  573.         final T z1 = field.getZero().newInstance(90.);
  574.         final T z2 = min(105.0, altKm);
  575.         T al = z2.divide(z1).log();
  576.         int n = 1 + (int) FastMath.floor(al.getReal() / R1);
  577.         T zr = al.divide(n).exp();
  578.         final T mb1 = mBar(z1);
  579.         final T tloc1 = localTemp(z1, tc);
  580.         T zend  = z1;
  581.         T sub2  = field.getZero();
  582.         T ain   = mb1.multiply(gravity(z1)).divide(tloc1);
  583.         T mb2   = field.getZero();
  584.         T tloc2 = field.getZero();
  585.         T z     = field.getZero();
  586.         T gravl = field.getZero();
  587.         for (int i = 0; i < n; ++i) {
  588.             z = zend;
  589.             zend = zr.multiply(z);
  590.             final T dz = zend.subtract(z).multiply(0.25);
  591.             T sum1 = ain.multiply(WT[0]);
  592.             for (int j = 1; j < 5; ++j) {
  593.                 z = z.add(dz);
  594.                 mb2   = mBar(z);
  595.                 tloc2 = localTemp(z, tc);
  596.                 gravl = gravity(z);
  597.                 ain   = mb2.multiply(gravl).divide(tloc2);
  598.                 sum1  = sum1.add(ain.multiply(WT[j]));
  599.             }
  600.             sub2 = sub2.add(dz.multiply(sum1));
  601.         }

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

  605.         // Equation (2)
  606.         final T anm = rho.multiply(AVOGAD);
  607.         final T an  = anm.divide(mb2);

  609.         // Equation (3)
  610.         T fact2  = anm.divide(28.960);
  611.         final T[] aln = MathArrays.buildArray(field, 6);
  612.         aln[0] = fact2.multiply(FRAC[0]).log();
  613.         aln[3] = fact2.multiply(FRAC[2]).log();
  614.         aln[4] = fact2.multiply(FRAC[3]).log();
  615.         // Equation (4)
  616.         aln[1] = fact2.multiply(1. + FRAC[1]).subtract(an).log();
  617.         aln[2] = an.subtract(fact2).multiply(2).log();

  619.         if (altKm.getReal() <= 105.0) {
  620.             temp[1] = tloc2;
  621.             // Put in negligible hydrogen for use in DO-LOOP 13
  622.             aln[5] = aln[4].subtract(25.0);
  623.         } else {
  624.             // Equation (6)
  625.             al   = min(500.0, altKm).divide(z).log();
  626.             n    = 1 + (int) FastMath.floor(al.getReal() / R2);
  627.             zr   = al.divide(n).exp();
  628.             sub2 = field.getZero();
  629.             ain  = gravl.divide(tloc2);

  631.             T tloc3 = field.getZero();
  632.             for (int i = 0; i < n; ++i) {
  633.                 z = zend;
  634.                 zend = zr.multiply(z);
  635.                 final T dz = zend.subtract(z).multiply(0.25);
  636.                 T sum1 = ain.multiply(WT[0]);
  637.                 for (int j = 1; j < 5; ++j) {
  638.                     z = z.add(dz);
  639.                     tloc3 = localTemp(z, tc);
  640.                     gravl = gravity(z);
  641.                     ain   = gravl.divide(tloc3);
  642.                     sum1  = sum1.add(ain.multiply(WT[j]));
  643.                 }
  644.                 sub2 = sub2.add(dz.multiply(sum1));
  645.             }

  647.             al = max(500.0, altKm).divide(z).log();
  648.             final double r = (altKm.getReal() > 500.0) ? R3 : R2;
  649.             n = 1 + (int) FastMath.floor(al.getReal() / r);
  650.             zr = al.divide(n).exp();
  651.             T sum3 = field.getZero();
  652.             T tloc4 = field.getZero();
  653.             for (int i = 0; i < n; ++i) {
  654.                 z = zend;
  655.                 zend = zr.multiply(z);
  656.                 final T dz = zend.subtract(z).multiply(0.25);
  657.                 T sum1 = ain.multiply(WT[0]);
  658.                 for (int j = 1; j < 5; ++j) {
  659.                     z = z.add(dz);
  660.                     tloc4 = localTemp(z, tc);
  661.                     gravl = gravity(z);
  662.                     ain   = gravl.divide(tloc4);
  663.                     sum1  = sum1.add(ain.multiply(WT[j]));
  664.                 }
  665.                 sum3 = sum3.add(dz.multiply(sum1));
  666.             }
  667.             final T altr;
  668.             final double hSign;
  669.             if (altKm.getReal() <= 500.) {
  670.                 temp[1] = tloc3;
  671.                 altr = tloc3.divide(tloc2).log();
  672.                 fact2 = sub2.divide(RSTAR);
  673.                 hSign = 1.0;
  674.             } else {
  675.                 temp[1] = tloc4;
  676.                 altr = tloc4.divide(tloc2).log();
  677.                 fact2 = sub2.add(sum3).divide(RSTAR);
  678.                 hSign = -1.0;
  679.             }
  680.             for (int i = 0; i < 5; ++i) {
  681.                 aln[i] = aln[i].subtract(altr.multiply(1.0 + ALPHA[i])).subtract(fact2.multiply(AMW[i]));
  682.             }

  684.             // Equation (7)
  685.             final T al10t5 = tinf.log10();
  686.             final T alnh5 = al10t5.multiply(5.5).subtract(39.40).multiply(al10t5).add(73.13);
  687.             aln[5] = alnh5.add(6.).multiply(LOG10).
  688.                      add(tloc4.divide(tloc3).log().add(sum3.multiply(AMW[5] / RSTAR)).multiply(hSign));
  689.         }

  691.         // Equation (24) - J70 Seasonal-Latitudinal Variation
  692.         T capPhi = dateMJD.subtract(36204.0).divide(365.2422);
  693.         capPhi = capPhi.subtract(FastMath.floor(capPhi.getReal()));
  694.         final int signum = (satLat.getReal() >= 0.) ? 1 : -1;
  695.         final T sinLat = satLat.sin();
  696.         final T hm90  = altKm.subtract(90.);
  697.         final T dlrsl = hm90.multiply(0.02).multiply(hm90.multiply(-0.045).exp()).
  698.                         multiply(capPhi.multiply(MathUtils.TWO_PI).add(1.72).sin()).
  699.                         multiply(signum).multiply(sinLat).multiply(sinLat);

  701.         // Equation (23) - Computes the semiannual variation
  702.         T dlrsa = field.getZero();
  703.         if (z.getReal() < 2000.0) {
  704.             // Use new semiannual model dLog(rho)
  705.             dlrsa = semian08(dayOfYear(dateMJD), altKm, f10B, s10B, xm10B);
  706.         }

  708.         // Sum the delta-log-rhos and apply to the number densities.
  709.         // In CIRA72 the following equation contains an actual sum,
  710.         // namely DLR = LOG10 * (DLRGM + DLRSA + DLRSL)
  711.         // However, for Jacchia 70, there is no DLRGM or DLRSA.
  712.         final T dlr = dlrsl.add(dlrsa).multiply(LOG10);
  713.         for (int i = 0; i < 6; ++i) {
  714.             aln[i] = aln[i].add(dlr);
  715.         }

  717.         // Compute mass-density and mean-molecular-weight and
  718.         // convert number density logs from natural to common.
  719.         T sumnm = field.getZero();
  720.         for (int i = 0; i < 6; ++i) {
  721.             sumnm = sumnm.add(aln[i].exp().multiply(AMW[i]));
  722.         }
  723.         rho = sumnm.divide(AVOGAD);

  725.         // Compute the high altitude exospheric density correction factor
  726.         T fex = field.getOne();
  727.         if (altKm.getReal() >= 1000.0 && altKm.getReal() < 1500.0) {
  728.             final T zeta = altKm.subtract(1000.).multiply(0.002);
  729.             final double f15c     = CHT[0] + CHT[1] * f10B + (CHT[2] + CHT[3] * f10B) * 1500.0;
  730.             final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
  731.             final double fex2     = 3.0 * f15c - f15cZeta - 3.0;
  732.             final double fex3     = f15cZeta - 2.0 * f15c + 2.0;
  733.             fex = fex.add(zeta.multiply(zeta).multiply(zeta.multiply(fex3).add(fex2)));
  734.         } else if (altKm.getReal() >= 1500.0) {
  735.             fex = altKm.multiply(CHT[3] * f10B).add(altKm.multiply(CHT[2])).add(CHT[0] + CHT[1] * f10B);
  736.         }

  738.         // Apply the exospheric density correction factor.
  739.         rho = rho.multiply(fex);

  741.         return rho;
  742.     }

  744.     /** Compute daily temperature correction for Jacchia-Bowman model.
  745.      * @param f10 solar flux index
  746.      * @param solarTime local solar time (hours in [0, 24[)
  747.      * @param satLat sat lat (radians)
  748.      * @param satAlt height (km)
  749.      * @return dTc correction
  750.      */
  751.     private static double dTc(final double f10, final double solarTime,
  752.                               final double satLat, final double satAlt) {
  753.         double dTc = 0.;
  754.         final double st = solarTime / 24.0;
  755.         final double cs = FastMath.cos(satLat);
  756.         final double fs = (f10 - 100.0) / 100.0;

  758.         // Calculates dTc according to height
  759.         if (satAlt >= 120 && satAlt <= 200) {
  760.             final double dtc200 = poly2CDTC(fs, st, cs);
  761.             final double dtc200dz = poly1CDTC(fs, st, cs);
  762.             final double cc = 3.0 * dtc200 - dtc200dz;
  763.             final double dd = dtc200 - cc;
  764.             final double zp = (satAlt - 120.0) / 80.0;
  765.             dTc = zp * zp * (cc + dd * zp);
  766.         } else if (satAlt > 200.0 && satAlt <= 240.0) {
  767.             final double h = (satAlt - 200.0) / 50.0;
  768.             dTc = poly1CDTC(fs, st, cs) * h + poly2CDTC(fs, st, cs);
  769.         } else if (satAlt > 240.0 && satAlt <= 300.0) {
  770.             final double h = 0.8;
  771.             final double bb = poly1CDTC(fs, st, cs);
  772.             final double aa = bb * h + poly2CDTC(fs, st, cs);
  773.             final double p2BDT = poly2BDTC(st);
  774.             final double dtc300 = poly1BDTC(fs, st, cs, 3 * p2BDT);
  775.             final double dtc300dz = cs * p2BDT;
  776.             final double cc = 3.0 * dtc300 - dtc300dz - 3.0 * aa - 2.0 * bb;
  777.             final double dd = dtc300 - aa - bb - cc;
  778.             final double zp = (satAlt - 240.0) / 60.0;
  779.             dTc = aa + zp * (bb + zp * (cc + zp * dd));
  780.         } else if (satAlt > 300.0 && satAlt <= 600.0) {
  781.             final double h = satAlt / 100.0;
  782.             dTc = poly1BDTC(fs, st, cs, h * poly2BDTC(st));
  783.         } else if (satAlt > 600.0 && satAlt <= 800.0) {
  784.             final double poly2 = poly2BDTC(st);
  785.             final double aa = poly1BDTC(fs, st, cs, 6 * poly2);
  786.             final double bb = cs * poly2;
  787.             final double cc = -(3.0 * aa + 4.0 * bb) / 4.0;
  788.             final double dd = (aa + bb) / 4.0;
  789.             final double zp = (satAlt - 600.0) / 100.0;
  790.             dTc = aa + zp * (bb + zp * (cc + zp * dd));
  791.         }

  793.         return dTc;
  794.     }

  796.     /** Compute daily temperature correction for Jacchia-Bowman model.
  797.      * @param f10 solar flux index
  798.      * @param solarTime local solar time (hours in [0, 24[)
  799.      * @param satLat sat lat (radians)
  800.      * @param satAlt height (km)
  801.      * @param <T> type of the filed elements
  802.      * @return dTc correction
  803.      */
  804.     private static <T extends CalculusFieldElement<T>> T dTc(final double f10, final T solarTime,
  805.                                                          final T satLat, final T satAlt) {
  806.         T dTc = solarTime.getField().getZero();
  807.         final T      st = solarTime.divide(24.0);
  808.         final T      cs = satLat.cos();
  809.         final double fs = (f10 - 100.0) / 100.0;

  811.         // Calculates dTc according to height
  812.         if (satAlt.getReal() >= 120 && satAlt.getReal() <= 200) {
  813.             final T dtc200   = poly2CDTC(fs, st, cs);
  814.             final T dtc200dz = poly1CDTC(fs, st, cs);
  815.             final T cc       = dtc200.multiply(3).subtract(dtc200dz);
  816.             final T dd       = dtc200.subtract(cc);
  817.             final T zp       = satAlt.subtract(120.0).divide(80.0);
  818.             dTc = zp.multiply(zp).multiply(cc.add(dd.multiply(zp)));
  819.         } else if (satAlt.getReal() > 200.0 && satAlt.getReal() <= 240.0) {
  820.             final T h = satAlt.subtract(200.0).divide(50.0);
  821.             dTc = poly1CDTC(fs, st, cs).multiply(h).add(poly2CDTC(fs, st, cs));
  822.         } else if (satAlt.getReal() > 240.0 && satAlt.getReal() <= 300.0) {
  823.             final T h = solarTime.getField().getZero().newInstance(0.8);
  824.             final T bb = poly1CDTC(fs, st, cs);
  825.             final T aa = bb.multiply(h).add(poly2CDTC(fs, st, cs));
  826.             final T p2BDT = poly2BDTC(st);
  827.             final T dtc300 = poly1BDTC(fs, st, cs, p2BDT.multiply(3));
  828.             final T dtc300dz = cs.multiply(p2BDT);
  829.             final T cc = dtc300.multiply(3).subtract(dtc300dz).subtract(aa.multiply(3)).subtract(bb.multiply(2));
  830.             final T dd = dtc300.subtract(aa).subtract(bb).subtract(cc);
  831.             final T zp = satAlt.subtract(240.0).divide(60.0);
  832.             dTc = aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
  833.         } else if (satAlt.getReal() > 300.0 && satAlt.getReal() <= 600.0) {
  834.             final T h = satAlt.divide(100.0);
  835.             dTc = poly1BDTC(fs, st, cs, h.multiply(poly2BDTC(st)));
  836.         } else if (satAlt.getReal() > 600.0 && satAlt.getReal() <= 800.0) {
  837.             final T poly2 = poly2BDTC(st);
  838.             final T aa = poly1BDTC(fs, st, cs, poly2.multiply(6));
  839.             final T bb = cs.multiply(poly2);
  840.             final T cc = aa.multiply(3).add(bb.multiply(4)).divide(-4.0);
  841.             final T dd = aa.add(bb).divide(4.0);
  842.             final T zp = satAlt.subtract(600.0).divide(100.0);
  843.             dTc = aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
  844.         }

  846.         return dTc;
  847.     }

  849.     /** Calculates first polynomial with CDTC array.
  850.      * @param fs scaled flux f10
  851.      * @param st local solar time in [0, 1[
  852.      * @param cs cosine of satLat
  853.      * @return the value of the polynomial
  854.      */
  855.     private static double poly1CDTC(final double fs, final double st, final double cs) {
  856.         return CDTC[0] +
  857.                 fs * (CDTC[1] + st * (CDTC[2] + st * (CDTC[3] + st * (CDTC[4] + st * (CDTC[5] + st * CDTC[6]))))) +
  858.                 cs * st * (CDTC[7] + st * (CDTC[8] + st * (CDTC[9] + st * (CDTC[10] + st * CDTC[11])))) +
  859.                 cs * (CDTC[12] + fs * (CDTC[13] + st * (CDTC[14] + st * CDTC[15])));
  860.     }

  862.     /** Calculates first polynomial with CDTC array.
  863.      * @param fs scaled flux f10
  864.      * @param st local solar time in [0, 1[
  865.      * @param cs cosine of satLat
  866.      * @param <T> type of the field elements
  867.      * @return the value of the polynomial
  868.      */
  869.     private static <T extends CalculusFieldElement<T>>  T poly1CDTC(final double fs, final T st, final T cs) {
  870.         return    st.multiply(CDTC[6]).
  871.               add(CDTC[5]).multiply(st).
  872.               add(CDTC[4]).multiply(st).
  873.               add(CDTC[3]).multiply(st).
  874.               add(CDTC[2]).multiply(st).
  875.               add(CDTC[1]).multiply(fs).
  876.               add(st.multiply(CDTC[11]).
  877.                   add(CDTC[10]).multiply(st).
  878.                   add(CDTC[ 9]).multiply(st).
  879.                   add(CDTC[ 8]).multiply(st).
  880.                   add(CDTC[7]).multiply(st).multiply(cs)).
  881.               add(st.multiply(CDTC[15]).
  882.                   add(CDTC[14]).multiply(st).
  883.                   add(CDTC[13]).multiply(fs).
  884.                   add(CDTC[12]).multiply(cs)).
  885.               add(CDTC[0]);
  886.     }

  888.     /** Calculates second polynomial with CDTC array.
  889.      * @param fs scaled flux f10
  890.      * @param st local solar time in [0, 1[
  891.      * @param cs cosine of satLat
  892.      * @return the value of the polynomial
  893.      */
  894.     private static double poly2CDTC(final double fs, final double st, final double cs) {
  895.         return CDTC[16] + st * cs * (CDTC[17] + st * (CDTC[18] + st * CDTC[19])) +
  896.                           fs * cs * (CDTC[20] + st * (CDTC[21] + st * CDTC[22]));
  897.     }

  899.     /** Calculates second polynomial with CDTC array.
  900.      * @param fs scaled flux f10
  901.      * @param st local solar time in [0, 1[
  902.      * @param cs cosine of satLat
  903.      * @param <T> type of the field elements
  904.      * @return the value of the polynomial
  905.      */
  906.     private static <T extends CalculusFieldElement<T>>  T poly2CDTC(final double fs, final T st, final T cs) {
  907.         return         st.multiply(CDTC[19]).
  908.                    add(CDTC[18]).multiply(st).
  909.                    add(CDTC[17]).multiply(cs).multiply(st).
  910.                add(    st.multiply(CDTC[22]).
  911.                    add(CDTC[21]).multiply(st).
  912.                    add(CDTC[20]).multiply(cs).multiply(fs)).
  913.                add(CDTC[16]);
  914.     }

  916.     /** Calculates first polynomial with BDTC array.
  917.      * @param fs scaled flux f10
  918.      * @param st local solar time in [0, 1[
  919.      * @param cs cosine of satLat
  920.      * @param hp scaled height * poly2BDTC
  921.      * @return the value of the polynomial
  922.      */
  923.     private static double poly1BDTC(final double fs, final double st, final double cs, final double hp) {
  924.         return BDTC[0] +
  925.                 fs * (BDTC[1] + st * (BDTC[2] + st * (BDTC[3] + st * (BDTC[4] + st * (BDTC[5] + st * BDTC[6]))))) +
  926.                 cs * (st * (BDTC[7] + st * (BDTC[8] + st * (BDTC[9] + st * (BDTC[10] + st * BDTC[11])))) + hp + BDTC[18]);
  927.     }

  929.     /** Calculates first polynomial with BDTC array.
  930.      * @param fs scaled flux f10
  931.      * @param st local solar time in [0, 1[
  932.      * @param cs cosine of satLat
  933.      * @param hp scaled height * poly2BDTC
  934.      * @param <T> type of the field elements
  935.      * @return the value of the polynomial
  936.      */
  937.     private static <T extends CalculusFieldElement<T>>  T poly1BDTC(final double fs, final T st, final T cs, final T hp) {
  938.         return     st.multiply(BDTC[6]).
  939.                add(BDTC[5]).multiply(st).
  940.                add(BDTC[4]).multiply(st).
  941.                add(BDTC[3]).multiply(st).
  942.                add(BDTC[2]).multiply(st).
  943.                add(BDTC[1]).multiply(fs).
  944.                add(    st.multiply(BDTC[11]).
  945.                    add(BDTC[10]).multiply(st).
  946.                    add(BDTC[ 9]).multiply(st).
  947.                    add(BDTC[ 8]).multiply(st).
  948.                    add(BDTC[ 7]).multiply(st).
  949.                    add(hp).add(BDTC[18]).multiply(cs)).
  950.                add(BDTC[0]);
  951.     }

  953.     /** Calculates second polynomial with BDTC array.
  954.      * @param st local solar time in [0, 1[
  955.      * @return the value of the polynomial
  956.      */
  957.     private static double poly2BDTC(final double st) {
  958.         return BDTC[12] + st * (BDTC[13] + st * (BDTC[14] + st * (BDTC[15] + st * (BDTC[16] + st * BDTC[17]))));
  959.     }

  961.     /** Calculates second polynomial with BDTC array.
  962.      * @param st local solar time in [0, 1[
  963.      * @param <T> type of the field elements
  964.      * @return the value of the polynomial
  965.      */
  966.     private static <T extends CalculusFieldElement<T>>  T poly2BDTC(final T st) {
  967.         return     st.multiply(BDTC[17]).
  968.                add(BDTC[16]).multiply(st).
  969.                add(BDTC[15]).multiply(st).
  970.                add(BDTC[14]).multiply(st).
  971.                add(BDTC[13]).multiply(st).
  972.                add(BDTC[12]);
  973.     }

  975.     /** Evaluates mean molecualr mass - Equation (1).
  976.      * @param z altitude (km)
  977.      * @return mean molecular mass
  978.      */
  979.     private static double mBar(final double z) {
  980.         final double dz = z - 100.;
  981.         double amb = CMB[6];
  982.         for (int i = 5; i >= 0; --i) {
  983.             amb = dz * amb + CMB[i];
  984.         }
  985.         return amb;
  986.     }

  988.     /** Evaluates mean molecualr mass - Equation (1).
  989.      * @param z altitude (km)
  990.      * @return mean molecular mass
  991.      * @param <T> type of the field elements
  992.      */
  993.     private static <T extends CalculusFieldElement<T>>  T mBar(final T z) {
  994.         final T dz = z.subtract(100.);
  995.         T amb = z.getField().getZero().newInstance(CMB[6]);
  996.         for (int i = 5; i >= 0; --i) {
  997.             amb = dz.multiply(amb).add(CMB[i]);
  998.         }
  999.         return amb;
  1000.     }

  1002.     /** Evaluates the local temperature, Eq. (10) or (13) depending on altitude.
  1003.      * @param z altitude
  1004.      * @param tc tc array
  1005.      * @return temperature profile
  1006.      */
  1007.     private static double localTemp(final double z, final double[] tc) {
  1008.         final double dz = z - 125.;
  1009.         if (dz <= 0.) {
  1010.             return ((-9.8204695e-6 * dz - 7.3039742e-4) * dz * dz + 1.0) * dz * tc[1] + tc[0];
  1011.         } else {
  1012.             return tc[0] + tc[2] * FastMath.atan(tc[3] * dz * (1 + 4.5e-6 * FastMath.pow(dz, 2.5)));
  1013.         }
  1014.     }

  1016.     /** Evaluates the local temperature, Eq. (10) or (13) depending on altitude.
  1017.      * @param z altitude
  1018.      * @param tc tc array
  1019.      * @return temperature profile
  1020.      * @param <T> type of the field elements
  1021.      */
  1022.     private static <T extends CalculusFieldElement<T>>  T localTemp(final T z, final T[] tc) {
  1023.         final T dz = z.subtract(125.);
  1024.         if (dz.getReal() <= 0.) {
  1025.             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]);
  1026.         } else {
  1027.             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]);
  1028.         }
  1029.     }

  1031.     /** Evaluates the gravity at the altitude - Equation (8).
  1032.      * @param z altitude (km)
  1033.      * @return the gravity (m/s2)
  1034.      */
  1035.     private static double gravity(final double z) {
  1036.         final double tmp = 1.0 + z / EARTH_RADIUS;
  1037.         return Constants.G0_STANDARD_GRAVITY / (tmp * tmp);
  1038.     }

  1040.     /** Evaluates the gravity at the altitude - Equation (8).
  1041.      * @param z altitude (km)
  1042.      * @return the gravity (m/s2)
  1043.      * @param <T> type of the field elements
  1044.      */
  1045.     private static <T extends CalculusFieldElement<T>>  T gravity(final T z) {
  1046.         final T tmp = z.divide(EARTH_RADIUS).add(1);
  1047.         return tmp.multiply(tmp).reciprocal().multiply(Constants.G0_STANDARD_GRAVITY);
  1048.     }

  1050.     /** Compute semi-annual variation (delta log(rho)).
  1051.      * @param doy day of year
  1052.      * @param alt height (km)
  1053.      * @param f10B average 81-day centered f10
  1054.      * @param s10B average 81-day centered s10
  1055.      * @param xm10B average 81-day centered xn10
  1056.      * @return semi-annual variation
  1057.      */
  1058.     private static double semian08(final double doy, final double alt,
  1059.                                    final double f10B, final double s10B, final double xm10B) {

  1061.         final double htz = alt / 1000.0;

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

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

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

  1072.         // SEMIANNUAL PHASE FUNCTION
  1073.         final double tau   = MathUtils.TWO_PI * (doy - 1.0) / 365;
  1074.         final SinCos sc1P = FastMath.sinCos(tau);
  1075.         final SinCos sc2P = SinCos.sum(sc1P, sc1P);
  1076.         final double gtz = GTM[0] + GTM[1] * sc1P.sin() + GTM[2] * sc1P.cos() + GTM[3] * sc2P.sin() + GTM[4] * sc2P.cos() +
  1077.                    fsmb * (GTM[5] + GTM[6] * sc1P.sin() + GTM[7] * sc1P.cos() + GTM[8] * sc2P.sin() + GTM[9] * sc2P.cos());

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

  1081.     }

  1083.     /** Compute semi-annual variation (delta log(rho)).
  1084.      * @param doy day of year
  1085.      * @param alt height (km)
  1086.      * @param f10B average 81-day centered f10
  1087.      * @param s10B average 81-day centered s10
  1088.      * @param xm10B average 81-day centered xn10
  1089.      * @return semi-annual variation
  1090.      * @param <T> type of the field elements
  1091.      */
  1092.     private static <T extends CalculusFieldElement<T>>  T semian08(final T doy, final T alt,
  1093.                                                                final double f10B, final double s10B, final double xm10B) {

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

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

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

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

  1106.         // SEMIANNUAL PHASE FUNCTION
  1107.         final T tau   = doy.subtract(1).divide(365).multiply(MathUtils.TWO_PI);
  1108.         final FieldSinCos<T> sc1P = FastMath.sinCos(tau);
  1109.         final FieldSinCos<T> sc2P = FieldSinCos.sum(sc1P, sc1P);
  1110.         final T gtz =           sc2P.cos().multiply(GTM[9]).
  1111.                             add(sc2P.sin().multiply(GTM[8])).
  1112.                             add(sc1P.cos().multiply(GTM[7])).
  1113.                             add(sc1P.sin().multiply(GTM[6])).
  1114.                             add(GTM[5]).multiply(fsmb).
  1115.                         add(    sc2P.cos().multiply(GTM[4]).
  1116.                             add(sc2P.sin().multiply(GTM[3])).
  1117.                             add(sc1P.cos().multiply(GTM[2])).
  1118.                             add(sc1P.sin().multiply(GTM[1])).
  1119.                             add(GTM[0]));

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

  1123.     }

  1125.     /** Compute day of year.
  1126.      * @param dateMJD Modified Julian date
  1127.      * @return the number days in year
  1128.      */
  1129.     private static double dayOfYear(final double dateMJD) {
  1130.         final double d1950 = dateMJD - 33281;

  1132.         int iyday = (int) d1950;
  1133.         final double frac = d1950 - iyday;
  1134.         iyday = iyday + 364;

  1136.         int itemp = iyday / 1461;

  1138.         iyday = iyday - itemp * 1461;
  1139.         itemp = iyday / 365;
  1140.         if (itemp >= 3) {
  1141.             itemp = 3;
  1142.         }
  1143.         iyday = iyday - 365 * itemp + 1;
  1144.         return iyday + frac;
  1145.     }

  1147.     /** Compute day of year.
  1148.      * @param dateMJD Modified Julian date
  1149.      * @param <T> type of the field elements
  1150.      * @return the number days in year
  1151.      */
  1152.     private static <T extends CalculusFieldElement<T>> T dayOfYear(final T dateMJD) {
  1153.         final T d1950 = dateMJD.subtract(33281);

  1155.         int iyday = (int) d1950.getReal();
  1156.         final T frac = d1950.subtract(iyday);
  1157.         iyday = iyday + 364;

  1159.         int itemp = iyday / 1461;

  1161.         iyday = iyday - itemp * 1461;
  1162.         itemp = iyday / 365;
  1163.         if (itemp >= 3) {
  1164.             itemp = 3;
  1165.         }
  1166.         iyday = iyday - 365 * itemp + 1;
  1167.         return frac.add(iyday);
  1168.     }

  1170.     // OUTPUT:

  1172.     /** Compute min of two values, one double and one field element.
  1173.      * @param d double value
  1174.      * @param f field element
  1175.      * @param <T> type of the field elements
  1176.      * @return min value
  1177.      */
  1178.     private <T extends CalculusFieldElement<T>> T min(final double d, final T f) {
  1179.         return (f.getReal() > d) ? f.getField().getZero().newInstance(d) : f;
  1180.     }

  1182.     /** Compute max of two values, one double and one field element.
  1183.      * @param d double value
  1184.      * @param f field element
  1185.      * @param <T> type of the field elements
  1186.      * @return max value
  1187.      */
  1188.     private <T extends CalculusFieldElement<T>> T max(final double d, final T f) {
  1189.         return (f.getReal() <= d) ? f.getField().getZero().newInstance(d) : f;
  1190.     }

  1192.     /** Get the local density.
  1193.      * @param date current date
  1194.      * @param position current position in frame
  1195.      * @param frame the frame in which is defined the position
  1196.      * @return local density (kg/m³)
  1197.      */
  1198.     public double getDensity(final AbsoluteDate date, final Vector3D position,
  1199.                              final Frame frame) {
  1200.         // check if data are available :
  1201.         if (date.compareTo(inputParams.getMaxDate()) > 0 ||
  1202.             date.compareTo(inputParams.getMinDate()) < 0) {
  1203.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  1204.                                       date, inputParams.getMinDate(), inputParams.getMaxDate());
  1205.         }

  1207.         // compute MJD date
  1208.         final DateTimeComponents dt = date.getComponents(utc);
  1209.         final double dateMJD = dt.getDate().getMJD() +
  1210.                 dt.getTime().getSecondsInLocalDay() / Constants.JULIAN_DAY;

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

  1215.         // compute sun position
  1216.         final Frame ecef = earth.getBodyFrame();
  1217.         final Vector3D sunPos = getSunPosition(date, ecef);
  1218.         final GeodeticPoint sunInBody = earth.transform(sunPos, ecef, date);

  1220.         return getDensity(dateMJD,
  1221.                           sunInBody.getLongitude(), sunInBody.getLatitude(),
  1222.                           inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude(),
  1223.                           inputParams.getF10(date), inputParams.getF10B(date),
  1224.                           inputParams.getS10(date), inputParams.getS10B(date),
  1225.                           inputParams.getXM10(date), inputParams.getXM10B(date),
  1226.                           inputParams.getY10(date), inputParams.getY10B(date),
  1227.                           inputParams.getDSTDTC(date));

  1229.     }

  1231.     /** {@inheritDoc} */
  1232.     @Override
  1233.     public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
  1234.                                                         final FieldVector3D<T> position,
  1235.                                                         final Frame frame) {

  1237.         // check if data are available :
  1238.         final AbsoluteDate dateD = date.toAbsoluteDate();
  1239.         if (dateD.compareTo(inputParams.getMaxDate()) > 0 ||
  1240.                         dateD.compareTo(inputParams.getMinDate()) < 0) {
  1241.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  1242.                                       dateD, inputParams.getMinDate(), inputParams.getMaxDate());
  1243.         }

  1245.         // compute MJD date
  1246.         final DateTimeComponents components = date.getComponents(utc);
  1247.         final T dateMJD = date
  1248.                 .durationFrom(new FieldAbsoluteDate<>(date.getField(), components, utc))
  1249.                 .add(components.getTime().getSecondsInLocalDay())
  1250.                 .divide(Constants.JULIAN_DAY)
  1251.                 .add(components.getDate().getMJD());

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

  1256.         // compute sun position
  1257.         final Frame ecef = earth.getBodyFrame();
  1258.         final FieldVector3D<T> sunPos = getSunPosition(date, ecef);
  1259.         final FieldGeodeticPoint<T> sunInBody = earth.transform(sunPos, ecef, date);

  1261.         return getDensity(dateMJD,
  1262.                           sunInBody.getLongitude(), sunInBody.getLatitude(),
  1263.                           inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude(),
  1264.                           inputParams.getF10(dateD), inputParams.getF10B(dateD),
  1265.                           inputParams.getS10(dateD), inputParams.getS10B(dateD),
  1266.                           inputParams.getXM10(dateD), inputParams.getXM10B(dateD),
  1267.                           inputParams.getY10(dateD), inputParams.getY10B(dateD),
  1268.                           inputParams.getDSTDTC(dateD));

  1270.     }

  1272. }
Created with JaCoCo 0.8.12.202403310830