AbstractJacchiaBowmanModel.java
/* Copyright 2002-2025 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.models.earth.atmosphere;
import java.util.stream.IntStream;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.MathUtils;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateTimeComponents;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.utils.Constants;
import org.orekit.utils.ExtendedPositionProvider;
/**
* Base class for Jacchia-Bowman atmospheric models.
* @see JB2006
* @see JB2008
* @author Bruce R Bowman (HQ AFSPC, Space Analysis Division), Feb 2006: FORTRAN routine
* @author Fabien Maussion (java translation)
* @author Bryan Cazabonne (field elements translation)
* @since 13.1
*/
public abstract class AbstractJacchiaBowmanModel extends AbstractSunInfluencedAtmosphere {
/** Minimum altitude (m) for JB models use. */
private static final double ALT_MIN = 90000.;
/** The alpha are the thermal diffusion coefficients in equation (6). */
private static final double[] ALPHA = {0, 0, 0, 0, -0.38};
/** Natural logarithm of 10.0. */
private static final double LOG10 = FastMath.log(10.);
/** Molecular weights in order: N2, O2, O, Ar, He and H. */
private static final double[] AMW = {28.0134, 31.9988, 15.9994, 39.9480, 4.0026, 1.00797};
/** Avogadro's number in mks units (molecules/kmol). */
private static final double AVOGAD = 6.02257e26;
/** The FRAC are the assumed sea-level volume fractions in order: N2, O2, Ar, and He. */
private static final double[] FRAC = {0.78110, 0.20955, 9.3400e-3, 1.2890e-5};
/** Universal gas-constant in mks units (joules/K/kmol). */
private static final double RSTAR = 8.31432;
/** Value used to establish height step sizes in the regime 90km to 105km. */
private static final double R1 = 0.010;
/** Value used to establish height step sizes in the regime 105km to 500km. */
private static final double R2 = 0.025;
/** Value used to establish height step sizes in the regime above 500km. */
private static final double R3 = 0.075;
/** Weights for the Newton-Cotes five-points quadrature formula. */
private static final double[] WT = {0.311111111111111, 1.422222222222222, 0.533333333333333, 1.422222222222222, 0.311111111111111};
/** Earth radius (km). */
private static final double EARTH_RADIUS = 6356.766;
/** DTC relative data. */
private static final double[] BDT_SUB = {-0.457512297e+01, -0.512114909e+01, -0.693003609e+02,
0.203716701e+03, 0.703316291e+03, -0.194349234e+04,
0.110651308e+04, -0.174378996e+03, 0.188594601e+04,
-0.709371517e+04, 0.922454523e+04, -0.384508073e+04,
-0.645841789e+01, 0.409703319e+02, -0.482006560e+03,
0.181870931e+04, -0.237389204e+04, 0.996703815e+03,
0.361416936e+02};
/** DTC relative data. */
private static final double[] CDT_SUB = {-0.155986211e+02, -0.512114909e+01, -0.693003609e+02,
0.203716701e+03, 0.703316291e+03, -0.194349234e+04,
0.110651308e+04, -0.220835117e+03, 0.143256989e+04,
-0.318481844e+04, 0.328981513e+04, -0.135332119e+04,
0.199956489e+02, -0.127093998e+02, 0.212825156e+02,
-0.275555432e+01, 0.110234982e+02, 0.148881951e+03,
-0.751640284e+03, 0.637876542e+03, 0.127093998e+02,
-0.212825156e+02, 0.275555432e+01};
/** Mbar polynomial coeffcients. */
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};
/** Coefficients for high altitude density correction. */
private static final double[] CHT = {0.22, -0.20e-02, 0.115e-02, -0.211e-05};
/** UTC time scale. */
private final TimeScale utc;
/** Earth body shape. */
private final BodyShape earth;
/** Earliest epoch of solar activity data. */
private final AbsoluteDate minDataEpoch;
/** Latest epoch of solar activity data. */
private final AbsoluteDate maxDataEpoch;
/**
* Constructor.
*
* @param sun position provider.
* @param utc UTC time scale. Used to computed the day fraction.
* @param earth the earth body shape
* @param minDataEpoch earliest epoch of solar activity data
* @param maxDataEpoch latest epoch of solar activity data
*/
protected AbstractJacchiaBowmanModel(final ExtendedPositionProvider sun, final TimeScale utc, final BodyShape earth,
final AbsoluteDate minDataEpoch, final AbsoluteDate maxDataEpoch) {
super(sun);
this.utc = utc;
this.earth = earth;
this.minDataEpoch = minDataEpoch;
this.maxDataEpoch = maxDataEpoch;
}
/** Get the UTC time scale.
* @return UTC time scale
*/
public TimeScale getUtc() {
return utc;
}
/** Get the Earth body shape.
* @return the Earth body shape
*/
public BodyShape getEarth() {
return earth;
}
/** {@inheritDoc}*/
@Override
public double getDensity(final AbsoluteDate date, final Vector3D position, final Frame frame) {
// Verify availability of data
if (date.compareTo(maxDataEpoch) > 0 || date.compareTo(minDataEpoch) < 0) {
throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, date, minDataEpoch, maxDataEpoch);
}
// compute geodetic position
final GeodeticPoint inBody = earth.transform(position, frame, date);
// compute sun position
final Frame ecef = getFrame();
final Vector3D sunPos = getSunPosition(date, ecef);
final GeodeticPoint sunInBody = earth.transform(sunPos, ecef, date);
return computeDensity(date, sunInBody.getLongitude(), sunInBody.getLatitude(), inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude());
}
/** {@inheritDoc}*/
@Override
public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position, final Frame frame) {
// Verify availability of data
final AbsoluteDate dateD = date.toAbsoluteDate();
if (dateD.compareTo(maxDataEpoch) > 0 || dateD.compareTo(minDataEpoch) < 0) {
throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, dateD, minDataEpoch, maxDataEpoch);
}
// compute geodetic position (km and °)
final FieldGeodeticPoint<T> inBody = earth.transform(position, frame, date);
// compute sun position
final Frame ecef = getFrame();
final FieldVector3D<T> sunPos = getSunPosition(date, ecef);
final FieldGeodeticPoint<T> sunInBody = earth.transform(sunPos, ecef, date);
return computeDensity(date,
sunInBody.getLongitude(), sunInBody.getLatitude(),
inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude());
}
/** {@inheritDoc} */
@Override
public Frame getFrame() {
return earth.getBodyFrame();
}
/** Computes the local density with initial entries.
* @param date computation epoch
* @param sunRA Right Ascension of Sun (radians)
* @param sunDecli Declination of Sun (radians)
* @param satLon Right Ascension of position (radians)
* @param satLat Geocentric latitude of position (radians)
* @param satAlt Height of position (m)
* @return total mass-Density at input position (kg/m³)
*/
protected double computeDensity(final AbsoluteDate date,
final double sunRA, final double sunDecli,
final double satLon, final double satLat, final double satAlt) {
if (satAlt < ALT_MIN) {
throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, satAlt, ALT_MIN);
}
final double altKm = satAlt / 1000.0;
final DateTimeComponents dt = date.getComponents(utc);
final double dateMJD = dt.getDate().getMJD() +
dt.getTime().getSecondsInLocalDay() / Constants.JULIAN_DAY;
// Equation (14)
// Temperature equation obtained using numerous satellites for the years from 1996 through 2004 when all new solar indices were available
final double tsubc = computeTc(date);
// Equation (15)
final double eta = 0.5 * FastMath.abs(satLat - sunDecli);
final double theta = 0.5 * FastMath.abs(satLat + sunDecli);
// Equation (16)
final double h = satLon - sunRA;
final double tau = h - 0.64577182 + 0.10471976 * FastMath.sin(h + 0.75049158);
final double solarTime = solarTimeHour(h);
// Equation (17)
final double tsubl = tSubL(eta, theta, tau, tsubc);
// Compute correction to dTc for local solar time and lat correction
final double dtclst = dTc(getF10(date), solarTime, satLat, altKm);
// Compute the local exospheric temperature.
final double tInf = computeTInf(date, tsubl, dtclst);
// Equation (9)
final double tsubx = tSubX(tInf);
// Equation (11)
final double gsubx = gSubX(tsubx);
// The TC array will be an argument in the call to "localTemp"
final double[] tc = tSubCArray(tsubx, gsubx, tInf);
// Equation (5)
final double z1 = 90.;
final double z2 = FastMath.min(altKm, 105.0);
double al = FastMath.log(z2 / z1);
int n = (int) FastMath.floor(al / R1) + 1;
double zr = FastMath.exp(al / n);
final double mb1 = mBar(z1);
final double tloc1 = localTemp(z1, tc);
double zend = z1;
double sub2 = 0.;
double ain = mb1 * gravity(z1) / tloc1;
double mb2 = 0;
double tloc2 = 0;
double z = 0;
double gravl = 0;
for (int i = 0; i < n; ++i) {
z = zend;
zend = zr * z;
final double dz = 0.25 * (zend - z);
double sum1 = WT[0] * ain;
for (int j = 1; j < 5; ++j) {
z += dz;
mb2 = mBar(z);
tloc2 = localTemp(z, tc);
gravl = gravity(z);
ain = mb2 * gravl / tloc2;
sum1 += WT[j] * ain;
}
sub2 += dz * sum1;
}
double rho = 3.46e-6 * mb2 * tloc1 / FastMath.exp(sub2 / RSTAR) / (mb1 * tloc2);
// Equation (2)
final double anm = AVOGAD * rho;
final double an = anm / mb2;
// Equation (3)
double fact2 = anm / 28.960;
final double[] aln = new double[6];
aln[0] = FastMath.log(FRAC[0] * fact2);
aln[3] = FastMath.log(FRAC[2] * fact2);
aln[4] = FastMath.log(FRAC[3] * fact2);
// Equation (4)
aln[1] = FastMath.log(fact2 * (1. + FRAC[1]) - an);
aln[2] = FastMath.log(2. * (an - fact2));
if (altKm <= 105.0) {
// Put in negligible hydrogen for use in DO-LOOP 13
aln[5] = aln[4] - 25.0;
}
else {
// Equation (6)
al = FastMath.log(FastMath.min(altKm, 500.0) / z);
n = 1 + (int) FastMath.floor(al / R2);
zr = FastMath.exp(al / n);
sub2 = 0.;
ain = gravl / tloc2;
double tloc3 = 0;
for (int i = 0; i < n; ++i) {
z = zend;
zend = zr * z;
final double dz = 0.25 * (zend - z);
double SUM1 = WT[0] * ain;
for (int j = 1; j < 5; ++j) {
z += dz;
tloc3 = localTemp(z, tc);
gravl = gravity(z);
ain = gravl / tloc3;
SUM1 += WT[j] * ain;
}
sub2 += dz * SUM1;
}
al = FastMath.log(FastMath.max(altKm, 500.0) / z);
final double r = (altKm > 500.0) ? R3 : R2;
n = 1 + (int) FastMath.floor(al / r);
zr = FastMath.exp(al / n);
double sum3 = 0.;
double tloc4 = 0;
for (int i = 0; i < n; ++i) {
z = zend;
zend = zr * z;
final double dz = 0.25 * (zend - z);
double sum1 = WT[0] * ain;
for (int j = 1; j < 5; ++j) {
z += dz;
tloc4 = localTemp(z, tc);
gravl = gravity(z);
ain = gravl / tloc4;
sum1 = sum1 + WT[j] * ain;
}
sum3 = sum3 + dz * sum1;
}
final double altr;
final double hSign;
if (altKm <= 500.) {
altr = FastMath.log(tloc3 / tloc2);
fact2 = sub2 / RSTAR;
hSign = 1.0;
}
else {
altr = FastMath.log(tloc4 / tloc2);
fact2 = (sub2 + sum3) / RSTAR;
hSign = -1.0;
}
for (int i = 0; i < 5; ++i) {
aln[i] = aln[i] - (1.0 + ALPHA[i]) * altr - fact2 * AMW[i];
}
// Equation (7) - Note that in CIRA72, AL10T5 = DLOG10(T500)
final double al10t5 = FastMath.log10(tInf);
final double alnh5 = (5.5 * al10t5 - 39.40) * al10t5 + 73.13;
aln[5] = LOG10 * (alnh5 + 6.) + hSign * (FastMath.log(tloc4 / tloc3) + sum3 * AMW[5] / RSTAR);
}
// Equation (24) - J70 Seasonal-Latitudinal Variation
final double dlrsl = dlrsl(altKm, dateMJD, satLat);
// Equation (23) - Computes the semiannual variation
double dlrsa = 0;
if (z < 2000.0) {
// Use new semiannual model DELTA LOG RHO
dlrsa = semian(date, dayOfYear(dateMJD), altKm);
}
// Sum the delta-log-rhos and apply to the number densities.
// In CIRA72 the following equation contains an actual sum,
// namely DLR = AL10 * (DLRGM + DLRSA + DLRSL)
// However, for Jacchia 70, there is no DLRGM or DLRSA.
final double dlr = LOG10 * (dlrsl + dlrsa);
for (int i = 0; i < 6; ++i) {
aln[i] += dlr;
}
// Compute mass-density and mean-molecular-weight and
// convert number density logs from natural to common.
final double sumnm = IntStream.range(0, 6).mapToDouble(i -> FastMath.exp(aln[i]) * AMW[i]).sum();
rho = sumnm / AVOGAD;
// Compute the high altitude exospheric density correction factor
final double fex = densityCorrectionFactor(altKm, getF10B(date));
// Apply the exospheric density correction factor.
rho *= fex;
return rho;
}
/** Computes the local density with initial entries.
* @param date computation epoch
* @param sunRA Right Ascension of Sun (radians)
* @param sunDecli Declination of Sun (radians)
* @param satLon Right Ascension of position (radians)
* @param satLat Geocentric latitude of position (radians)
* @param satAlt Height of position (m)
* @param <T> type of the elements
* @return total mass-Density at input position (kg/m³)
*/
protected <T extends CalculusFieldElement<T>> T computeDensity(final FieldAbsoluteDate<T> date,
final T sunRA, final T sunDecli,
final T satLon, final T satLat, final T satAlt) {
if (satAlt.getReal() < ALT_MIN) {
throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, satAlt, ALT_MIN);
}
final AbsoluteDate dateR = date.toAbsoluteDate();
final DateTimeComponents components = date.getComponents(utc);
final T dateMJD = date.durationFrom(new FieldAbsoluteDate<>(date.getField(), components, utc))
.add(components.getTime().getSecondsInLocalDay())
.divide(Constants.JULIAN_DAY)
.add(components.getDate().getMJD());
final Field<T> field = satAlt.getField();
final T altKm = satAlt.divide(1000.0);
// Equation (14) (Temperature equation)
final double tsubc = computeTc(dateR);
// Equation (15)
final T eta = FastMath.abs(satLat.subtract(sunDecli)).multiply(0.5);
final T theta = FastMath.abs(satLat.add(sunDecli)).multiply(0.5);
// Equation (16)
final T h = satLon.subtract(sunRA);
final T tau = h.subtract(0.64577182).add(h.add(0.75049158).sin().multiply(0.10471976));
final T solarTime = solarTimeHour(h);
// Equation (17)
final T tsubl = tSubL(eta, theta, tau, tsubc);
// Compute correction to dTc for local solar time and lat correction
final T dtclst = dTc(getF10(dateR), solarTime, satLat, altKm);
// Compute the local exospheric temperature.
final T tinf = computeTInf(dateR, tsubl, dtclst);
// Equation (9)
final T tsubx = tSubX(tinf);
// Equation (11)
final T gsubx = gSubX(tsubx);
// The TC array will be an argument in the call to "localTemp"
final T[] tc = tSubCArray(tsubx, gsubx, tinf, field);
// Equation (5)
final T z1 = field.getZero().newInstance(90.);
final T z2 = FastMath.min(altKm, 105.0);
T al = z2.divide(z1).log();
int n = 1 + (int) FastMath.floor(al.getReal() / R1);
T zr = al.divide(n).exp();
final T mb1 = mBar(z1);
final T tloc1 = localTemp(z1, tc);
T zend = z1;
T sub2 = field.getZero();
T ain = mb1.multiply(gravity(z1)).divide(tloc1);
T mb2 = field.getZero();
T tloc2 = field.getZero();
T z = field.getZero();
T gravl = field.getZero();
for (int i = 0; i < n; ++i) {
z = zend;
zend = zr.multiply(z);
final T dz = zend.subtract(z).multiply(0.25);
T sum1 = ain.multiply(WT[0]);
for (int j = 1; j < 5; ++j) {
z = z.add(dz);
mb2 = mBar(z);
tloc2 = localTemp(z, tc);
gravl = gravity(z);
ain = mb2.multiply(gravl).divide(tloc2);
sum1 = sum1.add(ain.multiply(WT[j]));
}
sub2 = sub2.add(dz.multiply(sum1));
}
T rho = mb2.multiply(3.46e-6).multiply(tloc1).divide(sub2.divide(RSTAR).exp().multiply(mb1.multiply(tloc2)));
// Equation (2)
final T anm = rho.multiply(AVOGAD);
final T an = anm.divide(mb2);
// Equation (3)
T fact2 = anm.divide(28.960);
final T[] aln = MathArrays.buildArray(field, 6);
aln[0] = fact2.multiply(FRAC[0]).log();
aln[3] = fact2.multiply(FRAC[2]).log();
aln[4] = fact2.multiply(FRAC[3]).log();
// Equation (4)
aln[1] = fact2.multiply(1. + FRAC[1]).subtract(an).log();
aln[2] = an.subtract(fact2).multiply(2).log();
if (altKm.getReal() <= 105.0) {
// Put in negligible hydrogen for use in DO-LOOP 13
aln[5] = aln[4].subtract(25.0);
}
else {
// Equation (6)
al = FastMath.min(altKm, 500.0).divide(z).log();
n = 1 + (int) FastMath.floor(al.getReal() / R2);
zr = al.divide(n).exp();
sub2 = field.getZero();
ain = gravl.divide(tloc2);
T tloc3 = field.getZero();
for (int i = 0; i < n; ++i) {
z = zend;
zend = zr.multiply(z);
final T dz = zend.subtract(z).multiply(0.25);
T sum1 = ain.multiply(WT[0]);
for (int j = 1; j < 5; ++j) {
z = z.add(dz);
tloc3 = localTemp(z, tc);
gravl = gravity(z);
ain = gravl.divide(tloc3);
sum1 = sum1.add(ain.multiply(WT[j]));
}
sub2 = sub2.add(dz.multiply(sum1));
}
al = FastMath.max(altKm, 500.0).divide(z).log();
final double r = (altKm.getReal() > 500.0) ? R3 : R2;
n = 1 + (int) FastMath.floor(al.getReal() / r);
zr = al.divide(n).exp();
T sum3 = field.getZero();
T tloc4 = field.getZero();
for (int i = 0; i < n; ++i) {
z = zend;
zend = zr.multiply(z);
final T dz = zend.subtract(z).multiply(0.25);
T sum1 = ain.multiply(WT[0]);
for (int j = 1; j < 5; ++j) {
z = z.add(dz);
tloc4 = localTemp(z, tc);
gravl = gravity(z);
ain = gravl.divide(tloc4);
sum1 = sum1.add(ain.multiply(WT[j]));
}
sum3 = sum3.add(dz.multiply(sum1));
}
final T altr;
final double hSign;
if (altKm.getReal() <= 500.) {
altr = tloc3.divide(tloc2).log();
fact2 = sub2.divide(RSTAR);
hSign = 1.0;
}
else {
altr = tloc4.divide(tloc2).log();
fact2 = sub2.add(sum3).divide(RSTAR);
hSign = -1.0;
}
for (int i = 0; i < 5; ++i) {
aln[i] = aln[i].subtract(altr.multiply(1.0 + ALPHA[i])).subtract(fact2.multiply(AMW[i]));
}
// Equation (7) - Note that in CIRA72, AL10T5 = DLOG10(T500)
final T al10t5 = tinf.log10();
final T alnh5 = al10t5.multiply(5.5).subtract(39.40).multiply(al10t5).add(73.13);
aln[5] = alnh5.add(6.).multiply(LOG10).add(tloc4.divide(tloc3).log().add(sum3.multiply(AMW[5] / RSTAR)).multiply(hSign));
}
// Equation (24) - J70 Seasonal-Latitudinal Variation
final T dlrsl = dlrsl(altKm, dateMJD, satLat);
// Equation (23) - Computes the semiannual variation
T dlrsa = field.getZero();
if (z.getReal() < 2000.0) {
// Use new semiannual model DELTA LOG RHO
dlrsa = semian(dateR, dayOfYear(dateMJD), altKm);
}
// Sum the delta-log-rhos and apply to the number densities.
// In CIRA72 the following equation contains an actual sum,
// namely DLR = AL10 * (DLRGM + DLRSA + DLRSL)
// However, for Jacchia 70, there is no DLRGM or DLRSA.
final T dlr = dlrsl.add(dlrsa).multiply(LOG10);
for (int i = 0; i < 6; ++i) {
aln[i] = aln[i].add(dlr);
}
// Compute mass-density and mean-molecular-weight and
// convert number density logs from natural to common.
T sumnm = field.getZero();
for (int i = 0; i < 6; ++i) {
sumnm = sumnm.add(aln[i].exp().multiply(AMW[i]));
}
rho = sumnm.divide(AVOGAD);
// Compute the high altitude exospheric density correction factor
final T fex = densityCorrectionFactor(altKm, getF10B(dateR), field);
// Apply the exospheric density correction factor.
rho = rho.multiply(fex);
return rho;
}
/** Computes the semi-annual variation (delta log(rho)).
* @param date computation epoch
* @param day day of year
* @param altKm height (km)
* @return semi-annual variation
*/
protected abstract double semian(AbsoluteDate date, double day, double altKm);
/**
* Computes the local exospheric temperature.
* @param date computation epoch
* @param tsubl exospheric temperature ("tsubl" term), given by Equation (17)
* @param dtclst correction to dTc for local solar time and lat correction
* @return the local exospheric temperature
*/
protected abstract double computeTInf(AbsoluteDate date, double tsubl, double dtclst);
/** Computes the temperature equation.
* @param date computation epoch
* @return the temperature equation
*/
protected abstract double computeTc(AbsoluteDate date);
/** Get the 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br> (Tabular time 1.0 day earlier).
* @param date computation epoch
* @return the 10.7-cm Solar flux from model input parameters
*/
protected abstract double getF10(AbsoluteDate date);
/** Get the 10.7-cm Solar Flux, averaged 81-day centered on the input time<br> (Tabular time 1.0 day earlier).
* @param date computation epoch
* @return the 10.7-cm Solar Flux, averaged 81-day centered on the input time
*/
protected abstract double getF10B(AbsoluteDate date);
/** Computes the semi-annual variation (delta log(rho)).
* @param date computation epoch
* @param day day of year
* @param altKm height (km)
* @param <T> type of the elements
* @return semi-annual variation
*/
protected abstract <T extends CalculusFieldElement<T>> T semian(AbsoluteDate date, T day, T altKm);
/**
* Computes the local exospheric temperature.
* @param date computation epoch
* @param tsubl exospheric temperature ("tsubl" term), given by Equation (17)
* @param dtclst correction to dTc for local solar time and lat correction
* @param <T> type of the elements
* @return the local exospheric temperature
*/
protected abstract <T extends CalculusFieldElement<T>> T computeTInf(AbsoluteDate date, T tsubl, T dtclst);
/** Evaluates the solar time in hours.
* @param h difference between the Right Ascension of position and the Right Ascension of Sun (radians)
* @return the solar time in hours
*/
private static double solarTimeHour(final double h) {
final double solTimeHour = FastMath.toDegrees(h + FastMath.PI) / 15.0;
if (solTimeHour >= 24) {
return solTimeHour - 24.;
}
if (solTimeHour < 0) {
return solTimeHour + 24.;
}
return solTimeHour;
}
/** Evaluates the solar time in hours.
* @param h difference between the Right Ascension of position and the Right Ascension of Sun (radians)
* @param <T> type of the field elements
* @return the solar time in hours
*/
private static <T extends CalculusFieldElement<T>> T solarTimeHour(final T h) {
final T solarTime = FastMath.toDegrees(h.add(FastMath.PI)).divide(15.0);
if (solarTime.getReal() >= 24) {
return solarTime.subtract(24);
}
if (solarTime.getReal() < 0) {
return solarTime.add(24);
}
return solarTime;
}
/** Evaluates the exospheric temperature ("tsubl" term), Equation (17).
* <p>
* This temperature corresponds to the exospheric temperature in low
* geomagnetic conditions.
* </p>
* @param eta "eta" term, Equation (15)
* @param theta "theta" term, Equation (15)
* @param tau "tau" term, Equation (16)
* @param tsubc exospheric temperature Tc, Equation (14)
* @return Tl temperature computed by Equation (17)
*/
private static double tSubL(final double eta, final double theta,
final double tau, final double tsubc) {
final double cosEta = FastMath.pow(FastMath.cos(eta), 2.5);
final double sinTheta = FastMath.pow(FastMath.sin(theta), 2.5);
final double cosTau = FastMath.abs(FastMath.cos(0.5 * tau));
final double df = sinTheta + (cosEta - sinTheta) * cosTau * cosTau * cosTau;
return tsubc * (1. + 0.31 * df);
}
/** Evaluates exospheric temperature ("tsubl" term), Equation (17).
* <p>
* This temperature corresponds to the exospheric temperature in low
* geomagnetic conditions.
* </p>
* @param eta "eta" term, Equation (15)
* @param theta "theta" term, Equation (15)
* @param tau "tau" term, Equation (16)
* @param tsubc exospheric temperature Tc, Equation (14)
* @param <T> type of the field elements
* @return Tl temperature computed by Equation (17)
*/
private static <T extends CalculusFieldElement<T>> T tSubL(final T eta, final T theta,
final T tau, final double tsubc) {
final T cos = eta.cos();
final T cosEta = cos.square().multiply(cos.sqrt());
final T sin = theta.sin();
final T sinTheta = sin.square().multiply(sin.sqrt());
final T cosTau = tau.multiply(0.5).cos().abs();
final T df = sinTheta.add(cosEta.subtract(sinTheta).multiply(cosTau).multiply(cosTau).multiply(cosTau));
return df.multiply(0.31).add(1).multiply(tsubc);
}
/** Evaluates the inflection temperature ("tsubx" term), Equation (9).
* <p>
* This temperature corresponds to the temperature at the inflexion altitude.
* At this altitude, the temperature profile has an inflection point.
* </p>
* @param tInf local exospheric temperature
* @return Tx temperature at inflection point, computed by Equation (9)
*/
private static double tSubX(final double tInf) {
return 444.3807 + 0.02385 * tInf - 392.8292 * FastMath.exp(-0.0021357 * tInf);
}
/** Evaluates the inflection temperature ("tsubx" term), Equation (9).
* <p>
* This temperature corresponds to the temperature at the inflexion altitude.
* At this altitude, the temperature profile has an inflection point.
* </p>
* @param tInf local exospheric temperature
* @param <T> type of the field elements
* @return Tx temperature at inflection point, computed by Equation (9)
*/
private static <T extends CalculusFieldElement<T>> T tSubX(final T tInf) {
return tInf.multiply(0.02385).add(444.3807).subtract(tInf.multiply(-0.0021357).exp().multiply(392.8292));
}
/** Evaluates the temperature gradient at the inflection point ("gsubx" term), Equation (11).
* @param tSubX temperature at inflection point
* @return the temperature gradient at the inflection point computed by Equation (11)
*/
private static double gSubX(final double tSubX) {
return 0.054285714 * (tSubX - 183.);
}
/** Evaluates the temperature gradient at the inflection point ("gsubx" term), Equation (11).
* @param tSubX temperature at inflection point
* @param <T> type of the field elements
* @return the temperature gradient at the inflection point computed by Equation (11)
*/
private static <T extends CalculusFieldElement<T>> T gSubX(final T tSubX) {
return tSubX.subtract(183.).multiply(0.054285714);
}
/** Evaluates the Tc array.
* <p>
* The Tc array will be used to evaluates {@link #localTemp(double, double[])}.
* </p>
* @param tSubX temperature at inflection point
* @param gSubX the temperature gradient at the inflection point
* @param tInf local exospheric temperature
* @return the Tc array
*/
private static double[] tSubCArray(final double tSubX, final double gSubX,
final double tInf) {
final double[] tc = new double[4];
tc[0] = tSubX;
tc[1] = gSubX;
tc[2] = (tInf - tSubX) / MathUtils.SEMI_PI;
tc[3] = gSubX / tc[2];
return tc;
}
/** Evaluates the Tc array.
* <p>
* The Tc array will be used to evaluates {@link #localTemp(CalculusFieldElement, CalculusFieldElement[])}
* </p>
* @param tSubX temperature at inflection point
* @param gSubX the temperature gradient at the inflection point
* @param tInf local exospheric temperature
* @param field field for the elements
* @param <T> type of the field elements
* @return the Tc array
*/
private static <T extends CalculusFieldElement<T>> T[] tSubCArray(final T tSubX, final T gSubX,
final T tInf, final Field<T> field) {
final T[] tc = MathArrays.buildArray(field, 4);
tc[0] = tSubX;
tc[1] = gSubX;
tc[2] = tInf.subtract(tSubX).divide(MathUtils.SEMI_PI);
tc[3] = gSubX.divide(tc[2]);
return tc;
}
/** Evaluates the local temperature, Equation (10) or (13) depending on altitude.
* @param z altitude
* @param tc tc array
* @return temperature profile
*/
private static double localTemp(final double z, final double[] tc) {
final double dz = z - 125;
if (dz <= 0) {
return ((-9.8204695e-6 * dz - 7.3039742e-4) * dz * dz + 1.0) * dz * tc[1] + tc[0];
}
else {
return tc[0] + tc[2] * FastMath.atan(tc[3] * dz * (1 + 4.5e-6 * FastMath.pow(dz, 2.5)));
}
}
/** Evaluates the local temperature, Equation (10) or (13) depending on altitude.
* @param z altitude
* @param tc tc array
* @return temperature profile
* @param <T> type of the field elements
*/
private static <T extends CalculusFieldElement<T>> T localTemp(final T z, final T[] tc) {
final T dz = z.subtract(125.);
if (dz.getReal() <= 0.) {
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]);
} else {
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]);
}
}
/** Evaluates mean molecualr mass, Equation (1).
* @param z altitude (km)
* @return mean molecular mass
*/
private static double mBar(final double z) {
final double dz = z - 100.;
double amb = CXAMB[6];
for (int i = 5; i >= 0; --i) {
amb = dz * amb + CXAMB[i];
}
return amb;
}
/** Evaluates mean molecualr mass, Equation (1).
* @param z altitude (km)
* @return mean molecular mass
* @param <T> type of the field elements
*/
private static <T extends CalculusFieldElement<T>> T mBar(final T z) {
final T dz = z.subtract(100.);
T amb = z.getField().getZero().newInstance(CXAMB[6]);
for (int i = 5; i >= 0; --i) {
amb = dz.multiply(amb).add(CXAMB[i]);
}
return amb;
}
/** Evaluates the gravity at the altitude, Equation (8).
* @param z altitude (km)
* @return the gravity field (m/s2)
*/
private static double gravity(final double z) {
final double temp = 1.0 + z / EARTH_RADIUS;
return Constants.G0_STANDARD_GRAVITY / (temp * temp);
}
/** Evaluates the gravity at the altitude, Equation (8).
* @param z altitude (km)
* @return the gravity (m/s2)
* @param <T> type of the field elements
*/
private static <T extends CalculusFieldElement<T>> T gravity(final T z) {
final T tmp = z.divide(EARTH_RADIUS).add(1);
return tmp.multiply(tmp).reciprocal().multiply(Constants.G0_STANDARD_GRAVITY);
}
/** Compute day of year.
* @param dateMJD Modified Julian date
* @return the number days in year
*/
private static double dayOfYear(final double dateMJD) {
final double d1950 = dateMJD - 33281;
int iyday = (int) d1950;
final double frac = d1950 - iyday;
iyday = iyday + 364;
int itemp = iyday / 1461;
iyday = iyday - itemp * 1461;
itemp = iyday / 365;
if (itemp >= 3) {
itemp = 3;
}
iyday = iyday - 365 * itemp + 1;
return iyday + frac;
}
/** Compute day of year.
* @param dateMJD Modified Julian date
* @param <T> type of the field elements
* @return the number days in year
*/
private static <T extends CalculusFieldElement<T>> T dayOfYear(final T dateMJD) {
final T d1950 = dateMJD.subtract(33281);
int iyday = (int) d1950.getReal();
final T frac = d1950.subtract(iyday);
iyday = iyday + 364;
int itemp = iyday / 1461;
iyday = iyday - itemp * 1461;
itemp = iyday / 365;
if (itemp >= 3) {
itemp = 3;
}
iyday = iyday - 365 * itemp + 1;
return frac.add(iyday);
}
/** Evaluates the J70 Seasonal-Latitudinal Variation, Equation (24).
* @param altKm satellite altitude (km)
* @param dateMJD Modified Julian date
* @param satLat Geocentric latitude of position (radians)
* @return the J70 Seasonal-Latitudinal Variation
*/
private static double dlrsl(final double altKm, final double dateMJD, final double satLat) {
final double capPhi = ((dateMJD - 36204.0) / 365.2422) % 1;
final int signum = (satLat >= 0) ? 1 : -1;
final double sinLat = FastMath.sin(satLat);
final double hm90 = altKm - 90.;
return 0.02 * hm90 * FastMath.exp(-0.045 * hm90) * signum * FastMath.sin(MathUtils.TWO_PI * capPhi + 1.72) * sinLat * sinLat;
}
/** Evaluates the J70 Seasonal-Latitudinal Variation, Equation (24).
* @param altKm satellite altitude (km)
* @param dateMJD Modified Julian date
* @param satLat Geocentric latitude of position (radians)
* @param <T> type of the elements
* @return the J70 Seasonal-Latitudinal Variation
*/
private static <T extends CalculusFieldElement<T>> T dlrsl(final T altKm, final T dateMJD, final T satLat) {
T capPhi = dateMJD.subtract(36204.0).divide(365.2422);
capPhi = capPhi.subtract(FastMath.floor(capPhi.getReal()));
final int signum = (satLat.getReal() >= 0.) ? 1 : -1;
final T sinLat = satLat.sin();
final T hm90 = altKm.subtract(90.);
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);
}
/** Evaluates the high altitude exospheric density correction factor.
* @param altKm satellite altitude in km
* @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time
* @return the high altitude exospheric density correction factor
*/
private static double densityCorrectionFactor(final double altKm, final double f10B) {
if (altKm >= 1000.0 && altKm < 1500.0) {
final double zeta = (altKm - 1000.) * 0.002;
final double f15c = CHT[0] + CHT[1] * f10B + CHT[2] * 1500.0 + CHT[3] * f10B * 1500.0;
final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
return 1.0 + zeta * zeta * (fex2 + fex3 * zeta);
}
if (altKm >= 1500.0) {
return CHT[0] + CHT[1] * f10B + CHT[2] * altKm + CHT[3] * f10B * altKm;
}
return 1.0;
}
/** Evaluates the high altitude exospheric density correction factor.
* @param altKm satellite altitude in km
* @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time
* @param field field of the elements
* @param <T> type of the field elements
* @return the high altitude exospheric density correction factor
*/
private static <T extends CalculusFieldElement<T>> T densityCorrectionFactor(final T altKm, final double f10B,
final Field<T> field) {
if (altKm.getReal() >= 1000.0 && altKm.getReal() < 1500.0) {
final T zeta = altKm.subtract(1000.).multiply(0.002);
final double f15c = CHT[0] + CHT[1] * f10B + CHT[2] * 1500.0 + CHT[3] * f10B * 1500.0;
final double f15cZeta = (CHT[2] + CHT[3] * f10B) * 500.0;
final double fex2 = 3.0 * f15c - f15cZeta - 3.0;
final double fex3 = f15cZeta - 2.0 * f15c + 2.0;
return field.getOne().add(zeta.multiply(zeta).multiply(zeta.multiply(fex3).add(fex2)));
}
if (altKm.getReal() >= 1500.0) {
return altKm.multiply(CHT[3] * f10B).add(altKm.multiply(CHT[2])).add(CHT[0] + CHT[1] * f10B);
}
return field.getOne();
}
/** Compute daily temperature correction for Jacchia-Bowman model.
* @param f10 solar flux index
* @param solarTime local solar time (hours in [0, 24[)
* @param satLat sat lat (radians)
* @param satAlt height (km)
* @return dTc correction
*/
private static double dTc(final double f10, final double solarTime,
final double satLat, final double satAlt) {
final double st = solarTime / 24.0;
final double cs = FastMath.cos(satLat);
final double fs = (f10 - 100.0) / 100.0;
// Calculates dTc according to height
if (satAlt >= 120 && satAlt <= 200) {
final double dtc200 = poly2CDTC(fs, st, cs);
final double dtc200dz = poly1CDTC(fs, st, cs);
final double cc = 3.0 * dtc200 - dtc200dz;
final double dd = dtc200 - cc;
final double zp = (satAlt - 120.0) / 80.0;
return zp * zp * (cc + dd * zp);
} else if (satAlt > 200.0 && satAlt <= 240.0) {
final double h = (satAlt - 200.0) / 50.0;
return poly1CDTC(fs, st, cs) * h + poly2CDTC(fs, st, cs);
} else if (satAlt > 240.0 && satAlt <= 300.0) {
final double h = 0.8;
final double bb = poly1CDTC(fs, st, cs);
final double aa = bb * h + poly2CDTC(fs, st, cs);
final double p2BDT = poly2BDTC(st);
final double dtc300 = poly1BDTC(fs, st, cs, 3 * p2BDT);
final double dtc300dz = cs * p2BDT;
final double cc = 3.0 * dtc300 - dtc300dz - 3.0 * aa - 2.0 * bb;
final double dd = dtc300 - aa - bb - cc;
final double zp = (satAlt - 240.0) / 60.0;
return aa + zp * (bb + zp * (cc + zp * dd));
} else if (satAlt > 300.0 && satAlt <= 600.0) {
final double h = satAlt / 100.0;
return poly1BDTC(fs, st, cs, h * poly2BDTC(st));
} else if (satAlt > 600.0 && satAlt <= 800.0) {
final double poly2 = poly2BDTC(st);
final double aa = poly1BDTC(fs, st, cs, 6 * poly2);
final double bb = cs * poly2;
final double cc = -(3.0 * aa + 4.0 * bb) / 4.0;
final double dd = (aa + bb) / 4.0;
final double zp = (satAlt - 600.0) / 100.0;
return aa + zp * (bb + zp * (cc + zp * dd));
}
return 0.;
}
/** Compute daily temperature correction for Jacchia-Bowman model.
* @param f10 solar flux index
* @param solarTime local solar time (hours in [0, 24[)
* @param satLat sat lat (radians)
* @param satAlt height (km)
* @param <T> type of the filed elements
* @return dTc correction
*/
private static <T extends CalculusFieldElement<T>> T dTc(final double f10, final T solarTime,
final T satLat, final T satAlt) {
final T st = solarTime.divide(24.0);
final T cs = satLat.cos();
final double fs = (f10 - 100.0) / 100.0;
// Calculates dTc according to height
if (satAlt.getReal() >= 120 && satAlt.getReal() <= 200) {
final T dtc200 = poly2CDTC(fs, st, cs);
final T dtc200dz = poly1CDTC(fs, st, cs);
final T cc = dtc200.multiply(3).subtract(dtc200dz);
final T dd = dtc200.subtract(cc);
final T zp = satAlt.subtract(120.0).divide(80.0);
return zp.multiply(zp).multiply(cc.add(dd.multiply(zp)));
} else if (satAlt.getReal() > 200.0 && satAlt.getReal() <= 240.0) {
final T h = satAlt.subtract(200.0).divide(50.0);
return poly1CDTC(fs, st, cs).multiply(h).add(poly2CDTC(fs, st, cs));
} else if (satAlt.getReal() > 240.0 && satAlt.getReal() <= 300.0) {
final T h = solarTime.getField().getZero().newInstance(0.8);
final T bb = poly1CDTC(fs, st, cs);
final T aa = bb.multiply(h).add(poly2CDTC(fs, st, cs));
final T p2BDT = poly2BDTC(st);
final T dtc300 = poly1BDTC(fs, st, cs, p2BDT.multiply(3));
final T dtc300dz = cs.multiply(p2BDT);
final T cc = dtc300.multiply(3).subtract(dtc300dz).subtract(aa.multiply(3)).subtract(bb.multiply(2));
final T dd = dtc300.subtract(aa).subtract(bb).subtract(cc);
final T zp = satAlt.subtract(240.0).divide(60.0);
return aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
} else if (satAlt.getReal() > 300.0 && satAlt.getReal() <= 600.0) {
final T h = satAlt.divide(100.0);
return poly1BDTC(fs, st, cs, h.multiply(poly2BDTC(st)));
} else if (satAlt.getReal() > 600.0 && satAlt.getReal() <= 800.0) {
final T poly2 = poly2BDTC(st);
final T aa = poly1BDTC(fs, st, cs, poly2.multiply(6));
final T bb = cs.multiply(poly2);
final T cc = aa.multiply(3).add(bb.multiply(4)).divide(-4.0);
final T dd = aa.add(bb).divide(4.0);
final T zp = satAlt.subtract(600.0).divide(100.0);
return aa.add(zp.multiply(bb.add(zp.multiply(cc.add(zp.multiply(dd))))));
}
return satLat.getField().getZero();
}
/** Calculates first polynomial with CDTC array.
* @param fs scaled flux f10
* @param st local solar time in [0, 1[
* @param cs cosine of satLat
* @return the value of the polynomial
*/
private static double poly1CDTC(final double fs, final double st, final double cs) {
return CDT_SUB[0] +
fs * (CDT_SUB[1] + st * (CDT_SUB[2] + st * (CDT_SUB[3] + st * (CDT_SUB[4] + st * (CDT_SUB[5] + st * CDT_SUB[6]))))) +
cs * st * (CDT_SUB[7] + st * (CDT_SUB[8] + st * (CDT_SUB[9] + st * (CDT_SUB[10] + st * CDT_SUB[11])))) +
cs * (CDT_SUB[12] + fs * (CDT_SUB[13] + st * (CDT_SUB[14] + st * CDT_SUB[15])));
}
/** Calculates first polynomial with CDTC array.
* @param fs scaled flux f10
* @param st local solar time in [0, 1[
* @param cs cosine of satLat
* @param <T> type of the field elements
* @return the value of the polynomial
*/
private static <T extends CalculusFieldElement<T>> T poly1CDTC(final double fs, final T st, final T cs) {
return st.multiply(CDT_SUB[6]).
add(CDT_SUB[5]).multiply(st).
add(CDT_SUB[4]).multiply(st).
add(CDT_SUB[3]).multiply(st).
add(CDT_SUB[2]).multiply(st).
add(CDT_SUB[1]).multiply(fs).
add(st.multiply(CDT_SUB[11]).
add(CDT_SUB[10]).multiply(st).
add(CDT_SUB[ 9]).multiply(st).
add(CDT_SUB[ 8]).multiply(st).
add(CDT_SUB[7]).multiply(st).multiply(cs)).
add(st.multiply(CDT_SUB[15]).
add(CDT_SUB[14]).multiply(st).
add(CDT_SUB[13]).multiply(fs).
add(CDT_SUB[12]).multiply(cs)).
add(CDT_SUB[0]);
}
/** Calculates second polynomial with CDTC array.
* @param fs scaled flux f10
* @param st local solar time in [0, 1[
* @param cs cosine of satLat
* @return the value of the polynomial
*/
private static double poly2CDTC(final double fs, final double st, final double cs) {
return CDT_SUB[16] + st * cs * (CDT_SUB[17] + st * (CDT_SUB[18] + st * CDT_SUB[19])) +
fs * cs * (CDT_SUB[20] + st * (CDT_SUB[21] + st * CDT_SUB[22]));
}
/** Calculates second polynomial with CDTC array.
* @param fs scaled flux f10
* @param st local solar time in [0, 1[
* @param cs cosine of satLat
* @param <T> type of the field elements
* @return the value of the polynomial
*/
private static <T extends CalculusFieldElement<T>> T poly2CDTC(final double fs, final T st, final T cs) {
return st.multiply(CDT_SUB[19]).
add(CDT_SUB[18]).multiply(st).
add(CDT_SUB[17]).multiply(cs).multiply(st).
add(st.multiply(CDT_SUB[22]).
add(CDT_SUB[21]).multiply(st).
add(CDT_SUB[20]).multiply(cs).multiply(fs)).
add(CDT_SUB[16]);
}
/** Calculates first polynomial with BDTC array.
* @param fs scaled flux f10
* @param st local solar time in [0, 1[
* @param cs cosine of satLat
* @param hp scaled height * poly2BDTC
* @return the value of the polynomial
*/
private static double poly1BDTC(final double fs, final double st, final double cs, final double hp) {
return BDT_SUB[0] +
fs * (BDT_SUB[1] + st * (BDT_SUB[2] + st * (BDT_SUB[3] + st * (BDT_SUB[4] + st * (BDT_SUB[5] + st * BDT_SUB[6]))))) +
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]);
}
/** Calculates first polynomial with BDTC array.
* @param fs scaled flux f10
* @param st local solar time in [0, 1[
* @param cs cosine of satLat
* @param hp scaled height * poly2BDTC
* @param <T> type of the field elements
* @return the value of the polynomial
*/
private static <T extends CalculusFieldElement<T>> T poly1BDTC(final double fs, final T st, final T cs, final T hp) {
return st.multiply(BDT_SUB[6]).
add(BDT_SUB[5]).multiply(st).
add(BDT_SUB[4]).multiply(st).
add(BDT_SUB[3]).multiply(st).
add(BDT_SUB[2]).multiply(st).
add(BDT_SUB[1]).multiply(fs).
add(st.multiply(BDT_SUB[11]).
add(BDT_SUB[10]).multiply(st).
add(BDT_SUB[ 9]).multiply(st).
add(BDT_SUB[ 8]).multiply(st).
add(BDT_SUB[ 7]).multiply(st).
add(hp).add(BDT_SUB[18]).multiply(cs)).
add(BDT_SUB[0]);
}
/** Calculates second polynomial with BDTC array.
* @param st local solar time in [0, 1[
* @return the value of the polynomial
*/
private static double poly2BDTC(final double st) {
return BDT_SUB[12] + st * (BDT_SUB[13] + st * (BDT_SUB[14] + st * (BDT_SUB[15] + st * (BDT_SUB[16] + st * BDT_SUB[17]))));
}
/** Calculates second polynomial with BDTC array.
* @param st local solar time in [0, 1[
* @param <T> type of the field elements
* @return the value of the polynomial
*/
private static <T extends CalculusFieldElement<T>> T poly2BDTC(final T st) {
return st.multiply(BDT_SUB[17]).
add(BDT_SUB[16]).multiply(st).
add(BDT_SUB[15]).multiply(st).
add(BDT_SUB[14]).multiply(st).
add(BDT_SUB[13]).multiply(st).
add(BDT_SUB[12]);
}
}