DTM2000.java
/* Copyright 2002-2018 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.forces.drag.atmosphere;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.util.Arrays;
import java.util.Calendar;
import java.util.GregorianCalendar;
import org.hipparchus.RealFieldElement;
import org.hipparchus.exception.DummyLocalizable;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.PVCoordinatesProvider;
/** This atmosphere model is the realization of the DTM-2000 model.
* <p>
* It is described in the paper: <br>
*
* <b>The DTM-2000 empirical thermosphere model with new data assimilation
* and constraints at lower boundary: accuracy and properties</b><br>
*
* <i>S. Bruinsma, G. Thuillier and F. Barlier</i> <br>
*
* Journal of Atmospheric and Solar-Terrestrial Physics 65 (2003) 1053–1070<br>
*
*</p>
*<p>
* This model provides dense output for altitudes beyond 120 km.
*</p>
*
* <p>
* The model needs geographical and time information to compute general values,
* but also needs space weather data : mean and instantaneous solar flux and
* geomagnetic indices.
* </p>
* <p>
* Mean solar flux is (for the moment) represented by the F10.7 indices. Instantaneous
* flux can be set to the mean value if the data is not available. Geomagnetic
* activity is represented by the Kp indice, which goes from 1 (very low activity) to
* 9 (high activity).
*
* <p>
* All these data can be found on the <a href="http://sec.noaa.gov/Data/index.html">
* NOAA (National Oceanic and Atmospheric Administration) website.</a>
* </p>
*
*
* @author R. Biancale, S. Bruinsma: original fortran routine
* @author Fabien Maussion (java translation)
*/
public class DTM2000 implements Atmosphere {
/** Identifier for hydrogen. */
public static final int HYDROGEN = 1;
/** Identifier for helium. */
public static final int HELIUM = 2;
/** Identifier for atomic oxygen. */
public static final int ATOMIC_OXYGEN = 3;
/** Identifier for molecular nitrogen. */
public static final int MOLECULAR_NITROGEN = 4;
/** Identifier for molecular oxygen. */
public static final int MOLECULAR_OXYGEN = 5;
/** Identifier for atomic nitrogen. */
public static final int ATOMIC_NITROGEN = 6;
/** Serializable UID. */
private static final long serialVersionUID = 20170705L;
// Constants :
/** Number of parameters. */
private static final int NLATM = 96;
/** Thermal diffusion coefficient. */
private static final double[] ALEFA = {
0, -0.40, -0.38, 0, 0, 0, 0
};
/** Atomic mass H, He, O, N2, O2, N. */
private static final double[] MA = {
0, 1, 4, 16, 28, 32, 14
};
/** Atomic mass H, He, O, N2, O2, N. */
private static final double[] VMA = {
0, 1.6606e-24, 6.6423e-24, 26.569e-24, 46.4958e-24, 53.1381e-24, 23.2479e-24
};
/** Polar Earth radius. */
private static final double RE = 6356.77;
/** Reference altitude. */
private static final double ZLB0 = 120.0;
/** Cosine of the latitude of the magnetic pole (79N, 71W). */
private static final double CPMG = .19081;
/** Sine of the latitude of the magnetic pole (79N, 71W). */
private static final double SPMG = .98163;
/** Longitude (in radians) of the magnetic pole (79N, 71W). */
private static final double XLMG = -1.2392;
/** Gravity acceleration at 120 km altitude. */
private static final double GSURF = 980.665;
/** Universal gas constant. */
private static final double RGAS = 831.4;
/** 2 * π / 365. */
private static final double ROT = 0.017214206;
/** 2 * rot. */
private static final double ROT2 = 0.034428412;
/** Resources text file. */
private static final String DTM2000 = "/assets/org/orekit/dtm_2000.txt";
// CHECKSTYLE: stop JavadocVariable check
/** Elements coefficients. */
private static double[] tt = null;
private static double[] h = null;
private static double[] he = null;
private static double[] o = null;
private static double[] az2 = null;
private static double[] o2 = null;
private static double[] az = null;
private static double[] t0 = null;
private static double[] tp = null;
/** Sun position. */
private PVCoordinatesProvider sun;
/** External data container. */
private DTM2000InputParameters inputParams;
/** Earth body shape. */
private BodyShape earth;
/** Simple constructor for independent computation.
* @param parameters the solar and magnetic activity data
* @param sun the sun position
* @param earth the earth body shape
* @exception OrekitException if some resource file reading error occurs
*/
public DTM2000(final DTM2000InputParameters parameters,
final PVCoordinatesProvider sun, final BodyShape earth)
throws OrekitException {
synchronized (DTM2000.class) {
// lazy reading of model coefficients
if (tt == null) {
readcoefficients();
}
}
this.earth = earth;
this.sun = sun;
this.inputParams = parameters;
}
/** {@inheritDoc} */
public Frame getFrame() {
return earth.getBodyFrame();
}
/** Get the local density with initial entries.
* @param day day of year
* @param alti altitude in meters
* @param lon local longitude (rad)
* @param lat local latitude (rad)
* @param hl local solar time in rad (O hr = 0 rad)
* @param f instantaneous solar flux (F10.7)
* @param fbar mean solar flux (F10.7)
* @param akp3 3 hrs geomagnetic activity index (1-9)
* @param akp24 Mean of last 24 hrs geomagnetic activity index (1-9)
* @return the local density (kg/m³)
* @exception OrekitException if altitude is outside of supported range
*/
public double getDensity(final int day,
final double alti, final double lon, final double lat,
final double hl, final double f, final double fbar,
final double akp3, final double akp24)
throws OrekitException {
final double threshold = 120000;
if (alti < threshold) {
throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
alti, threshold);
}
final Computation result = new Computation(day, alti / 1000, lon, lat, hl,
new double[] {
0, f, 0
}, new double[] {
0, fbar, 0
}, new double[] {
0, akp3, 0, akp24, 0
});
return result.ro * 1000;
}
/** Get the local density with initial entries.
* @param day day of year
* @param alti altitude in meters
* @param lon local longitude (rad)
* @param lat local latitude (rad)
* @param hl local solar time in rad (O hr = 0 rad)
* @param f instantaneous solar flux (F10.7)
* @param fbar mean solar flux (F10.7)
* @param akp3 3 hrs geomagnetic activity index (1-9)
* @param akp24 Mean of last 24 hrs geomagnetic activity index (1-9)
* @param <T> type of the field elements
* @return the local density (kg/m³)
* @exception OrekitException if altitude is outside of supported range
* @since 9.0
*/
public <T extends RealFieldElement<T>> T getDensity(final int day,
final T alti, final T lon, final T lat,
final T hl, final double f, final double fbar,
final double akp3, final double akp24)
throws OrekitException {
final double threshold = 120000;
if (alti.getReal() < threshold) {
throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
alti, threshold);
}
final FieldComputation<T> result = new FieldComputation<>(day, alti.divide(1000), lon, lat, hl,
new double[] {
0, f, 0
}, new double[] {
0, fbar, 0
}, new double[] {
0, akp3, 0, akp24, 0
});
return result.ro.multiply(1000);
}
/** Store the DTM model elements coefficients in internal arrays.
* @exception OrekitException if some resource file reading error occurs
*/
private static void readcoefficients() throws OrekitException {
final int size = NLATM + 1;
tt = new double[size];
h = new double[size];
he = new double[size];
o = new double[size];
az2 = new double[size];
o2 = new double[size];
az = new double[size];
t0 = new double[size];
tp = new double[size];
final InputStream in = DTM2000.class.getResourceAsStream(DTM2000);
if (in == null) {
throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_RESOURCE, DTM2000);
}
BufferedReader r = null;
try {
r = new BufferedReader(new InputStreamReader(in, "UTF-8"));
r.readLine();
r.readLine();
for (String line = r.readLine(); line != null; line = r.readLine()) {
final int num = Integer.parseInt(line.substring(0, 4).replace(' ', '0'));
line = line.substring(4);
tt[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
h[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
he[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
o[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
az2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
o2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
az[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
t0[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
line = line.substring(13 + 9);
tp[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
}
} catch (IOException ioe) {
throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
} finally {
if (r != null) {
try {
r.close();
} catch (IOException ioe) {
throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
}
}
}
}
/** Get the local density.
* @param date current date
* @param position current position in frame
* @param frame the frame in which is defined the position
* @return local density (kg/m³)
* @exception OrekitException if date is out of range of solar activity model
* or if some frame conversion cannot be performed
*/
public double getDensity(final AbsoluteDate date, final Vector3D position,
final Frame frame)
throws OrekitException {
// check if data are available :
if ((date.compareTo(inputParams.getMaxDate()) > 0) ||
(date.compareTo(inputParams.getMinDate()) < 0)) {
throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
date, inputParams.getMinDate(), inputParams.getMaxDate());
}
// compute day number in current year
final Calendar cal = new GregorianCalendar();
cal.setTime(date.toDate(TimeScalesFactory.getUTC()));
final int day = cal.get(Calendar.DAY_OF_YEAR);
//position in ECEF so we only have to do the transform once
final Frame ecef = earth.getBodyFrame();
final Vector3D pEcef = frame.getTransformTo(ecef, date)
.transformPosition(position);
// compute geodetic position
final GeodeticPoint inBody = earth.transform(pEcef, ecef, date);
final double alti = inBody.getAltitude();
final double lon = inBody.getLongitude();
final double lat = inBody.getLatitude();
// compute local solar time
final Vector3D sunPos = sun.getPVCoordinates(date, ecef).getPosition();
final double hl = FastMath.PI + FastMath.atan2(
sunPos.getX() * pEcef.getY() - sunPos.getY() * pEcef.getX(),
sunPos.getX() * pEcef.getX() + sunPos.getY() * pEcef.getY());
// get current solar activity data and compute
return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(date),
inputParams.getMeanFlux(date), inputParams.getThreeHourlyKP(date),
inputParams.get24HoursKp(date));
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> T
getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position,
final Frame frame)
throws OrekitException {
// check if data are available :
final AbsoluteDate dateD = date.toAbsoluteDate();
if ((dateD.compareTo(inputParams.getMaxDate()) > 0) ||
(dateD.compareTo(inputParams.getMinDate()) < 0)) {
throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
dateD, inputParams.getMinDate(), inputParams.getMaxDate());
}
// compute day number in current year
final Calendar cal = new GregorianCalendar();
cal.setTime(date.toDate(TimeScalesFactory.getUTC()));
final int day = cal.get(Calendar.DAY_OF_YEAR);
//position in ECEF so we only have to do the transform once
final Frame ecef = earth.getBodyFrame();
final FieldVector3D<T> pEcef = frame.getTransformTo(ecef, date).transformPosition(position);
// compute geodetic position
final FieldGeodeticPoint<T> inBody = earth.transform(pEcef, ecef, date);
final T alti = inBody.getAltitude();
final T lon = inBody.getLongitude();
final T lat = inBody.getLatitude();
// compute local solar time
final Vector3D sunPos = sun.getPVCoordinates(dateD, ecef).getPosition();
final T y = pEcef.getY().multiply(sunPos.getX()).subtract(pEcef.getX().multiply(sunPos.getY()));
final T x = pEcef.getX().multiply(sunPos.getX()).add(pEcef.getY().multiply(sunPos.getY()));
final T hl = y.atan2(x).add(FastMath.PI);
// get current solar activity data and compute
return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(dateD),
inputParams.getMeanFlux(dateD), inputParams.getThreeHourlyKP(dateD),
inputParams.get24HoursKp(dateD));
}
/** Local holder for intermediate results ensuring the model is reentrant. */
private static class Computation {
/** Number of days in current year. */
private final int day;
/** Instant solar flux. f[1] = instantaneous flux; f[2] = 0. (not used). */
private final double[] f;
/** Mean solar flux. fbar[1] = mean flux; fbar[2] = 0. (not used). */
private final double[] fbar;
/** Kp coefficients.
* <ul>
* <li>akp[1] = 3-hourly kp</li>
* <li>akp[2] = 0 (not used)</li>
* <li>akp[3] = mean kp of last 24 hours</li>
* <li>akp[4] = 0 (not used)</li>
* </ul>
*/
private final double[] akp;
/** Cosine of the longitude. */
private final double clfl;
/** Sine of the longitude. */
private final double slfl;
/** Total density (g/cm3). */
private final double ro;
// CHECKSTYLE: stop JavadocVariable check
/** Legendre coefficients. */
private final double p10;
private final double p20;
private final double p30;
private final double p40;
private final double p50;
private final double p60;
private final double p11;
private final double p21;
private final double p31;
private final double p41;
private final double p51;
private final double p22;
private final double p32;
private final double p42;
private final double p52;
private final double p62;
private final double p33;
private final double p10mg;
private final double p20mg;
private final double p40mg;
/** Local time intermediate values. */
private final double hl0;
private final double ch;
private final double sh;
private final double c2h;
private final double s2h;
private final double c3h;
private final double s3h;
/** Simple constructor.
* @param day day of year
* @param altiKM altitude <em>in kilometers</em>
* @param lon local longitude (rad)
* @param lat local latitude (rad)
* @param hl local solar time in rad (O hr = 0 rad)
* @param f instantaneous solar flux (F10.7)
* @param fbar mean solar flux (F10.7)
* @param akp geomagnetic activity index
*/
Computation(final int day,
final double altiKM, final double lon, final double lat,
final double hl, final double[] f, final double[] fbar,
final double[] akp) {
this.day = day;
this.f = f;
this.fbar = fbar;
this.akp = akp;
// compute Legendre polynomials wrt geographic pole
final double c = FastMath.sin(lat);
final double c2 = c * c;
final double c4 = c2 * c2;
final double s = FastMath.cos(lat);
final double s2 = s * s;
p10 = c;
p20 = 1.5 * c2 - 0.5;
p30 = c * (2.5 * c2 - 1.5);
p40 = 4.375 * c4 - 3.75 * c2 + 0.375;
p50 = c * (7.875 * c4 - 8.75 * c2 + 1.875);
p60 = (5.5 * c * p50 - 2.5 * p40) / 3.0;
p11 = s;
p21 = 3.0 * c * s;
p31 = s * (7.5 * c2 - 1.5);
p41 = c * s * (17.5 * c2 - 7.5);
p51 = s * (39.375 * c4 - 26.25 * c2 + 1.875);
p22 = 3.0 * s2;
p32 = 15.0 * c * s2;
p42 = s2 * (52.5 * c2 - 7.5);
p52 = 3.0 * c * p42 - 2.0 * p32;
p62 = 2.75 * c * p52 - 1.75 * p42;
p33 = 15.0 * s * s2;
// compute Legendre polynomials wrt magnetic pole (79N, 71W)
final double clmlmg = FastMath.cos(lon - XLMG);
final double cmg = s * CPMG * clmlmg + c * SPMG;
final double cmg2 = cmg * cmg;
final double cmg4 = cmg2 * cmg2;
p10mg = cmg;
p20mg = 1.5 * cmg2 - 0.5;
p40mg = 4.375 * cmg4 - 3.75 * cmg2 + 0.375;
clfl = FastMath.cos(lon);
slfl = FastMath.sin(lon);
// local time
hl0 = hl;
ch = FastMath.cos(hl0);
sh = FastMath.sin(hl0);
c2h = ch * ch - sh * sh;
s2h = 2.0 * ch * sh;
c3h = c2h * ch - s2h * sh;
s3h = s2h * ch + c2h * sh;
final double zlb = ZLB0; // + dzlb ??
final double[] dtt = new double[tt.length];
final double[] dh = new double[tt.length];
final double[] dhe = new double[tt.length];
final double[] dox = new double[tt.length];
final double[] daz2 = new double[tt.length];
final double[] do2 = new double[tt.length];
final double[] daz = new double[tt.length];
final double[] dt0 = new double[tt.length];
final double[] dtp = new double[tt.length];
Arrays.fill(dtt, Double.NaN);
Arrays.fill(dh, Double.NaN);
Arrays.fill(dhe, Double.NaN);
Arrays.fill(dox, Double.NaN);
Arrays.fill(daz2, Double.NaN);
Arrays.fill(do2, Double.NaN);
Arrays.fill(daz, Double.NaN);
Arrays.fill(dt0, Double.NaN);
Arrays.fill(dtp, Double.NaN);
// compute function g(l) / tinf, t120, tp120
int kleq = 1;
final double gdelt = gFunction(tt, dtt, 1, kleq);
dtt[1] = 1.0 + gdelt;
final double tinf = tt[1] * dtt[1];
kleq = 0; // equinox
if ((day < 59) || (day > 284)) {
kleq = -1; // north winter
}
if ((day > 99) && (day < 244)) {
kleq = 1; // north summer
}
final double gdelt0 = gFunction(t0, dt0, 0, kleq);
dt0[1] = (t0[1] + gdelt0) / t0[1];
final double t120 = t0[1] + gdelt0;
final double gdeltp = gFunction(tp, dtp, 0, kleq);
dtp[1] = (tp[1] + gdeltp) / tp[1];
final double tp120 = tp[1] + gdeltp;
// compute n(z) concentrations: H, He, O, N2, O2, N
final double sigma = tp120 / (tinf - t120);
final double dzeta = (RE + zlb) / (RE + altiKM);
final double zeta = (altiKM - zlb) * dzeta;
final double sigzeta = sigma * zeta;
final double expsz = FastMath.exp(-sigzeta);
final double tz = tinf - (tinf - t120) * expsz;
final double[] dbase = new double[7];
kleq = 1;
final double gdelh = gFunction(h, dh, 0, kleq);
dh[1] = FastMath.exp(gdelh);
dbase[1] = h[1] * dh[1];
final double gdelhe = gFunction(he, dhe, 0, kleq);
dhe[1] = FastMath.exp(gdelhe);
dbase[2] = he[1] * dhe[1];
final double gdelo = gFunction(o, dox, 1, kleq);
dox[1] = FastMath.exp(gdelo);
dbase[3] = o[1] * dox[1];
final double gdelaz2 = gFunction(az2, daz2, 1, kleq);
daz2[1] = FastMath.exp(gdelaz2);
dbase[4] = az2[1] * daz2[1];
final double gdelo2 = gFunction(o2, do2, 1, kleq);
do2[1] = FastMath.exp(gdelo2);
dbase[5] = o2[1] * do2[1];
final double gdelaz = gFunction(az, daz, 1, kleq);
daz[1] = FastMath.exp(gdelaz);
dbase[6] = az[1] * daz[1];
final double zlbre = 1.0 + zlb / RE;
final double glb = (GSURF / (zlbre * zlbre)) / (sigma * RGAS * tinf);
final double t120tz = t120 / tz;
// Partial densities in (g/cm3).
// d(1) = hydrogen
// d(2) = helium
// d(3) = atomic oxygen
// d(4) = molecular nitrogen
// d(5) = molecular oxygen
// d(6) = atomic nitrogen
double tmpro = 0;
for (int i = 1; i <= 6; i++) {
final double gamma = MA[i] * glb;
final double upapg = 1.0 + ALEFA[i] + gamma;
final double fzI = FastMath.pow(t120tz, upapg) * FastMath.exp(-sigzeta * gamma);
// concentrations of H, He, O, N2, O2, N (particles/cm³)
final double ccI = dbase[i] * fzI;
// contribution of densities of H, He, O, N2, O2, N (g/cm³)
tmpro += ccI * VMA[i];
}
this.ro = tmpro;
}
/** Computation of function G.
* @param a vector of coefficients for computation
* @param da vector of partial derivatives
* @param ff0 coefficient flag (1 for Ox, Az, He, T°; 0 for H and tp120)
* @param kle_eq season indicator flag (summer, winter, equinox)
* @return value of G
*/
private double gFunction(final double[] a, final double[] da,
final int ff0, final int kle_eq) {
final double[] fmfb = new double[3];
final double[] fbm150 = new double[3];
// latitude terms
da[2] = p20;
da[3] = p40;
da[74] = p10;
double a74 = a[74];
double a77 = a[77];
double a78 = a[78];
if (kle_eq == -1) {
// winter
a74 = -a74;
a77 = -a77;
a78 = -a78;
}
if (kle_eq == 0 ) {
// equinox
a74 = semestrialCorrection(a74);
a77 = semestrialCorrection(a77);
a78 = semestrialCorrection(a78);
}
da[77] = p30;
da[78] = p50;
da[79] = p60;
// flux terms
fmfb[1] = f[1] - fbar[1];
fmfb[2] = f[2] - fbar[2];
fbm150[1] = fbar[1] - 150.0;
fbm150[2] = fbar[2];
da[4] = fmfb[1];
da[6] = fbm150[1];
da[4] = da[4] + a[70] * fmfb[2];
da[6] = da[6] + a[71] * fbm150[2];
da[70] = fmfb[2] * (a[4] + 2.0 * a[5] * da[4] + a[82] * p10 +
a[83] * p20 + a[84] * p30);
da[71] = fbm150[2] * (a[6] + 2.0 * a[69] * da[6] + a[85] * p10 +
a[86] * p20 + a[87] * p30);
da[5] = da[4] * da[4];
da[69] = da[6] * da[6];
da[82] = da[4] * p10;
da[83] = da[4] * p20;
da[84] = da[4] * p30;
da[85] = da[6] * p20;
da[86] = da[6] * p30;
da[87] = da[6] * p40;
// Kp terms
final int ikp = 62;
final int ikpm = 67;
final double c2fi = 1.0 - p10mg * p10mg;
final double dkp = akp[1] + (a[ikp] + c2fi * a[ikp + 1]) * akp[2];
double dakp = a[7] + a[8] * p20mg + a[68] * p40mg +
2.0 * dkp * (a[60] + a[61] * p20mg +
a[75] * 2.0 * dkp * dkp);
da[ikp] = dakp * akp[2];
da[ikp + 1] = da[ikp] * c2fi;
final double dkpm = akp[3] + a[ikpm] * akp[4];
final double dakpm = a[64] + a[65] * p20mg + a[72] * p40mg +
2.0 * dkpm * (a[66] + a[73] * p20mg +
a[76] * 2.0 * dkpm * dkpm);
da[ikpm] = dakpm * akp[4];
da[7] = dkp;
da[8] = p20mg * dkp;
da[68] = p40mg * dkp;
da[60] = dkp * dkp;
da[61] = p20mg * da[60];
da[75] = da[60] * da[60];
da[64] = dkpm;
da[65] = p20mg * dkpm;
da[72] = p40mg * dkpm;
da[66] = dkpm * dkpm;
da[73] = p20mg * da[66];
da[76] = da[66] * da[66];
// non-periodic g(l) function
double f0 = a[4] * da[4] + a[5] * da[5] + a[6] * da[6] +
a[69] * da[69] + a[82] * da[82] + a[83] * da[83] +
a[84] * da[84] + a[85] * da[85] + a[86] * da[86] +
a[87] * da[87];
final double f1f = 1.0 + f0 * ff0;
f0 = f0 + a[2] * da[2] + a[3] * da[3] + a74 * da[74] +
a77 * da[77] + a[7] * da[7] + a[8] * da[8] +
a[60] * da[60] + a[61] * da[61] + a[68] * da[68] +
a[64] * da[64] + a[65] * da[65] + a[66] * da[66] +
a[72] * da[72] + a[73] * da[73] + a[75] * da[75] +
a[76] * da[76] + a78 * da[78] + a[79] * da[79];
// termes annuels symetriques en latitude
da[9] = FastMath.cos(ROT * (day - a[11]));
da[10] = p20 * da[9];
// termes semi-annuels symetriques en latitude
da[12] = FastMath.cos(ROT2 * (day - a[14]));
da[13] = p20 * da[12];
// termes annuels non symetriques en latitude
final double coste = FastMath.cos(ROT * (day - a[18]));
da[15] = p10 * coste;
da[16] = p30 * coste;
da[17] = p50 * coste;
// terme semi-annuel non symetrique en latitude
final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
da[19] = p10 * cos2te;
da[39] = p30 * cos2te;
da[59] = p50 * cos2te;
// termes diurnes [et couples annuel]
da[21] = p11 * ch;
da[22] = p31 * ch;
da[23] = p51 * ch;
da[24] = da[21] * coste;
da[25] = p21 * ch * coste;
da[26] = p11 * sh;
da[27] = p31 * sh;
da[28] = p51 * sh;
da[29] = da[26] * coste;
da[30] = p21 * sh * coste;
// termes semi-diurnes [et couples annuel]
da[31] = p22 * c2h;
da[37] = p42 * c2h;
da[32] = p32 * c2h * coste;
da[33] = p22 * s2h;
da[38] = p42 * s2h;
da[34] = p32 * s2h * coste;
da[88] = p32 * c2h;
da[89] = p32 * s2h;
da[90] = p52 * c2h;
da[91] = p52 * s2h;
double a88 = a[88];
double a89 = a[89];
double a90 = a[90];
double a91 = a[91];
if (kle_eq == -1) { //hiver
a88 = -a88;
a89 = -a89;
a90 = -a90;
a91 = -a91;
}
if (kle_eq == 0) { //equinox
a88 = semestrialCorrection(a88);
a89 = semestrialCorrection(a89);
a90 = semestrialCorrection(a90);
a91 = semestrialCorrection(a91);
}
da[92] = p62 * c2h;
da[93] = p62 * s2h;
// termes ter-diurnes
da[35] = p33 * c3h;
da[36] = p33 * s3h;
// fonction g[l] periodique
double fp = a[9] * da[9] + a[10] * da[10] + a[12] * da[12] + a[13] * da[13] +
a[15] * da[15] + a[16] * da[16] + a[17] * da[17] + a[19] * da[19] +
a[21] * da[21] + a[22] * da[22] + a[23] * da[23] + a[24] * da[24] +
a[25] * da[25] + a[26] * da[26] + a[27] * da[27] + a[28] * da[28] +
a[29] * da[29] + a[30] * da[30] + a[31] * da[31] + a[32] * da[32] +
a[33] * da[33] + a[34] * da[34] + a[35] * da[35] + a[36] * da[36] +
a[37] * da[37] + a[38] * da[38] + a[39] * da[39] + a[59] * da[59] +
a88 * da[88] + a89 * da[89] + a90 * da[90] + a91 * da[91] +
a[92] * da[92] + a[93] * da[93];
// termes d'activite magnetique
da[40] = p10 * coste * dkp;
da[41] = p30 * coste * dkp;
da[42] = p50 * coste * dkp;
da[43] = p11 * ch * dkp;
da[44] = p31 * ch * dkp;
da[45] = p51 * ch * dkp;
da[46] = p11 * sh * dkp;
da[47] = p31 * sh * dkp;
da[48] = p51 * sh * dkp;
// fonction g[l] periodique supplementaire
fp += a[40] * da[40] + a[41] * da[41] + a[42] * da[42] + a[43] * da[43] +
a[44] * da[44] + a[45] * da[45] + a[46] * da[46] + a[47] * da[47] +
a[48] * da[48];
dakp = (a[40] * p10 + a[41] * p30 + a[42] * p50) * coste +
(a[43] * p11 + a[44] * p31 + a[45] * p51) * ch +
(a[46] * p11 + a[47] * p31 + a[48] * p51) * sh;
da[ikp] += dakp * akp[2];
da[ikp + 1] = da[ikp] + dakp * c2fi * akp[2];
// termes de longitude
da[49] = p11 * clfl;
da[50] = p21 * clfl;
da[51] = p31 * clfl;
da[52] = p41 * clfl;
da[53] = p51 * clfl;
da[54] = p11 * slfl;
da[55] = p21 * slfl;
da[56] = p31 * slfl;
da[57] = p41 * slfl;
da[58] = p51 * slfl;
// fonction g[l] periodique supplementaire
fp += a[49] * da[49] + a[50] * da[50] + a[51] * da[51] + a[52] * da[52] +
a[53] * da[53] + a[54] * da[54] + a[55] * da[55] + a[56] * da[56] +
a[57] * da[57] + a[58] * da[58];
// fonction g(l) totale (couplage avec le flux)
return f0 + fp * f1f;
}
/** Apply a correction coefficient to the given parameter.
* @param param the parameter to correct
* @return the corrected parameter
*/
private double semestrialCorrection(final double param) {
final int debeq_pr = 59;
final int debeq_au = 244;
final double result;
if (day >= 100) {
final double xmult = (day - debeq_au) / 40.0;
result = param - 2.0 * param * xmult;
} else {
final double xmult = (day - debeq_pr) / 40.0;
result = 2.0 * param * xmult - param;
}
return result;
}
}
/** Local holder for intermediate results ensuring the model is reentrant.
* @param <T> type of the field elements
*/
private static class FieldComputation<T extends RealFieldElement<T>> {
/** Number of days in current year. */
private int day;
/** Instant solar flux. f[1] = instantaneous flux; f[2] = 0. (not used). */
private double[] f = new double[3];
/** Mean solar flux. fbar[1] = mean flux; fbar[2] = 0. (not used). */
private double[] fbar = new double[3];
/** Kp coefficients.
* <ul>
* <li>akp[1] = 3-hourly kp</li>
* <li>akp[2] = 0 (not used)</li>
* <li>akp[3] = mean kp of last 24 hours</li>
* <li>akp[4] = 0 (not used)</li>
* </ul>
*/
private double[] akp = new double[5];
/** Cosine of the longitude. */
private final T clfl;
/** Sine of the longitude. */
private final T slfl;
/** Total density (g/cm3). */
private final T ro;
// CHECKSTYLE: stop JavadocVariable check
/** Legendre coefficients. */
private final T p10;
private final T p20;
private final T p30;
private final T p40;
private final T p50;
private final T p60;
private final T p11;
private final T p21;
private final T p31;
private final T p41;
private final T p51;
private final T p22;
private final T p32;
private final T p42;
private final T p52;
private final T p62;
private final T p33;
private final T p10mg;
private final T p20mg;
private final T p40mg;
/** Local time intermediate values. */
private final T hl0;
private final T ch;
private final T sh;
private final T c2h;
private final T s2h;
private final T c3h;
private final T s3h;
/** Simple constructor.
* @param day day of year
* @param altiKM altitude <em>in kilometers</em>
* @param lon local longitude (rad)
* @param lat local latitude (rad)
* @param hl local solar time in rad (O hr = 0 rad)
* @param f instantaneous solar flux (F10.7)
* @param far mean solar flux (F10.7)
* @param akp geomagnetic activity index
*/
FieldComputation(final int day,
final T altiKM, final T lon, final T lat,
final T hl, final double[] f, final double[] far,
final double[] akp) {
this.day = day;
this.f = f;
this.fbar = far;
this.akp = akp;
// compute Legendre polynomials wrt geographic pole
final T c = lat.sin();
final T c2 = c.multiply(c);
final T c4 = c2.multiply(c2);
final T s = lat.cos();
final T s2 = s.multiply(s);
p10 = c;
p20 = c2.multiply(1.5).subtract(0.5);
p30 = c.multiply(c2.multiply(2.5).subtract(1.5));
p40 = c4.multiply(4.375).subtract(c2.multiply(3.75)).add(0.375);
p50 = c.multiply(c4.multiply(7.875).subtract(c2.multiply(8.75)).add(1.875));
p60 = (c.multiply(5.5).multiply(p50).subtract(p40.multiply(2.5))).divide(3.0);
p11 = s;
p21 = c.multiply(3.0).multiply(s);
p31 = s.multiply(c2.multiply(7.5).subtract(1.5));
p41 = c.multiply(s).multiply(c2.multiply(17.5).subtract(7.5));
p51 = s.multiply(c4.multiply(39.375).subtract(c2.multiply(26.25)).add(1.875));
p22 = s2.multiply(3.0);
p32 = c.multiply(15.0).multiply(s2);
p42 = s2.multiply(c2.multiply(52.5).subtract(7.5));
p52 = c.multiply(3.0).multiply(p42).subtract(p32.multiply(2.0));
p62 = c.multiply(2.75).multiply(p52).subtract(p42.multiply(1.75));
p33 = s.multiply(15.0).multiply(s2);
// compute Legendre polynomials wrt magnetic pole (79N, 71W)
final T clmlmg = lon.subtract(XLMG).cos();
final T cmg = s.multiply(CPMG).multiply(clmlmg).add(c.multiply(SPMG));
final T cmg2 = cmg.multiply(cmg);
final T cmg4 = cmg2.multiply(cmg2);
p10mg = cmg;
p20mg = cmg2.multiply(1.5).subtract(0.5);
p40mg = cmg4.multiply(4.375).subtract(cmg2.multiply(3.75)).add(0.375);
clfl = lon.cos();
slfl = lon.sin();
// local time
hl0 = hl;
ch = hl0.cos();
sh = hl0.sin();
c2h = ch.multiply(ch).subtract(sh.multiply(sh));
s2h = ch.multiply(sh).multiply(2);
c3h = c2h.multiply(ch).subtract(s2h.multiply(sh));
s3h = s2h.multiply(ch).add(c2h.multiply(sh));
final double zlb = ZLB0; // + dzlb ??
final T[] dtt = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] dh = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] dhe = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] dox = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] daz2 = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] do2 = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] daz = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] dt0 = MathArrays.buildArray(altiKM.getField(), tt.length);
final T[] dtp = MathArrays.buildArray(altiKM.getField(), tt.length);
// compute function g(l) / tinf, t120, tp120
int kleq = 1;
final T gdelt = gFunction(tt, dtt, 1, kleq);
dtt[1] = gdelt.add(1);
final T tinf = dtt[1].multiply(tt[1]);
kleq = 0; // equinox
if ((day < 59) || (day > 284)) {
kleq = -1; // north winter
}
if ((day > 99) && (day < 244)) {
kleq = 1; // north summer
}
final T gdelt0 = gFunction(t0, dt0, 0, kleq);
dt0[1] = gdelt0.add(t0[1]).divide(t0[1]);
final T t120 = gdelt0.add(t0[1]);
final T gdeltp = gFunction(tp, dtp, 0, kleq);
dtp[1] = gdeltp.add(tp[1]).divide(tp[1]);
final T tp120 = gdeltp.add(tp[1]);
// compute n(z) concentrations: H, He, O, N2, O2, N
final T sigma = tp120.divide(tinf.subtract(t120));
final T dzeta = altiKM.add(RE).reciprocal().multiply(zlb + RE);
final T zeta = altiKM.subtract(zlb).multiply(dzeta);
final T sigzeta = sigma.multiply(zeta);
final T expsz = sigzeta.negate().exp();
final T tz = tinf.subtract(tinf.subtract(t120).multiply(expsz));
final T[] dbase = MathArrays.buildArray(altiKM.getField(), 7);
kleq = 1;
final T gdelh = gFunction(h, dh, 0, kleq);
dh[1] = gdelh.exp();
dbase[1] = dh[1].multiply(h[1]);
final T gdelhe = gFunction(he, dhe, 0, kleq);
dhe[1] = gdelhe.exp();
dbase[2] = dhe[1].multiply(he[1]);
final T gdelo = gFunction(o, dox, 1, kleq);
dox[1] = gdelo.exp();
dbase[3] = dox[1].multiply(o[1]);
final T gdelaz2 = gFunction(az2, daz2, 1, kleq);
daz2[1] = gdelaz2.exp();
dbase[4] = daz2[1].multiply(az2[1]);
final T gdelo2 = gFunction(o2, do2, 1, kleq);
do2[1] = gdelo2.exp();
dbase[5] = do2[1].multiply(o2[1]);
final T gdelaz = gFunction(az, daz, 1, kleq);
daz[1] = gdelaz.exp();
dbase[6] = daz[1].multiply(az[1]);
final double zlbre = 1.0 + zlb / RE;
final T glb = sigma.multiply(RGAS).multiply(tinf).reciprocal().multiply(GSURF / (zlbre * zlbre));
final T t120tz = t120.divide(tz);
// Partial densities in (g/cm3).
// d(1) = hydrogen
// d(2) = helium
// d(3) = atomic oxygen
// d(4) = molecular nitrogen
// d(5) = molecular oxygen
// d(6) = atomic nitrogen
T tmpro = altiKM.getField().getZero();
for (int i = 1; i <= 6; i++) {
final T gamma = glb.multiply(MA[i]);
final T upapg = gamma.add(1.0 + ALEFA[i]);
final T fzI = t120tz.pow(upapg).multiply(sigzeta.negate().multiply(gamma).exp());
// concentrations of H, He, O, N2, O2, N (particles/cm³)
final T ccI = dbase[i].multiply(fzI);
// contribution of densities of H, He, O, N2, O2, N (g/cm³)
tmpro = tmpro.add(ccI.multiply(VMA[i]));
}
this.ro = tmpro;
}
/** Computation of function G.
* @param a vector of coefficients for computation
* @param da vector of partial derivatives
* @param ff0 coefficient flag (1 for Ox, Az, He, T°; 0 for H and tp120)
* @param kle_eq season indicator flag (summer, winter, equinox)
* @return value of G
*/
private T gFunction(final double[] a, final T[] da,
final int ff0, final int kle_eq) {
final T zero = da[0].getField().getZero();
final double[] fmfb = new double[3];
final double[] fbm150 = new double[3];
// latitude terms
da[2] = p20;
da[3] = p40;
da[74] = p10;
double a74 = a[74];
double a77 = a[77];
double a78 = a[78];
if (kle_eq == -1) {
// winter
a74 = -a74;
a77 = -a77;
a78 = -a78;
}
if (kle_eq == 0 ) {
// equinox
a74 = semestrialCorrection(a74);
a77 = semestrialCorrection(a77);
a78 = semestrialCorrection(a78);
}
da[77] = p30;
da[78] = p50;
da[79] = p60;
// flux terms
fmfb[1] = f[1] - fbar[1];
fmfb[2] = f[2] - fbar[2];
fbm150[1] = fbar[1] - 150.0;
fbm150[2] = fbar[2];
da[4] = zero.add(fmfb[1]);
da[6] = zero.add(fbm150[1]);
da[4] = da[4].add(a[70] * fmfb[2]);
da[6] = da[6].add(a[71] * fbm150[2]);
da[70] = da[4].multiply(a[ 5]).multiply(2).
add(p10.multiply(a[82])).
add(p20.multiply(a[83])).
add(p30.multiply(a[84])).
add(a[4]).
multiply(fmfb[2]);
da[71] = da[6].multiply(a[69]).multiply(2).
add(p10.multiply(a[85])).
add(p20.multiply(a[86])).
add(p30.multiply(a[87])).
add(a[6]).
multiply(fbm150[2]);
da[5] = da[4].multiply(da[4]);
da[69] = da[6].multiply(da[6]);
da[82] = da[4].multiply(p10);
da[83] = da[4].multiply(p20);
da[84] = da[4].multiply(p30);
da[85] = da[6].multiply(p20);
da[86] = da[6].multiply(p30);
da[87] = da[6].multiply(p40);
// Kp terms
final int ikp = 62;
final int ikpm = 67;
final T c2fi = p10mg.multiply(p10mg).negate().add(1);
final T dkp = c2fi.multiply(a[ikp + 1]).add(a[ikp]).multiply(akp[2]).add(akp[1]);
T dakp = p20mg.multiply(a[8]).add(p40mg.multiply(a[68])).
add(p20mg.multiply(a[61]).add(dkp.multiply(dkp).multiply(2 * a[75]).add(a[60])).multiply(dkp.multiply(2))).
add(a[7]);
da[ikp] = dakp.multiply(akp[2]);
da[ikp + 1] = da[ikp].multiply(c2fi);
final double dkpm = akp[3] + a[ikpm] * akp[4];
final T dakpm = p20mg.multiply(a[65]).add(p40mg.multiply(a[72])).
add(p20mg.multiply(a[73]).add(a[66] + a[76] * 2.0 * dkpm * dkpm).multiply( 2.0 * dkpm)).
add(a[64]);
da[ikpm] = dakpm.multiply(akp[4]);
da[7] = dkp;
da[8] = p20mg.multiply(dkp);
da[68] = p40mg.multiply(dkp);
da[60] = dkp.multiply(dkp);
da[61] = p20mg.multiply(da[60]);
da[75] = da[60].multiply(da[60]);
da[64] = zero.add(dkpm);
da[65] = p20mg.multiply(dkpm);
da[72] = p40mg.multiply(dkpm);
da[66] = zero.add(dkpm * dkpm);
da[73] = p20mg.multiply(da[66]);
da[76] = da[66].multiply(da[66]);
// non-periodic g(l) function
T f0 = da[4].multiply(a[4]).
add(da[5].multiply(a[5])).
add(da[6].multiply(a[6])).
add(da[69].multiply(a[69])).
add(da[82].multiply(a[82])).
add(da[83].multiply(a[83])).
add(da[84].multiply(a[84])).
add(da[85].multiply(a[85])).
add(da[86].multiply(a[86])).
add(da[87].multiply(a[87]));
final T f1f = f0.multiply(ff0).add(1);
f0 = f0.
add(da[2].multiply(a[2])).
add(da[3].multiply(a[3])).
add(da[7].multiply(a[7])).
add(da[8].multiply(a[8])).
add(da[60].multiply(a[60])).
add(da[61].multiply(a[61])).
add(da[68].multiply(a[68])).
add(da[64].multiply(a[64])).
add(da[65].multiply(a[65])).
add(da[66].multiply(a[66])).
add(da[72].multiply(a[72])).
add(da[73].multiply(a[73])).
add(da[74].multiply(a74)).
add(da[75].multiply(a[75])).
add(da[76].multiply(a[76])).
add(da[77].multiply(a77)).
add(da[78].multiply(a78)).
add(da[79].multiply(a[79]));
// termes annuels symetriques en latitude
da[9] = zero.add(FastMath.cos(ROT * (day - a[11])));
da[10] = p20.multiply(da[9]);
// termes semi-annuels symetriques en latitude
da[12] = zero.add(FastMath.cos(ROT2 * (day - a[14])));
da[13] = p20.multiply(da[12]);
// termes annuels non symetriques en latitude
final double coste = FastMath.cos(ROT * (day - a[18]));
da[15] = p10.multiply(coste);
da[16] = p30.multiply(coste);
da[17] = p50.multiply(coste);
// terme semi-annuel non symetrique en latitude
final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
da[19] = p10.multiply(cos2te);
da[39] = p30.multiply(cos2te);
da[59] = p50.multiply(cos2te);
// termes diurnes [et couples annuel]
da[21] = p11.multiply(ch);
da[22] = p31.multiply(ch);
da[23] = p51.multiply(ch);
da[24] = da[21].multiply(coste);
da[25] = p21.multiply(ch).multiply(coste);
da[26] = p11.multiply(sh);
da[27] = p31.multiply(sh);
da[28] = p51.multiply(sh);
da[29] = da[26].multiply(coste);
da[30] = p21.multiply(sh).multiply(coste);
// termes semi-diurnes [et couples annuel]
da[31] = p22.multiply(c2h);
da[37] = p42.multiply(c2h);
da[32] = p32.multiply(c2h).multiply(coste);
da[33] = p22.multiply(s2h);
da[38] = p42.multiply(s2h);
da[34] = p32.multiply(s2h).multiply(coste);
da[88] = p32.multiply(c2h);
da[89] = p32.multiply(s2h);
da[90] = p52.multiply(c2h);
da[91] = p52.multiply(s2h);
double a88 = a[88];
double a89 = a[89];
double a90 = a[90];
double a91 = a[91];
if (kle_eq == -1) { //hiver
a88 = -a88;
a89 = -a89;
a90 = -a90;
a91 = -a91;
}
if (kle_eq == 0) { //equinox
a88 = semestrialCorrection(a88);
a89 = semestrialCorrection(a89);
a90 = semestrialCorrection(a90);
a91 = semestrialCorrection(a91);
}
da[92] = p62.multiply(c2h);
da[93] = p62.multiply(s2h);
// termes ter-diurnes
da[35] = p33.multiply(c3h);
da[36] = p33.multiply(s3h);
// fonction g[l] periodique
T fp = da[ 9].multiply(a[ 9]) .add(da[10].multiply(a[10])).add(da[12].multiply(a[12])).add(da[13].multiply(a[13])).
add(da[15].multiply(a[15])).add(da[16].multiply(a[16])).add(da[17].multiply(a[17])).add(da[19].multiply(a[19])).
add(da[21].multiply(a[21])).add(da[22].multiply(a[22])).add(da[23].multiply(a[23])).add(da[24].multiply(a[24])).
add(da[25].multiply(a[25])).add(da[26].multiply(a[26])).add(da[27].multiply(a[27])).add(da[28].multiply(a[28])).
add(da[29].multiply(a[29])).add(da[30].multiply(a[30])).add(da[31].multiply(a[31])).add(da[32].multiply(a[32])).
add(da[33].multiply(a[33])).add(da[34].multiply(a[34])).add(da[35].multiply(a[35])).add(da[36].multiply(a[36])).
add(da[37].multiply(a[37])).add(da[38].multiply(a[38])).add(da[39].multiply(a[39])).add(da[59].multiply(a[59])).
add(da[88].multiply(a88)) .add(da[89].multiply(a89) ).add(da[90].multiply(a90) ).add(da[91].multiply(a91) ).
add(da[92].multiply(a[92])).add(da[93].multiply(a[93]));
// termes d'activite magnetique
da[40] = p10.multiply(coste).multiply(dkp);
da[41] = p30.multiply(coste).multiply(dkp);
da[42] = p50.multiply(coste).multiply(dkp);
da[43] = p11.multiply(ch).multiply(dkp);
da[44] = p31.multiply(ch).multiply(dkp);
da[45] = p51.multiply(ch).multiply(dkp);
da[46] = p11.multiply(sh).multiply(dkp);
da[47] = p31.multiply(sh).multiply(dkp);
da[48] = p51.multiply(sh).multiply(dkp);
// fonction g[l] periodique supplementaire
fp = fp.
add(da[40].multiply(a[40])).
add(da[41].multiply(a[41])).
add(da[42].multiply(a[42])).
add(da[43].multiply(a[43])).
add(da[44].multiply(a[44])).
add(da[45].multiply(a[45])).
add(da[46].multiply(a[46])).
add(da[47].multiply(a[47])).
add(da[48].multiply(a[48]));
dakp = p10.multiply(a[40]).add(p30.multiply(a[41])).add(p50.multiply(a[42])).multiply(coste).
add(p11.multiply(a[40]).add(p31.multiply(a[44])).add(p51.multiply(a[45])).multiply(ch)).
add(p11.multiply(a[46]).add(p31.multiply(a[47])).add(p51.multiply(a[48])).multiply(sh));
da[ikp] = da[ikp].add(dakp.multiply(akp[2]));
da[ikp + 1] = da[ikp].add(dakp.multiply(c2fi).multiply(akp[2]));
// termes de longitude
da[49] = p11.multiply(clfl);
da[50] = p21.multiply(clfl);
da[51] = p31.multiply(clfl);
da[52] = p41.multiply(clfl);
da[53] = p51.multiply(clfl);
da[54] = p11.multiply(slfl);
da[55] = p21.multiply(slfl);
da[56] = p31.multiply(slfl);
da[57] = p41.multiply(slfl);
da[58] = p51.multiply(slfl);
// fonction g[l] periodique supplementaire
fp = fp.
add(da[49].multiply(a[49])).
add(da[50].multiply(a[50])).
add(da[51].multiply(a[51])).
add(da[52].multiply(a[52])).
add(da[53].multiply(a[53])).
add(da[54].multiply(a[54])).
add(da[55].multiply(a[55])).
add(da[56].multiply(a[56])).
add(da[57].multiply(a[57])).
add(da[58].multiply(a[58]));
// fonction g(l) totale (couplage avec le flux)
return f0.add(fp.multiply(f1f));
}
/** Apply a correction coefficient to the given parameter.
* @param param the parameter to correct
* @return the corrected parameter
*/
private double semestrialCorrection(final double param) {
final int debeq_pr = 59;
final int debeq_au = 244;
final double result;
if (day >= 100) {
final double xmult = (day - debeq_au) / 40.0;
result = param - 2.0 * param * xmult;
} else {
final double xmult = (day - debeq_pr) / 40.0;
result = 2.0 * param * xmult - param;
}
return result;
}
}
}