DTM2000.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 org.hipparchus.CalculusFieldElement;
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.FieldSinCos;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.SinCos;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
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.TimeScale;
import org.orekit.utils.ExtendedPositionProvider;

import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.nio.charset.StandardCharsets;
import java.util.Arrays;

/** 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 index, which goes from 1 (very low activity) to
 * 9 (high activity).
 *
 * <p>
 * All these data can be found on the <a href="https://www.noaa.gov/">
 * 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 extends AbstractSunInfluencedAtmosphere {
    // 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;

    /** External data container. */
    private DTM2000InputParameters inputParams;

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

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

    /** Simple constructor for independent computation.
     *
     * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
     *
     * @param parameters the solar and magnetic activity data
     * @param sun the sun position
     * @param earth the earth body shape
     * @see #DTM2000(DTM2000InputParameters, ExtendedPositionProvider, BodyShape, TimeScale)
     */
    @DefaultDataContext
    public DTM2000(final DTM2000InputParameters parameters,
                   final ExtendedPositionProvider sun, final BodyShape earth) {
        this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
    }

    /** Simple constructor for independent computation.
     * @param parameters the solar and magnetic activity data
     * @param sun the sun position
     * @param earth the earth body shape
     * @param utc UTC time scale.
     * @since 10.1
     */
    public DTM2000(final DTM2000InputParameters parameters,
                   final ExtendedPositionProvider sun,
                   final BodyShape earth,
                   final TimeScale utc) {
        super(sun);

        synchronized (DTM2000.class) {
            // lazy reading of model coefficients
            if (tt == null) {
                readcoefficients();
            }
        }

        this.earth = earth;
        this.inputParams = parameters;

        this.utc = utc;
    }

    /** {@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³)
     */
    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) {
        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³)
          * @since 9.0
     */
    public <T extends CalculusFieldElement<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) {
        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.
     */
    private static void readcoefficients() {

        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];

        try (InputStream in = checkNull(DTM2000.class.getResourceAsStream(DTM2000));
             BufferedReader r = new BufferedReader(new InputStreamReader(in, StandardCharsets.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()));
        }
    }

    /** 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³)
     */
    public double getDensity(final AbsoluteDate date, final Vector3D position,
                             final Frame frame) {

        // 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 int day = date.getComponents(utc).getDate().getDayOfYear();
        //position in ECEF so we only have to do the transform once
        final Frame ecef = earth.getBodyFrame();
        final Vector3D pEcef = frame.getStaticTransformTo(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 = getSunPosition(date, ecef);
        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 CalculusFieldElement<T>> T
        getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position,
                   final Frame frame) {
        // 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 int day = date.getComponents(utc).getDate().getDayOfYear();
        // position in ECEF so we only have to do the transform once
        final Frame ecef = earth.getBodyFrame();
        final FieldVector3D<T> pEcef = frame.getStaticTransformTo(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 FieldVector3D<T> sunPos = getSunPosition(date, ecef);
        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(y.getPi());

        // 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));
    }

    /**
     * Helper method to check for null resources. Throws an exception if {@code
     * stream} is null.
     *
     * @param stream loaded from the class resources.
     * @return {@code stream}.
     */
    private static InputStream checkNull(final InputStream stream) {
        if (stream == null) {
            throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_RESOURCE, DTM2000);
        }
        return stream;
    }

    /** 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;

            // Sine and cosine of local latitude and longitude
            final SinCos scLat = FastMath.sinCos(lat);
            final SinCos scLon = FastMath.sinCos(lon);

            // compute Legendre polynomials wrt geographic pole
            final double c = scLat.sin();
            final double c2 = c * c;
            final double c4 = c2 * c2;
            final double s = scLat.cos();
            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 = scLon.cos();
            slfl = scLon.sin();

            // local time
            hl0 = hl;
            final SinCos scHlo = FastMath.sinCos(hl0);
            ch  = scHlo.cos();
            sh  = scHlo.sin();
            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.exp(FastMath.log(t120tz) * upapg - 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 CalculusFieldElement<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;

            // Sine and cosine of local latitude and longitude
            final FieldSinCos<T> scLat = FastMath.sinCos(lat);
            final FieldSinCos<T> scLon = FastMath.sinCos(lon);

            // compute Legendre polynomials wrt geographic pole
            final T c = scLat.sin();
            final T c2 = c.square();
            final T c4 = c2.square();
            final T s = scLat.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 = scLon.cos();
            slfl = scLon.sin();

            // local time
            hl0 = hl;
            final FieldSinCos<T> scHlo = FastMath.sinCos(hl0);
            ch  = scHlo.cos();
            sh  = scHlo.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.log().multiply(upapg).subtract(sigzeta.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.newInstance(fmfb[1]);
            da[6]     = zero.newInstance(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.square().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.square();
            da[61]   = p20mg.multiply(da[60]);
            da[75]   = da[60].square();
            da[64]   = zero.newInstance(dkpm);
            da[65]   = p20mg.multiply(dkpm);
            da[72]   = p40mg.multiply(dkpm);
            da[66]   = zero.newInstance(dkpm * dkpm);
            da[73]   = p20mg.multiply(da[66]);
            da[76]   = da[66].square();

            // 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.newInstance(FastMath.cos(ROT * (day - a[11])));
            da[10] = p20.multiply(da[9]);
//          termes semi-annuels symetriques en latitude
            da[12] = zero.newInstance(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;
        }

    }

}