AbstractJacchiaBowmanModel.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 java.util.stream.IntStream;
  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.Field;
  22. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  23. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  24. import org.hipparchus.util.FastMath;
  25. import org.hipparchus.util.MathArrays;
  26. import org.hipparchus.util.MathUtils;
  27. import org.orekit.bodies.BodyShape;
  28. import org.orekit.bodies.FieldGeodeticPoint;
  29. import org.orekit.bodies.GeodeticPoint;
  30. import org.orekit.errors.OrekitException;
  31. import org.orekit.errors.OrekitMessages;
  32. import org.orekit.frames.Frame;
  33. import org.orekit.time.AbsoluteDate;
  34. import org.orekit.time.DateTimeComponents;
  35. import org.orekit.time.FieldAbsoluteDate;
  36. import org.orekit.time.TimeScale;
  37. import org.orekit.utils.Constants;
  38. import org.orekit.utils.ExtendedPositionProvider;

  40. /**
  41.  * Base class for Jacchia-Bowman atmospheric models.
  42.  * @see JB2006
  43.  * @see JB2008
  44.  * @author Bruce R Bowman (HQ AFSPC, Space Analysis Division), Feb 2006: FORTRAN routine
  45.  * @author Fabien Maussion (java translation)
  46.  * @author Bryan Cazabonne (field elements translation)
  47.  * @since 13.1
  48.  */
  49. public abstract class AbstractJacchiaBowmanModel extends AbstractSunInfluencedAtmosphere {

  51.     /** Minimum altitude (m) for JB models use. */
  52.     private static final double ALT_MIN = 90000.;

  54.     /** The alpha are the thermal diffusion coefficients in equation (6). */
  55.     private static final double[] ALPHA = {0, 0, 0, 0, -0.38};

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

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

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

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

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

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

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

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

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

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

  87.     /** DTC relative data. */
  88.     private static final double[] BDT_SUB = {-0.457512297e+01, -0.512114909e+01, -0.693003609e+02,
  89.                                              0.203716701e+03, 0.703316291e+03, -0.194349234e+04,
  90.                                              0.110651308e+04, -0.174378996e+03, 0.188594601e+04,
  91.                                              -0.709371517e+04, 0.922454523e+04, -0.384508073e+04,
  92.                                              -0.645841789e+01, 0.409703319e+02, -0.482006560e+03,
  93.                                              0.181870931e+04, -0.237389204e+04, 0.996703815e+03,
  94.                                              0.361416936e+02};

  96.     /** DTC relative data.  */
  97.     private static final double[] CDT_SUB = {-0.155986211e+02, -0.512114909e+01, -0.693003609e+02,
  98.                                              0.203716701e+03, 0.703316291e+03, -0.194349234e+04,
  99.                                              0.110651308e+04, -0.220835117e+03, 0.143256989e+04,
  100.                                              -0.318481844e+04, 0.328981513e+04, -0.135332119e+04,
  101.                                              0.199956489e+02, -0.127093998e+02, 0.212825156e+02,
  102.                                              -0.275555432e+01, 0.110234982e+02, 0.148881951e+03,
  103.                                              -0.751640284e+03, 0.637876542e+03, 0.127093998e+02,
  104.                                              -0.212825156e+02, 0.275555432e+01};

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

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

  112.     /** UTC time scale. */
  113.     private final TimeScale utc;

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

  118.     /** Earliest epoch of solar activity data. */
  119.     private final AbsoluteDate minDataEpoch;

  121.     /** Latest epoch of solar activity data. */
  122.     private final AbsoluteDate maxDataEpoch;

  124.     /**
  125.      * Constructor.
  126.      *
  127.      * @param sun          position provider.
  128.      * @param utc          UTC time scale. Used to computed the day fraction.
  129.      * @param earth        the earth body shape
  130.      * @param minDataEpoch earliest epoch of solar activity data
  131.      * @param maxDataEpoch latest epoch of solar activity data
  132.      */
  133.     protected AbstractJacchiaBowmanModel(final ExtendedPositionProvider sun, final TimeScale utc, final BodyShape earth,
  134.                                          final AbsoluteDate minDataEpoch, final AbsoluteDate maxDataEpoch) {
  135.         super(sun);
  136.         this.utc          = utc;
  137.         this.earth        = earth;
  138.         this.minDataEpoch = minDataEpoch;
  139.         this.maxDataEpoch = maxDataEpoch;
  140.     }

  142.     /** Get the UTC time scale.
  143.      * @return UTC time scale
  144.      */
  145.     public TimeScale getUtc() {
  146.         return utc;
  147.     }

  149.     /** Get the Earth body shape.
  150.      * @return the Earth body shape
  151.      */
  152.     public BodyShape getEarth() {
  153.         return earth;
  154.     }

  156.     /** {@inheritDoc}*/
  157.     @Override
  158.     public double getDensity(final AbsoluteDate date, final Vector3D position, final Frame frame) {

  160.         // Verify availability of data
  161.         if (date.compareTo(maxDataEpoch) > 0 || date.compareTo(minDataEpoch) < 0) {
  162.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, date, minDataEpoch, maxDataEpoch);
  163.         }

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

  168.         // compute sun position
  169.         final Frame ecef = getFrame();
  170.         final Vector3D sunPos = getSunPosition(date, ecef);
  171.         final GeodeticPoint sunInBody = earth.transform(sunPos, ecef, date);
  172.         return computeDensity(date, sunInBody.getLongitude(), sunInBody.getLatitude(), inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude());
  173.     }

  175.     /** {@inheritDoc}*/
  176.     @Override
  177.     public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position, final Frame frame) {

  179.         // Verify availability of data
  180.         final AbsoluteDate dateD = date.toAbsoluteDate();
  181.         if (dateD.compareTo(maxDataEpoch) > 0 || dateD.compareTo(minDataEpoch) < 0) {
  182.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, dateD, minDataEpoch, maxDataEpoch);
  183.         }

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

  188.         // compute sun position
  189.         final Frame ecef = getFrame();
  190.         final FieldVector3D<T> sunPos = getSunPosition(date, ecef);
  191.         final FieldGeodeticPoint<T> sunInBody = earth.transform(sunPos, ecef, date);
  192.         return computeDensity(date,
  193.                               sunInBody.getLongitude(), sunInBody.getLatitude(),
  194.                               inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude());
  195.     }

  197.     /** {@inheritDoc} */
  198.     @Override
  199.     public Frame getFrame() {
  200.         return earth.getBodyFrame();
  201.     }

  203.     /** Computes the local density with initial entries.
  204.      * @param date computation epoch
  205.      * @param sunRA Right Ascension of Sun (radians)
  206.      * @param sunDecli Declination of Sun (radians)
  207.      * @param satLon Right Ascension of position (radians)
  208.      * @param satLat Geocentric latitude of position (radians)
  209.      * @param satAlt Height of position (m)
  210.      * @return total mass-Density at input position (kg/m³)
  211.      */
  212.     protected double computeDensity(final AbsoluteDate date,
  213.                                     final double sunRA, final double sunDecli,
  214.                                     final double satLon, final double satLat, final double satAlt) {

  216.         if (satAlt < ALT_MIN) {
  217.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, satAlt, ALT_MIN);
  218.         }
  219.         final double altKm = satAlt / 1000.0;

  221.         final DateTimeComponents dt = date.getComponents(utc);
  222.         final double dateMJD = dt.getDate().getMJD() +
  223.                                dt.getTime().getSecondsInLocalDay() / Constants.JULIAN_DAY;

  225.         // Equation (14)
  226.         // Temperature equation obtained using numerous satellites for the years from 1996 through 2004 when all new solar indices were available
  227.         final double tsubc = computeTc(date);

  229.         // Equation (15)
  230.         final double eta = 0.5 * FastMath.abs(satLat - sunDecli);
  231.         final double theta = 0.5 * FastMath.abs(satLat + sunDecli);

  233.         // Equation (16)
  234.         final double h   = satLon - sunRA;
  235.         final double tau = h - 0.64577182 + 0.10471976 * FastMath.sin(h + 0.75049158);
  236.         final double solarTime = solarTimeHour(h);

  238.         // Equation (17)
  239.         final double tsubl = tSubL(eta, theta, tau, tsubc);

  241.         // Compute correction to dTc for local solar time and lat correction
  242.         final double dtclst = dTc(getF10(date), solarTime, satLat, altKm);

  244.         // Compute the local exospheric temperature.
  245.         final double tInf = computeTInf(date, tsubl, dtclst);

  247.         // Equation (9)
  248.         final double tsubx = tSubX(tInf);

  250.         // Equation (11)
  251.         final double gsubx = gSubX(tsubx);

  253.         // The TC array will be an argument in the call to "localTemp"
  254.         final double[] tc = tSubCArray(tsubx, gsubx, tInf);

  256.         // Equation (5)
  257.         final double z1    = 90.;
  258.         final double z2    = FastMath.min(altKm, 105.0);
  259.         double al          = FastMath.log(z2 / z1);
  260.         int n              = (int) FastMath.floor(al / R1) + 1;
  261.         double zr          = FastMath.exp(al / n);
  262.         final double mb1   = mBar(z1);
  263.         final double tloc1 = localTemp(z1, tc);
  264.         double zend        = z1;
  265.         double sub2        = 0.;
  266.         double ain         = mb1 * gravity(z1) / tloc1;
  267.         double mb2         = 0;
  268.         double tloc2       = 0;
  269.         double z           = 0;
  270.         double gravl       = 0;

  272.         for (int i = 0; i < n; ++i) {
  273.             z = zend;
  274.             zend = zr * z;
  275.             final double dz = 0.25 * (zend - z);
  276.             double sum1 = WT[0] * ain;
  277.             for (int j = 1; j < 5; ++j) {
  278.                 z += dz;
  279.                 mb2   = mBar(z);
  280.                 tloc2 = localTemp(z, tc);
  281.                 gravl = gravity(z);
  282.                 ain   = mb2 * gravl / tloc2;
  283.                 sum1 += WT[j] * ain;
  284.             }
  285.             sub2 += dz * sum1;
  286.         }
  287.         double rho = 3.46e-6 * mb2 * tloc1 / FastMath.exp(sub2 / RSTAR) / (mb1 * tloc2);

  289.         // Equation (2)
  290.         final double anm = AVOGAD * rho;
  291.         final double an = anm / mb2;

  293.         // Equation (3)
  294.         double fact2 = anm / 28.960;
  295.         final double[] aln = new double[6];
  296.         aln[0] = FastMath.log(FRAC[0] * fact2);
  297.         aln[3] = FastMath.log(FRAC[2] * fact2);
  298.         aln[4] = FastMath.log(FRAC[3] * fact2);
  299.         // Equation (4)
  300.         aln[1] = FastMath.log(fact2 * (1. + FRAC[1]) - an);
  301.         aln[2] = FastMath.log(2. * (an - fact2));

  303.         if (altKm <= 105.0) {
  304.             // Put in negligible hydrogen for use in DO-LOOP 13
  305.             aln[5] = aln[4] - 25.0;
  306.         }
  307.         else {
  308.             // Equation (6)
  309.             al = FastMath.log(FastMath.min(altKm, 500.0) / z);
  310.             n = 1 + (int) FastMath.floor(al / R2);
  311.             zr = FastMath.exp(al / n);
  312.             sub2 = 0.;
  313.             ain = gravl / tloc2;

  315.             double tloc3 = 0;
  316.             for (int i = 0; i < n; ++i) {
  317.                 z = zend;
  318.                 zend = zr * z;
  319.                 final double dz = 0.25 * (zend - z);
  320.                 double SUM1 = WT[0] * ain;
  321.                 for (int j = 1; j < 5; ++j) {
  322.                     z += dz;
  323.                     tloc3 = localTemp(z, tc);
  324.                     gravl = gravity(z);
  325.                     ain = gravl / tloc3;
  326.                     SUM1 += WT[j] * ain;
  327.                 }
  328.                 sub2 += dz * SUM1;
  329.             }

  331.             al = FastMath.log(FastMath.max(altKm, 500.0) / z);
  332.             final double r = (altKm > 500.0) ? R3 : R2;
  333.             n = 1 + (int) FastMath.floor(al / r);
  334.             zr = FastMath.exp(al / n);
  335.             double sum3 = 0.;
  336.             double tloc4 = 0;
  337.             for (int i = 0; i < n; ++i) {
  338.                 z = zend;
  339.                 zend = zr * z;
  340.                 final double dz = 0.25 * (zend - z);
  341.                 double sum1 = WT[0] * ain;
  342.                 for (int j = 1; j < 5; ++j) {
  343.                     z += dz;
  344.                     tloc4 = localTemp(z, tc);
  345.                     gravl = gravity(z);
  346.                     ain = gravl / tloc4;
  347.                     sum1 = sum1 + WT[j] * ain;
  348.                 }
  349.                 sum3 = sum3 + dz * sum1;
  350.             }
  351.             final double altr;
  352.             final double hSign;
  353.             if (altKm <= 500.) {
  354.                 altr = FastMath.log(tloc3 / tloc2);
  355.                 fact2 = sub2 / RSTAR;
  356.                 hSign = 1.0;
  357.             }
  358.             else {
  359.                 altr = FastMath.log(tloc4 / tloc2);
  360.                 fact2 = (sub2 + sum3) / RSTAR;
  361.                 hSign = -1.0;
  362.             }
  363.             for (int i = 0; i < 5; ++i) {
  364.                 aln[i] = aln[i] - (1.0 + ALPHA[i]) * altr - fact2 * AMW[i];
  365.             }

  367.             // Equation (7) - Note that in CIRA72, AL10T5 = DLOG10(T500)
  368.             final double al10t5 = FastMath.log10(tInf);
  369.             final double alnh5 = (5.5 * al10t5 - 39.40) * al10t5 + 73.13;
  370.             aln[5] = LOG10 * (alnh5 + 6.) + hSign * (FastMath.log(tloc4 / tloc3) + sum3 * AMW[5] / RSTAR);
  371.         }

  373.         // Equation (24)  - J70 Seasonal-Latitudinal Variation
  374.         final double dlrsl = dlrsl(altKm, dateMJD, satLat);

  376.         // Equation (23) - Computes the semiannual variation
  377.         double dlrsa = 0;
  378.         if (z < 2000.0) {
  379.             // Use new semiannual model DELTA LOG RHO
  380.             dlrsa = semian(date, dayOfYear(dateMJD), altKm);
  381.         }

  383.         // Sum the delta-log-rhos and apply to the number densities.
  384.         // In CIRA72 the following equation contains an actual sum,
  385.         // namely DLR = AL10 * (DLRGM + DLRSA + DLRSL)
  386.         // However, for Jacchia 70, there is no DLRGM or DLRSA.
  387.         final double dlr = LOG10 * (dlrsl + dlrsa);
  388.         for (int i = 0; i < 6; ++i) {
  389.             aln[i] += dlr;
  390.         }

  392.         // Compute mass-density and mean-molecular-weight and
  393.         // convert number density logs from natural to common.
  394.         final double sumnm = IntStream.range(0, 6).mapToDouble(i -> FastMath.exp(aln[i]) * AMW[i]).sum();
  395.         rho = sumnm / AVOGAD;

  397.         // Compute the high altitude exospheric density correction factor
  398.         final double fex = densityCorrectionFactor(altKm, getF10B(date));

  400.         // Apply the exospheric density correction factor.
  401.         rho *= fex;

  403.         return rho;
  404.     }

  406.     /** Computes the local density with initial entries.
  407.      * @param date computation epoch
  408.      * @param sunRA Right Ascension of Sun (radians)
  409.      * @param sunDecli Declination of Sun (radians)
  410.      * @param satLon Right Ascension of position (radians)
  411.      * @param satLat Geocentric latitude of position (radians)
  412.      * @param satAlt Height of position (m)
  413.      * @param <T> type of the elements
  414.      * @return total mass-Density at input position (kg/m³)
  415.      */
  416.     protected <T extends CalculusFieldElement<T>> T computeDensity(final FieldAbsoluteDate<T> date,
  417.                                                                    final T sunRA, final T sunDecli,
  418.                                                                    final T satLon, final T satLat, final T satAlt) {

  420.         if (satAlt.getReal() < ALT_MIN) {
  421.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, satAlt, ALT_MIN);
  422.         }

  424.         final AbsoluteDate dateR = date.toAbsoluteDate();
  425.         final DateTimeComponents components = date.getComponents(utc);
  426.         final T dateMJD = date.durationFrom(new FieldAbsoluteDate<>(date.getField(), components, utc))
  427.                               .add(components.getTime().getSecondsInLocalDay())
  428.                               .divide(Constants.JULIAN_DAY)
  429.                               .add(components.getDate().getMJD());

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

  434.         // Equation (14) (Temperature equation)
  435.         final double tsubc = computeTc(dateR);

  437.         // Equation (15)
  438.         final T eta   = FastMath.abs(satLat.subtract(sunDecli)).multiply(0.5);
  439.         final T theta = FastMath.abs(satLat.add(sunDecli)).multiply(0.5);

  441.         // Equation (16)
  442.         final T h = satLon.subtract(sunRA);
  443.         final T tau = h.subtract(0.64577182).add(h.add(0.75049158).sin().multiply(0.10471976));
  444.         final T solarTime = solarTimeHour(h);

  446.         // Equation (17)
  447.         final T tsubl = tSubL(eta, theta, tau, tsubc);

  449.         // Compute correction to dTc for local solar time and lat correction
  450.         final T dtclst = dTc(getF10(dateR), solarTime, satLat, altKm);

  452.         // Compute the local exospheric temperature.
  453.         final T tinf = computeTInf(dateR, tsubl, dtclst);

  455.         // Equation (9)
  456.         final T tsubx = tSubX(tinf);

  458.         // Equation (11)
  459.         final T gsubx = gSubX(tsubx);

  461.         // The TC array will be an argument in the call to "localTemp"
  462.         final T[] tc = tSubCArray(tsubx, gsubx, tinf, field);

  464.         // Equation (5)
  465.         final T z1    = field.getZero().newInstance(90.);
  466.         final T z2    = FastMath.min(altKm, 105.0);
  467.         T al          = z2.divide(z1).log();
  468.         int n         = 1 + (int) FastMath.floor(al.getReal() / R1);
  469.         T zr          = al.divide(n).exp();
  470.         final T mb1   = mBar(z1);
  471.         final T tloc1 = localTemp(z1, tc);
  472.         T zend        = z1;
  473.         T sub2        = field.getZero();
  474.         T ain         = mb1.multiply(gravity(z1)).divide(tloc1);
  475.         T mb2         = field.getZero();
  476.         T tloc2       = field.getZero();
  477.         T z           = field.getZero();
  478.         T gravl       = field.getZero();

  480.         for (int i = 0; i < n; ++i) {
  481.             z = zend;
  482.             zend = zr.multiply(z);
  483.             final T dz = zend.subtract(z).multiply(0.25);
  484.             T sum1 = ain.multiply(WT[0]);
  485.             for (int j = 1; j < 5; ++j) {
  486.                 z = z.add(dz);
  487.                 mb2   = mBar(z);
  488.                 tloc2 = localTemp(z, tc);
  489.                 gravl = gravity(z);
  490.                 ain   = mb2.multiply(gravl).divide(tloc2);
  491.                 sum1  = sum1.add(ain.multiply(WT[j]));
  492.             }
  493.             sub2 = sub2.add(dz.multiply(sum1));
  494.         }
  495.         T rho = mb2.multiply(3.46e-6).multiply(tloc1).divide(sub2.divide(RSTAR).exp().multiply(mb1.multiply(tloc2)));

  497.         // Equation (2)
  498.         final T anm = rho.multiply(AVOGAD);
  499.         final T an  = anm.divide(mb2);

  501.         // Equation (3)
  502.         T fact2  = anm.divide(28.960);
  503.         final T[] aln = MathArrays.buildArray(field, 6);
  504.         aln[0] = fact2.multiply(FRAC[0]).log();
  505.         aln[3] = fact2.multiply(FRAC[2]).log();
  506.         aln[4] = fact2.multiply(FRAC[3]).log();

  508.         // Equation (4)
  509.         aln[1] = fact2.multiply(1. + FRAC[1]).subtract(an).log();
  510.         aln[2] = an.subtract(fact2).multiply(2).log();

  512.         if (altKm.getReal() <= 105.0) {
  513.             // Put in negligible hydrogen for use in DO-LOOP 13
  514.             aln[5] = aln[4].subtract(25.0);
  515.         }
  516.         else {
  517.             // Equation (6)
  518.             al   = FastMath.min(altKm, 500.0).divide(z).log();
  519.             n    = 1 + (int) FastMath.floor(al.getReal() / R2);
  520.             zr   = al.divide(n).exp();
  521.             sub2 = field.getZero();
  522.             ain  = gravl.divide(tloc2);

  524.             T tloc3 = field.getZero();
  525.             for (int i = 0; i < n; ++i) {
  526.                 z = zend;
  527.                 zend = zr.multiply(z);
  528.                 final T dz = zend.subtract(z).multiply(0.25);
  529.                 T sum1 = ain.multiply(WT[0]);
  530.                 for (int j = 1; j < 5; ++j) {
  531.                     z     = z.add(dz);
  532.                     tloc3 = localTemp(z, tc);
  533.                     gravl = gravity(z);
  534.                     ain   = gravl.divide(tloc3);
  535.                     sum1  = sum1.add(ain.multiply(WT[j]));
  536.                 }
  537.                 sub2 = sub2.add(dz.multiply(sum1));
  538.             }

  540.             al = FastMath.max(altKm, 500.0).divide(z).log();
  541.             final double r = (altKm.getReal() > 500.0) ? R3 : R2;
  542.             n = 1 + (int) FastMath.floor(al.getReal() / r);
  543.             zr = al.divide(n).exp();
  544.             T sum3 = field.getZero();
  545.             T tloc4 = field.getZero();
  546.             for (int i = 0; i < n; ++i) {
  547.                 z = zend;
  548.                 zend = zr.multiply(z);
  549.                 final T dz = zend.subtract(z).multiply(0.25);
  550.                 T sum1 = ain.multiply(WT[0]);
  551.                 for (int j = 1; j < 5; ++j) {
  552.                     z = z.add(dz);
  553.                     tloc4 = localTemp(z, tc);
  554.                     gravl = gravity(z);
  555.                     ain   = gravl.divide(tloc4);
  556.                     sum1  = sum1.add(ain.multiply(WT[j]));
  557.                 }
  558.                 sum3 = sum3.add(dz.multiply(sum1));
  559.             }
  560.             final T altr;
  561.             final double hSign;
  562.             if (altKm.getReal() <= 500.) {
  563.                 altr = tloc3.divide(tloc2).log();
  564.                 fact2 = sub2.divide(RSTAR);
  565.                 hSign = 1.0;
  566.             }
  567.             else {
  568.                 altr = tloc4.divide(tloc2).log();
  569.                 fact2 = sub2.add(sum3).divide(RSTAR);
  570.                 hSign = -1.0;
  571.             }
  572.             for (int i = 0; i < 5; ++i) {
  573.                 aln[i] = aln[i].subtract(altr.multiply(1.0 + ALPHA[i])).subtract(fact2.multiply(AMW[i]));
  574.             }

  576.             // Equation (7) - Note that in CIRA72, AL10T5 = DLOG10(T500)
  577.             final T al10t5 = tinf.log10();
  578.             final T alnh5 = al10t5.multiply(5.5).subtract(39.40).multiply(al10t5).add(73.13);
  579.             aln[5] = alnh5.add(6.).multiply(LOG10).add(tloc4.divide(tloc3).log().add(sum3.multiply(AMW[5] / RSTAR)).multiply(hSign));
  580.         }

  582.         // Equation (24)  - J70 Seasonal-Latitudinal Variation
  583.         final T dlrsl = dlrsl(altKm, dateMJD, satLat);

  585.         // Equation (23) - Computes the semiannual variation
  586.         T dlrsa = field.getZero();
  587.         if (z.getReal() < 2000.0) {
  588.             // Use new semiannual model DELTA LOG RHO
  589.             dlrsa = semian(dateR, dayOfYear(dateMJD), altKm);
  590.         }

  592.         // Sum the delta-log-rhos and apply to the number densities.
  593.         // In CIRA72 the following equation contains an actual sum,
  594.         // namely DLR = AL10 * (DLRGM + DLRSA + DLRSL)
  595.         // However, for Jacchia 70, there is no DLRGM or DLRSA.
  596.         final T dlr = dlrsl.add(dlrsa).multiply(LOG10);
  597.         for (int i = 0; i < 6; ++i) {
  598.             aln[i] = aln[i].add(dlr);
  599.         }

  601.         // Compute mass-density and mean-molecular-weight and
  602.         // convert number density logs from natural to common.
  603.         T sumnm = field.getZero();
  604.         for (int i = 0; i < 6; ++i) {
  605.             sumnm = sumnm.add(aln[i].exp().multiply(AMW[i]));
  606.         }
  607.         rho = sumnm.divide(AVOGAD);

  609.         // Compute the high altitude exospheric density correction factor
  610.         final T fex = densityCorrectionFactor(altKm, getF10B(dateR), field);

  612.         // Apply the exospheric density correction factor.
  613.         rho = rho.multiply(fex);

  615.         return rho;
  616.     }

  618.     /** Computes the semi-annual variation (delta log(rho)).
  619.      * @param date computation epoch
  620.      * @param day day of year
  621.      * @param altKm height (km)
  622.      * @return semi-annual variation
  623.      */
  624.     protected abstract double semian(AbsoluteDate date, double day, double altKm);

  626.     /**
  627.      * Computes the local exospheric temperature.
  628.      * @param date computation epoch
  629.      * @param tsubl exospheric temperature ("tsubl" term), given by Equation (17)
  630.      * @param dtclst correction to dTc for local solar time and lat correction
  631.      * @return the local exospheric temperature
  632.      */
  633.     protected abstract double computeTInf(AbsoluteDate date, double tsubl, double dtclst);

  635.     /** Computes the temperature equation.
  636.      * @param date computation epoch
  637.      * @return the temperature equation
  638.      */
  639.     protected abstract double computeTc(AbsoluteDate date);

  641.     /** Get the 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br> (Tabular time 1.0 day earlier).
  642.      * @param date computation epoch
  643.      * @return the 10.7-cm Solar flux from model input parameters
  644.      */
  645.     protected abstract double getF10(AbsoluteDate date);

  647.     /** Get the 10.7-cm Solar Flux, averaged 81-day centered on the input time<br> (Tabular time 1.0 day earlier).
  648.      * @param date computation epoch
  649.      * @return the 10.7-cm Solar Flux, averaged 81-day centered on the input time
  650.      */
  651.     protected abstract double getF10B(AbsoluteDate date);

  653.     /** Computes the semi-annual variation (delta log(rho)).
  654.      * @param date computation epoch
  655.      * @param day day of year
  656.      * @param altKm height (km)
  657.      * @param <T> type of the elements
  658.      * @return semi-annual variation
  659.      */
  660.     protected abstract <T extends CalculusFieldElement<T>> T semian(AbsoluteDate date, T day, T altKm);

  662.     /**
  663.      * Computes the local exospheric temperature.
  664.      * @param date computation epoch
  665.      * @param tsubl exospheric temperature ("tsubl" term), given by Equation (17)
  666.      * @param dtclst correction to dTc for local solar time and lat correction
  667.      * @param <T> type of the elements
  668.      * @return the local exospheric temperature
  669.      */
  670.     protected abstract <T extends CalculusFieldElement<T>> T computeTInf(AbsoluteDate date, T tsubl, T dtclst);

  672.     /** Evaluates the solar time in hours.
  673.      * @param h difference between the Right Ascension of position and the Right Ascension of Sun (radians)
  674.      * @return the solar time in hours
  675.      */
  676.     private static double solarTimeHour(final double h) {
  677.         final double solTimeHour = FastMath.toDegrees(h + FastMath.PI) / 15.0;
  678.         if (solTimeHour >= 24) {
  679.             return solTimeHour - 24.;
  680.         }
  681.         if (solTimeHour < 0) {
  682.             return solTimeHour + 24.;
  683.         }
  684.         return solTimeHour;
  685.     }

  687.     /** Evaluates the solar time in hours.
  688.      * @param h difference between the Right Ascension of position and the Right Ascension of Sun (radians)
  689.      * @param <T> type of the field elements
  690.      * @return the solar time in hours
  691.      */
  692.     private static <T extends CalculusFieldElement<T>> T solarTimeHour(final T h) {
  693.         final T solarTime = FastMath.toDegrees(h.add(FastMath.PI)).divide(15.0);
  694.         if (solarTime.getReal() >= 24) {
  695.             return solarTime.subtract(24);
  696.         }
  697.         if (solarTime.getReal() < 0) {
  698.             return solarTime.add(24);
  699.         }
  700.         return solarTime;
  701.     }

  703.     /** Evaluates the exospheric temperature ("tsubl" term), Equation (17).
  704.      * <p>
  705.      * This temperature corresponds to the exospheric temperature in low
  706.      * geomagnetic conditions.
  707.      * </p>
  708.      * @param eta "eta" term, Equation (15)
  709.      * @param theta "theta" term, Equation (15)
  710.      * @param tau "tau" term, Equation (16)
  711.      * @param tsubc exospheric temperature Tc, Equation (14)
  712.      * @return Tl temperature computed by Equation (17)
  713.      */
  714.     private static double tSubL(final double eta, final double theta,
  715.                                final double tau, final double tsubc) {
  716.         final double cosEta   = FastMath.pow(FastMath.cos(eta), 2.5);
  717.         final double sinTheta = FastMath.pow(FastMath.sin(theta), 2.5);
  718.         final double cosTau   = FastMath.abs(FastMath.cos(0.5 * tau));
  719.         final double df       = sinTheta + (cosEta - sinTheta) * cosTau * cosTau * cosTau;
  720.         return tsubc * (1. + 0.31 * df);
  721.     }

  723.     /** Evaluates exospheric temperature ("tsubl" term), Equation (17).
  724.      * <p>
  725.      * This temperature corresponds to the exospheric temperature in low
  726.      * geomagnetic conditions.
  727.      * </p>
  728.      * @param eta "eta" term, Equation (15)
  729.      * @param theta "theta" term, Equation (15)
  730.      * @param tau "tau" term, Equation (16)
  731.      * @param tsubc exospheric temperature Tc, Equation (14)
  732.      * @param <T> type of the field elements
  733.      * @return Tl temperature computed by Equation (17)
  734.      */
  735.     private static <T extends CalculusFieldElement<T>>  T tSubL(final T eta, final T theta,
  736.                                                                final T tau, final double tsubc) {
  737.         final T cos      = eta.cos();
  738.         final T cosEta   = cos.square().multiply(cos.sqrt());
  739.         final T sin      = theta.sin();
  740.         final T sinTheta = sin.square().multiply(sin.sqrt());
  741.         final T cosTau   = tau.multiply(0.5).cos().abs();
  742.         final T df       = sinTheta.add(cosEta.subtract(sinTheta).multiply(cosTau).multiply(cosTau).multiply(cosTau));
  743.         return df.multiply(0.31).add(1).multiply(tsubc);
  744.     }

  746.     /** Evaluates the inflection temperature ("tsubx" term), Equation (9).
  747.      * <p>
  748.      * This temperature corresponds to the temperature at the inflexion altitude.
  749.      * At this altitude, the temperature profile has an inflection point.
  750.      * </p>
  751.      * @param tInf local exospheric temperature
  752.      * @return Tx temperature at inflection point, computed by Equation (9)
  753.      */
  754.     private static double tSubX(final double tInf) {
  755.         return 444.3807 + 0.02385 * tInf - 392.8292 * FastMath.exp(-0.0021357 * tInf);
  756.     }

  758.     /** Evaluates the inflection temperature ("tsubx" term), Equation (9).
  759.      * <p>
  760.      * This temperature corresponds to the temperature at the inflexion altitude.
  761.      * At this altitude, the temperature profile has an inflection point.
  762.      * </p>
  763.      * @param tInf local exospheric temperature
  764.      * @param <T> type of the field elements
  765.      * @return Tx temperature at inflection point, computed by Equation (9)
  766.      */
  767.     private static <T extends CalculusFieldElement<T>>  T tSubX(final T tInf) {
  768.         return tInf.multiply(0.02385).add(444.3807).subtract(tInf.multiply(-0.0021357).exp().multiply(392.8292));
  769.     }

  771.     /** Evaluates the temperature gradient at the inflection point ("gsubx" term), Equation (11).
  772.      * @param tSubX temperature at inflection point
  773.      * @return the temperature gradient at the inflection point computed by Equation (11)
  774.      */
  775.     private static double gSubX(final double tSubX) {
  776.         return 0.054285714 * (tSubX - 183.);
  777.     }

  779.     /** Evaluates the temperature gradient at the inflection point ("gsubx" term), Equation (11).
  780.      * @param tSubX temperature at inflection point
  781.      * @param <T> type of the field elements
  782.      * @return the temperature gradient at the inflection point computed by Equation (11)
  783.      */
  784.     private static <T extends CalculusFieldElement<T>> T gSubX(final T tSubX) {
  785.         return tSubX.subtract(183.).multiply(0.054285714);
  786.     }

  788.     /** Evaluates the Tc array.
  789.      * <p>
  790.      * The Tc array will be used to evaluates {@link #localTemp(double, double[])}.
  791.      * </p>
  792.      * @param tSubX temperature at inflection point
  793.      * @param gSubX the temperature gradient at the inflection point
  794.      * @param tInf local exospheric temperature
  795.      * @return the Tc array
  796.      */
  797.     private static double[] tSubCArray(final double tSubX, final double gSubX,
  798.                                       final double tInf) {
  799.         final double[] tc = new double[4];
  800.         tc[0] = tSubX;
  801.         tc[1] = gSubX;
  802.         tc[2] = (tInf - tSubX) / MathUtils.SEMI_PI;
  803.         tc[3] = gSubX / tc[2];
  804.         return tc;
  805.     }

  807.     /** Evaluates the Tc array.
  808.      * <p>
  809.      * The Tc array will be used to evaluates {@link #localTemp(CalculusFieldElement, CalculusFieldElement[])}
  810.      * </p>
  811.      * @param tSubX temperature at inflection point
  812.      * @param gSubX the temperature gradient at the inflection point
  813.      * @param tInf local exospheric temperature
  814.      * @param field field for the elements
  815.      * @param <T> type of the field elements
  816.      * @return the Tc array
  817.      */
  818.     private static <T extends CalculusFieldElement<T>> T[] tSubCArray(final T tSubX, final T gSubX,
  819.                                                                      final T tInf, final Field<T> field) {
  820.         final T[] tc = MathArrays.buildArray(field, 4);
  821.         tc[0] = tSubX;
  822.         tc[1] = gSubX;
  823.         tc[2] = tInf.subtract(tSubX).divide(MathUtils.SEMI_PI);
  824.         tc[3] = gSubX.divide(tc[2]);
  825.         return tc;
  826.     }

  828.     /** Evaluates the local temperature, Equation (10) or (13) depending on altitude.
  829.      * @param z  altitude
  830.      * @param tc tc array
  831.      * @return temperature profile
  832.      */
  833.     private static double localTemp(final double z, final double[] tc) {
  834.         final double dz = z - 125;
  835.         if (dz <= 0) {
  836.             return ((-9.8204695e-6 * dz - 7.3039742e-4) * dz * dz + 1.0) * dz * tc[1] + tc[0];
  837.         }
  838.         else {
  839.             return tc[0] + tc[2] * FastMath.atan(tc[3] * dz * (1 + 4.5e-6 * FastMath.pow(dz, 2.5)));
  840.         }
  841.     }

  843.     /** Evaluates the local temperature, Equation (10) or (13) depending on altitude.
  844.      * @param z altitude
  845.      * @param tc tc array
  846.      * @return temperature profile
  847.      * @param <T> type of the field elements
  848.      */
  849.     private static <T extends CalculusFieldElement<T>> T localTemp(final T z, final T[] tc) {
  850.         final T dz = z.subtract(125.);
  851.         if (dz.getReal() <= 0.) {
  852.             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]);
  853.         } else {
  854.             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]);
  855.         }
  856.     }

  858.     /** Evaluates mean molecualr mass, Equation (1).
  859.      * @param z altitude (km)
  860.      * @return mean molecular mass
  861.      */
  862.     private static double mBar(final double z) {
  863.         final double dz = z - 100.;
  864.         double amb = CXAMB[6];
  865.         for (int i = 5; i >= 0; --i) {
  866.             amb = dz * amb + CXAMB[i];
  867.         }
  868.         return amb;
  869.     }

  871.     /** Evaluates mean molecualr mass, Equation (1).
  872.      * @param z altitude (km)
  873.      * @return mean molecular mass
  874.      * @param <T> type of the field elements
  875.      */
  876.     private static <T extends CalculusFieldElement<T>>  T mBar(final T z) {
  877.         final T dz = z.subtract(100.);
  878.         T amb = z.getField().getZero().newInstance(CXAMB[6]);
  879.         for (int i = 5; i >= 0; --i) {
  880.             amb = dz.multiply(amb).add(CXAMB[i]);
  881.         }
  882.         return amb;
  883.     }

  885.     /** Evaluates the gravity at the altitude, Equation (8).
  886.      * @param z altitude (km)
  887.      * @return the gravity field (m/s2)
  888.      */
  889.     private static double gravity(final double z) {
  890.         final double temp = 1.0 + z / EARTH_RADIUS;
  891.         return Constants.G0_STANDARD_GRAVITY / (temp * temp);
  892.     }

  894.     /** Evaluates the gravity at the altitude, Equation (8).
  895.      * @param z altitude (km)
  896.      * @return the gravity (m/s2)
  897.      * @param <T> type of the field elements
  898.      */
  899.     private static <T extends CalculusFieldElement<T>> T gravity(final T z) {
  900.         final T tmp = z.divide(EARTH_RADIUS).add(1);
  901.         return tmp.multiply(tmp).reciprocal().multiply(Constants.G0_STANDARD_GRAVITY);
  902.     }

  904.     /** Compute day of year.
  905.      * @param dateMJD Modified Julian date
  906.      * @return the number days in year
  907.      */
  908.     private static double dayOfYear(final double dateMJD) {
  909.         final double d1950 = dateMJD - 33281;

  911.         int iyday = (int) d1950;
  912.         final double frac = d1950 - iyday;
  913.         iyday = iyday + 364;

  915.         int itemp = iyday / 1461;

  917.         iyday = iyday - itemp * 1461;
  918.         itemp = iyday / 365;
  919.         if (itemp >= 3) {
  920.             itemp = 3;
  921.         }
  922.         iyday = iyday - 365 * itemp + 1;
  923.         return iyday + frac;
  924.     }

  926.     /** Compute day of year.
  927.      * @param dateMJD Modified Julian date
  928.      * @param <T> type of the field elements
  929.      * @return the number days in year
  930.      */
  931.     private static <T extends CalculusFieldElement<T>> T dayOfYear(final T dateMJD) {
  932.         final T d1950 = dateMJD.subtract(33281);

  934.         int iyday = (int) d1950.getReal();
  935.         final T frac = d1950.subtract(iyday);
  936.         iyday = iyday + 364;

  938.         int itemp = iyday / 1461;

  940.         iyday = iyday - itemp * 1461;
  941.         itemp = iyday / 365;
  942.         if (itemp >= 3) {
  943.             itemp = 3;
  944.         }
  945.         iyday = iyday - 365 * itemp + 1;
  946.         return frac.add(iyday);
  947.     }

  949.     /** Evaluates the J70 Seasonal-Latitudinal Variation, Equation (24).
  950.      * @param altKm satellite altitude (km)
  951.      * @param dateMJD Modified Julian date
  952.      * @param satLat Geocentric latitude of position (radians)
  953.      * @return the J70 Seasonal-Latitudinal Variation
  954.      */
  955.     private static double dlrsl(final double altKm, final double dateMJD, final double satLat) {
  956.         final double capPhi = ((dateMJD - 36204.0) / 365.2422) % 1;
  957.         final int signum = (satLat >= 0) ? 1 : -1;
  958.         final double sinLat = FastMath.sin(satLat);
  959.         final double hm90  = altKm - 90.;
  960.         return 0.02 * hm90 * FastMath.exp(-0.045 * hm90) * signum * FastMath.sin(MathUtils.TWO_PI * capPhi + 1.72) * sinLat * sinLat;
  961.     }

  963.     /** Evaluates the J70 Seasonal-Latitudinal Variation, Equation (24).
  964.      * @param altKm satellite altitude (km)
  965.      * @param dateMJD Modified Julian date
  966.      * @param satLat Geocentric latitude of position (radians)
  967.      * @param <T> type of the elements
  968.      * @return the J70 Seasonal-Latitudinal Variation
  969.      */
  970.     private static <T extends CalculusFieldElement<T>> T dlrsl(final T altKm, final T dateMJD, final T satLat) {
  971.         T capPhi = dateMJD.subtract(36204.0).divide(365.2422);
  972.         capPhi = capPhi.subtract(FastMath.floor(capPhi.getReal()));
  973.         final int signum = (satLat.getReal() >= 0.) ? 1 : -1;
  974.         final T sinLat = satLat.sin();
  975.         final T hm90  = altKm.subtract(90.);
  976.         return hm90.multiply(0.02).multiply(hm90.multiply(-0.045).exp()).multiply(capPhi.multiply(MathUtils.TWO_PI).add(1.72).sin()).multiply(signum).multiply(sinLat).multiply(sinLat);
  977.     }

  979.     /** Evaluates the high altitude exospheric density correction factor.
  980.      * @param altKm satellite altitude in km
  981.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time
  982.      * @return the high altitude exospheric density correction factor
  983.      */
  984.     private static double densityCorrectionFactor(final double altKm, final double f10B) {
  985.         if (altKm >= 1000.0 && altKm < 1500.0) {
  986.             final double zeta = (altKm - 1000.) * 0.002;
  987.             final double f15c = CHT[0] + CHT[1] * f10B + CHT[2] * 1500.0 + CHT[3] * f10B * 1500.0;
  988.             final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
  989.             final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
  990.             final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
  991.             return 1.0 + zeta * zeta * (fex2 + fex3 * zeta);
  992.         }
  993.         if (altKm >= 1500.0) {
  994.             return CHT[0] + CHT[1] * f10B + CHT[2] * altKm + CHT[3] * f10B * altKm;
  995.         }
  996.         return 1.0;
  997.     }

  999.     /** Evaluates the high altitude exospheric density correction factor.
  1000.      * @param altKm satellite altitude in km
  1001.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time
  1002.      * @param field field of the elements
  1003.      * @param <T> type of the field elements
  1004.      * @return the high altitude exospheric density correction factor
  1005.      */
  1006.     private static <T extends CalculusFieldElement<T>> T densityCorrectionFactor(final T altKm, final double f10B,
  1007.                                                                                 final Field<T> field) {
  1008.         if (altKm.getReal() >= 1000.0 && altKm.getReal() < 1500.0) {
  1009.             final T zeta = altKm.subtract(1000.).multiply(0.002);
  1010.             final double f15c = CHT[0] + CHT[1] * f10B + CHT[2] * 1500.0 + CHT[3] * f10B * 1500.0;
  1011.             final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
  1012.             final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
  1013.             final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
  1014.             return field.getOne().add(zeta.multiply(zeta).multiply(zeta.multiply(fex3).add(fex2)));
  1015.         }
  1016.         if (altKm.getReal() >= 1500.0) {
  1017.             return altKm.multiply(CHT[3] * f10B).add(altKm.multiply(CHT[2])).add(CHT[0] + CHT[1] * f10B);
  1018.         }
  1019.         return field.getOne();
  1020.     }

  1022.     /** Compute daily temperature correction for Jacchia-Bowman model.
  1023.      * @param f10 solar flux index
  1024.      * @param solarTime local solar time (hours in [0, 24[)
  1025.      * @param satLat sat lat (radians)
  1026.      * @param satAlt height (km)
  1027.      * @return dTc correction
  1028.      */
  1029.     private static double dTc(final double f10, final double solarTime,
  1030.                               final double satLat, final double satAlt) {
  1031.         final double st = solarTime / 24.0;
  1032.         final double cs = FastMath.cos(satLat);
  1033.         final double fs = (f10 - 100.0) / 100.0;

  1035.         // Calculates dTc according to height
  1036.         if (satAlt >= 120 && satAlt <= 200) {
  1037.             final double dtc200   = poly2CDTC(fs, st, cs);
  1038.             final double dtc200dz = poly1CDTC(fs, st, cs);
  1039.             final double cc       = 3.0 * dtc200 - dtc200dz;
  1040.             final double dd       = dtc200 - cc;
  1041.             final double zp       = (satAlt - 120.0) / 80.0;
  1042.             return zp * zp * (cc + dd * zp);

  1044.         } else if (satAlt > 200.0 && satAlt <= 240.0) {
  1045.             final double h = (satAlt - 200.0) / 50.0;
  1046.             return poly1CDTC(fs, st, cs) * h + poly2CDTC(fs, st, cs);

  1048.         } else if (satAlt > 240.0 && satAlt <= 300.0) {
  1049.             final double h        = 0.8;
  1050.             final double bb       = poly1CDTC(fs, st, cs);
  1051.             final double aa       = bb * h + poly2CDTC(fs, st, cs);
  1052.             final double p2BDT    = poly2BDTC(st);
  1053.             final double dtc300   = poly1BDTC(fs, st, cs, 3 * p2BDT);
  1054.             final double dtc300dz = cs * p2BDT;
  1055.             final double cc       = 3.0 * dtc300 - dtc300dz - 3.0 * aa - 2.0 * bb;
  1056.             final double dd       = dtc300 - aa - bb - cc;
  1057.             final double zp       = (satAlt - 240.0) / 60.0;
  1058.             return aa + zp * (bb + zp * (cc + zp * dd));

  1060.         } else if (satAlt > 300.0 && satAlt <= 600.0) {
  1061.             final double h = satAlt / 100.0;
  1062.             return poly1BDTC(fs, st, cs, h * poly2BDTC(st));

  1064.         } else if (satAlt > 600.0 && satAlt <= 800.0) {
  1065.             final double poly2 = poly2BDTC(st);
  1066.             final double aa    = poly1BDTC(fs, st, cs, 6 * poly2);
  1067.             final double bb    = cs * poly2;
  1068.             final double cc    = -(3.0 * aa + 4.0 * bb) / 4.0;
  1069.             final double dd    = (aa + bb) / 4.0;
  1070.             final double zp    = (satAlt - 600.0) / 100.0;
  1071.             return aa + zp * (bb + zp * (cc + zp * dd));

  1073.         }
  1074.         return 0.;
  1075.     }

  1077.     /** Compute daily temperature correction for Jacchia-Bowman model.
  1078.      * @param f10 solar flux index
  1079.      * @param solarTime local solar time (hours in [0, 24[)
  1080.      * @param satLat sat lat (radians)
  1081.      * @param satAlt height (km)
  1082.      * @param <T> type of the filed elements
  1083.      * @return dTc correction
  1084.      */
  1085.     private static <T extends CalculusFieldElement<T>> T dTc(final double f10, final T solarTime,
  1086.                                                              final T satLat, final T satAlt) {
  1087.         final T      st = solarTime.divide(24.0);
  1088.         final T      cs = satLat.cos();
  1089.         final double fs = (f10 - 100.0) / 100.0;

  1091.         // Calculates dTc according to height
  1092.         if (satAlt.getReal() >= 120 && satAlt.getReal() <= 200) {
  1093.             final T dtc200   = poly2CDTC(fs, st, cs);
  1094.             final T dtc200dz = poly1CDTC(fs, st, cs);
  1095.             final T cc       = dtc200.multiply(3).subtract(dtc200dz);
  1096.             final T dd       = dtc200.subtract(cc);
  1097.             final T zp       = satAlt.subtract(120.0).divide(80.0);
  1098.             return zp.multiply(zp).multiply(cc.add(dd.multiply(zp)));

  1100.         } else if (satAlt.getReal() > 200.0 && satAlt.getReal() <= 240.0) {
  1101.             final T h = satAlt.subtract(200.0).divide(50.0);
  1102.             return poly1CDTC(fs, st, cs).multiply(h).add(poly2CDTC(fs, st, cs));

  1104.         } else if (satAlt.getReal() > 240.0 && satAlt.getReal() <= 300.0) {
  1105.             final T h        = solarTime.getField().getZero().newInstance(0.8);
  1106.             final T bb       = poly1CDTC(fs, st, cs);
  1107.             final T aa       = bb.multiply(h).add(poly2CDTC(fs, st, cs));
  1108.             final T p2BDT    = poly2BDTC(st);
  1109.             final T dtc300   = poly1BDTC(fs, st, cs, p2BDT.multiply(3));
  1110.             final T dtc300dz = cs.multiply(p2BDT);
  1111.             final T cc       = dtc300.multiply(3).subtract(dtc300dz).subtract(aa.multiply(3)).subtract(bb.multiply(2));
  1112.             final T dd       = dtc300.subtract(aa).subtract(bb).subtract(cc);
  1113.             final T zp       = satAlt.subtract(240.0).divide(60.0);
  1114.             return aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));

  1116.         } else if (satAlt.getReal() > 300.0 && satAlt.getReal() <= 600.0) {
  1117.             final T h = satAlt.divide(100.0);
  1118.             return poly1BDTC(fs, st, cs, h.multiply(poly2BDTC(st)));

  1120.         } else if (satAlt.getReal() > 600.0 && satAlt.getReal() <= 800.0) {
  1121.             final T poly2 = poly2BDTC(st);
  1122.             final T aa    = poly1BDTC(fs, st, cs, poly2.multiply(6));
  1123.             final T bb    = cs.multiply(poly2);
  1124.             final T cc    = aa.multiply(3).add(bb.multiply(4)).divide(-4.0);
  1125.             final T dd    = aa.add(bb).divide(4.0);
  1126.             final T zp    = satAlt.subtract(600.0).divide(100.0);
  1127.             return aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));

  1129.         }
  1130.         return satLat.getField().getZero();
  1131.     }

  1133.     /** Calculates first polynomial with CDTC array.
  1134.      * @param fs scaled flux f10
  1135.      * @param st local solar time in [0, 1[
  1136.      * @param cs cosine of satLat
  1137.      * @return the value of the polynomial
  1138.      */
  1139.     private static double poly1CDTC(final double fs, final double st, final double cs) {
  1140.         return CDT_SUB[0] +
  1141.                fs * (CDT_SUB[1] + st * (CDT_SUB[2] + st * (CDT_SUB[3] + st * (CDT_SUB[4] + st * (CDT_SUB[5] + st * CDT_SUB[6]))))) +
  1142.                cs * st * (CDT_SUB[7] + st * (CDT_SUB[8] + st * (CDT_SUB[9] + st * (CDT_SUB[10] + st * CDT_SUB[11])))) +
  1143.                cs * (CDT_SUB[12] + fs * (CDT_SUB[13] + st * (CDT_SUB[14] + st * CDT_SUB[15])));
  1144.     }

  1146.     /** Calculates first polynomial with CDTC array.
  1147.      * @param fs scaled flux f10
  1148.      * @param st local solar time in [0, 1[
  1149.      * @param cs cosine of satLat
  1150.      * @param <T> type of the field elements
  1151.      * @return the value of the polynomial
  1152.      */
  1153.     private static <T extends CalculusFieldElement<T>>  T poly1CDTC(final double fs, final T st, final T cs) {
  1154.         return st.multiply(CDT_SUB[6]).
  1155.                add(CDT_SUB[5]).multiply(st).
  1156.                add(CDT_SUB[4]).multiply(st).
  1157.                add(CDT_SUB[3]).multiply(st).
  1158.                add(CDT_SUB[2]).multiply(st).
  1159.                add(CDT_SUB[1]).multiply(fs).
  1160.                add(st.multiply(CDT_SUB[11]).
  1161.                add(CDT_SUB[10]).multiply(st).
  1162.                add(CDT_SUB[ 9]).multiply(st).
  1163.                add(CDT_SUB[ 8]).multiply(st).
  1164.                add(CDT_SUB[7]).multiply(st).multiply(cs)).
  1165.                add(st.multiply(CDT_SUB[15]).
  1166.                add(CDT_SUB[14]).multiply(st).
  1167.                add(CDT_SUB[13]).multiply(fs).
  1168.                add(CDT_SUB[12]).multiply(cs)).
  1169.                add(CDT_SUB[0]);
  1170.     }

  1172.     /** Calculates second polynomial with CDTC array.
  1173.      * @param fs scaled flux f10
  1174.      * @param st local solar time in [0, 1[
  1175.      * @param cs cosine of satLat
  1176.      * @return the value of the polynomial
  1177.      */
  1178.     private static double poly2CDTC(final double fs, final double st, final double cs) {
  1179.         return CDT_SUB[16] + st * cs * (CDT_SUB[17] + st * (CDT_SUB[18] + st * CDT_SUB[19])) +
  1180.                fs * cs * (CDT_SUB[20] + st * (CDT_SUB[21] + st * CDT_SUB[22]));
  1181.     }

  1183.     /** Calculates second polynomial with CDTC array.
  1184.      * @param fs scaled flux f10
  1185.      * @param st local solar time in [0, 1[
  1186.      * @param cs cosine of satLat
  1187.      * @param <T> type of the field elements
  1188.      * @return the value of the polynomial
  1189.      */
  1190.     private static <T extends CalculusFieldElement<T>>  T poly2CDTC(final double fs, final T st, final T cs) {
  1191.         return st.multiply(CDT_SUB[19]).
  1192.                add(CDT_SUB[18]).multiply(st).
  1193.                add(CDT_SUB[17]).multiply(cs).multiply(st).
  1194.                add(st.multiply(CDT_SUB[22]).
  1195.                add(CDT_SUB[21]).multiply(st).
  1196.                add(CDT_SUB[20]).multiply(cs).multiply(fs)).
  1197.                add(CDT_SUB[16]);
  1198.     }

  1200.     /** Calculates first polynomial with BDTC array.
  1201.      * @param fs scaled flux f10
  1202.      * @param st local solar time in [0, 1[
  1203.      * @param cs cosine of satLat
  1204.      * @param hp scaled height * poly2BDTC
  1205.      * @return the value of the polynomial
  1206.      */
  1207.     private static double poly1BDTC(final double fs, final double st, final double cs, final double hp) {
  1208.         return BDT_SUB[0] +
  1209.                fs * (BDT_SUB[1] + st * (BDT_SUB[2] + st * (BDT_SUB[3] + st * (BDT_SUB[4] + st * (BDT_SUB[5] + st * BDT_SUB[6]))))) +
  1210.                cs * (st * (BDT_SUB[7] + st * (BDT_SUB[8] + st * (BDT_SUB[9] + st * (BDT_SUB[10] + st * BDT_SUB[11])))) + hp + BDT_SUB[18]);
  1211.     }

  1213.     /** Calculates first polynomial with BDTC array.
  1214.      * @param fs scaled flux f10
  1215.      * @param st local solar time in [0, 1[
  1216.      * @param cs cosine of satLat
  1217.      * @param hp scaled height * poly2BDTC
  1218.      * @param <T> type of the field elements
  1219.      * @return the value of the polynomial
  1220.      */
  1221.     private static <T extends CalculusFieldElement<T>>  T poly1BDTC(final double fs, final T st, final T cs, final T hp) {
  1222.         return st.multiply(BDT_SUB[6]).
  1223.                add(BDT_SUB[5]).multiply(st).
  1224.                add(BDT_SUB[4]).multiply(st).
  1225.                add(BDT_SUB[3]).multiply(st).
  1226.                add(BDT_SUB[2]).multiply(st).
  1227.                add(BDT_SUB[1]).multiply(fs).
  1228.                add(st.multiply(BDT_SUB[11]).
  1229.                add(BDT_SUB[10]).multiply(st).
  1230.                add(BDT_SUB[ 9]).multiply(st).
  1231.                add(BDT_SUB[ 8]).multiply(st).
  1232.                add(BDT_SUB[ 7]).multiply(st).
  1233.                add(hp).add(BDT_SUB[18]).multiply(cs)).
  1234.                add(BDT_SUB[0]);
  1235.     }

  1237.     /** Calculates second polynomial with BDTC array.
  1238.      * @param st local solar time in [0, 1[
  1239.      * @return the value of the polynomial
  1240.      */
  1241.     private static double poly2BDTC(final double st) {
  1242.         return BDT_SUB[12] + st * (BDT_SUB[13] + st * (BDT_SUB[14] + st * (BDT_SUB[15] + st * (BDT_SUB[16] + st * BDT_SUB[17]))));
  1243.     }

  1245.         /** Calculates second polynomial with BDTC array.
  1246.      * @param st local solar time in [0, 1[
  1247.      * @param <T> type of the field elements
  1248.      * @return the value of the polynomial
  1249.      */
  1250.     private static <T extends CalculusFieldElement<T>>  T poly2BDTC(final T st) {
  1251.         return st.multiply(BDT_SUB[17]).
  1252.                add(BDT_SUB[16]).multiply(st).
  1253.                add(BDT_SUB[15]).multiply(st).
  1254.                add(BDT_SUB[14]).multiply(st).
  1255.                add(BDT_SUB[13]).multiply(st).
  1256.                add(BDT_SUB[12]);
  1257.     }

  1259. }
Created with JaCoCo 0.8.12.202403310830