DTM2000.java

  1. /* Copyright 2002-2022 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.models.earth.atmosphere;

  18. import java.io.BufferedReader;
  19. import java.io.IOException;
  20. import java.io.InputStream;
  21. import java.io.InputStreamReader;
  22. import java.nio.charset.StandardCharsets;
  23. import java.util.Arrays;

  24. import org.hipparchus.CalculusFieldElement;
  25. import org.hipparchus.exception.DummyLocalizable;
  26. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  27. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.FieldSinCos;
  30. import org.hipparchus.util.MathArrays;
  31. import org.hipparchus.util.SinCos;
  32. import org.orekit.annotation.DefaultDataContext;
  33. import org.orekit.bodies.BodyShape;
  34. import org.orekit.bodies.FieldGeodeticPoint;
  35. import org.orekit.bodies.GeodeticPoint;
  36. import org.orekit.data.DataContext;
  37. import org.orekit.errors.OrekitException;
  38. import org.orekit.errors.OrekitMessages;
  39. import org.orekit.frames.Frame;
  40. import org.orekit.time.AbsoluteDate;
  41. import org.orekit.time.FieldAbsoluteDate;
  42. import org.orekit.time.TimeScale;
  43. import org.orekit.utils.PVCoordinatesProvider;

  44. /** This atmosphere model is the realization of the DTM-2000 model.
  45.  * <p>
  46.  * It is described in the paper: <br>
  47.  *
  48.  * <b>The DTM-2000 empirical thermosphere model with new data assimilation
  49.  *  and constraints at lower boundary: accuracy and properties</b><br>
  50.  *
  51.  * <i>S. Bruinsma, G. Thuillier and F. Barlier</i> <br>
  52.  *
  53.  * Journal of Atmospheric and Solar-Terrestrial Physics 65 (2003) 1053–1070<br>
  54.  *
  55.  *</p>
  56.  *<p>
  57.  * This model provides dense output for altitudes beyond 120 km.
  58.  *</p>
  59.  *
  60.  * <p>
  61.  * The model needs geographical and time information to compute general values,
  62.  * but also needs space weather data : mean and instantaneous solar flux and
  63.  * geomagnetic indices.
  64.  * </p>
  65.  * <p>
  66.  * Mean solar flux is (for the moment) represented by the F10.7 indices. Instantaneous
  67.  * flux can be set to the mean value if the data is not available. Geomagnetic
  68.  * activity is represented by the Kp indice, which goes from 1 (very low activity) to
  69.  * 9 (high activity).
  70.  *
  71.  * <p>
  72.  * All these data can be found on the <a href="https://www.noaa.gov/">
  73.  * NOAA (National Oceanic and Atmospheric Administration) website.</a>
  74.  * </p>
  75.  *
  76.  *
  77.  * @author R. Biancale, S. Bruinsma: original fortran routine
  78.  * @author Fabien Maussion (java translation)
  79.  */
  80. public class DTM2000 implements Atmosphere {

  81.     /** Identifier for hydrogen. */
  82.     public static final int HYDROGEN = 1;

  83.     /** Identifier for helium. */
  84.     public static final int HELIUM = 2;

  85.     /** Identifier for atomic oxygen. */
  86.     public static final int ATOMIC_OXYGEN = 3;

  87.     /** Identifier for molecular nitrogen. */
  88.     public static final int MOLECULAR_NITROGEN = 4;

  89.     /** Identifier for molecular oxygen. */
  90.     public static final int MOLECULAR_OXYGEN = 5;

  91.     /** Identifier for atomic nitrogen. */
  92.     public static final int ATOMIC_NITROGEN = 6;

  93.     /** Serializable UID. */
  94.     private static final long serialVersionUID = 20170705L;

  95.     // Constants :

  96.     /** Number of parameters. */
  97.     private static final int NLATM = 96;

  98.     /** Thermal diffusion coefficient. */
  99.     private static final double[] ALEFA = {
  100.         0, -0.40, -0.38, 0, 0, 0, 0
  101.     };

  102.     /** Atomic mass  H, He, O, N2, O2, N. */
  103.     private static final double[] MA = {
  104.         0, 1, 4, 16, 28, 32, 14
  105.     };

  106.     /** Atomic mass  H, He, O, N2, O2, N. */
  107.     private static final double[] VMA = {
  108.         0, 1.6606e-24, 6.6423e-24, 26.569e-24, 46.4958e-24, 53.1381e-24, 23.2479e-24
  109.     };

  110.     /** Polar Earth radius. */
  111.     private static final double RE = 6356.77;

  112.     /** Reference altitude. */
  113.     private static final double ZLB0 = 120.0;

  114.     /** Cosine of the latitude of the magnetic pole (79N, 71W). */
  115.     private static final double CPMG = .19081;

  116.     /** Sine of the latitude of the magnetic pole (79N, 71W). */
  117.     private static final double SPMG = .98163;

  118.     /** Longitude (in radians) of the magnetic pole (79N, 71W). */
  119.     private static final double XLMG = -1.2392;

  120.     /** Gravity acceleration at 120 km altitude. */
  121.     private static final double GSURF = 980.665;

  122.     /** Universal gas constant. */
  123.     private static final double RGAS = 831.4;

  124.     /** 2 * π / 365. */
  125.     private static final double ROT = 0.017214206;

  126.     /** 2 * rot. */
  127.     private static final double ROT2 = 0.034428412;

  128.     /** Resources text file. */
  129.     private static final String DTM2000 = "/assets/org/orekit/dtm_2000.txt";

  130.     // CHECKSTYLE: stop JavadocVariable check

  131.     /** Elements coefficients. */
  132.     private static double[] tt   = null;
  133.     private static double[] h    = null;
  134.     private static double[] he   = null;
  135.     private static double[] o    = null;
  136.     private static double[] az2  = null;
  137.     private static double[] o2   = null;
  138.     private static double[] az   = null;
  139.     private static double[] t0   = null;
  140.     private static double[] tp   = null;

  141.     /** Sun position. */
  142.     private PVCoordinatesProvider sun;

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

  145.     /** Earth body shape. */
  146.     private BodyShape earth;

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

  149.     /** Simple constructor for independent computation.
  150.      *
  151.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  152.      *
  153.      * @param parameters the solar and magnetic activity data
  154.      * @param sun the sun position
  155.      * @param earth the earth body shape
  156.      * @see #DTM2000(DTM2000InputParameters, PVCoordinatesProvider, BodyShape, TimeScale)
  157.      */
  158.     @DefaultDataContext
  159.     public DTM2000(final DTM2000InputParameters parameters,
  160.                    final PVCoordinatesProvider sun, final BodyShape earth) {
  161.         this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
  162.     }

  163.     /** Simple constructor for independent computation.
  164.      * @param parameters the solar and magnetic activity data
  165.      * @param sun the sun position
  166.      * @param earth the earth body shape
  167.      * @param utc UTC time scale.
  168.      * @since 10.1
  169.      */
  170.     public DTM2000(final DTM2000InputParameters parameters,
  171.                    final PVCoordinatesProvider sun,
  172.                    final BodyShape earth,
  173.                    final TimeScale utc) {

  174.         synchronized (DTM2000.class) {
  175.             // lazy reading of model coefficients
  176.             if (tt == null) {
  177.                 readcoefficients();
  178.             }
  179.         }

  180.         this.earth = earth;
  181.         this.sun = sun;
  182.         this.inputParams = parameters;

  183.         this.utc = utc;
  184.     }

  185.     /** {@inheritDoc} */
  186.     public Frame getFrame() {
  187.         return earth.getBodyFrame();
  188.     }

  189.     /** Get the local density with initial entries.
  190.      * @param day day of year
  191.      * @param alti altitude in meters
  192.      * @param lon local longitude (rad)
  193.      * @param lat local latitude (rad)
  194.      * @param hl local solar time in rad (O hr = 0 rad)
  195.      * @param f instantaneous solar flux (F10.7)
  196.      * @param fbar mean solar flux (F10.7)
  197.      * @param akp3 3 hrs geomagnetic activity index (1-9)
  198.      * @param akp24 Mean of last 24 hrs geomagnetic activity index (1-9)
  199.      * @return the local density (kg/m³)
  200.      */
  201.     public double getDensity(final int day,
  202.                              final double alti, final double lon, final double lat,
  203.                              final double hl, final double f, final double fbar,
  204.                              final double akp3, final double akp24) {
  205.         final double threshold = 120000;
  206.         if (alti < threshold) {
  207.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
  208.                                       alti, threshold);
  209.         }
  210.         final Computation result = new Computation(day, alti / 1000, lon, lat, hl,
  211.                                                    new double[] {
  212.                                                        0, f, 0
  213.                                                    }, new double[] {
  214.                                                        0, fbar, 0
  215.                                                    }, new double[] {
  216.                                                        0, akp3, 0, akp24, 0
  217.                                                    });
  218.         return result.ro * 1000;
  219.     }

  220.     /** Get the local density with initial entries.
  221.      * @param day day of year
  222.      * @param alti altitude in meters
  223.      * @param lon local longitude (rad)
  224.      * @param lat local latitude (rad)
  225.      * @param hl local solar time in rad (O hr = 0 rad)
  226.      * @param f instantaneous solar flux (F10.7)
  227.      * @param fbar mean solar flux (F10.7)
  228.      * @param akp3 3 hrs geomagnetic activity index (1-9)
  229.      * @param akp24 Mean of last 24 hrs geomagnetic activity index (1-9)
  230.      * @param <T> type of the field elements
  231.      * @return the local density (kg/m³)
  232.           * @since 9.0
  233.      */
  234.     public <T extends CalculusFieldElement<T>> T getDensity(final int day,
  235.                                                         final T alti, final T lon, final T lat,
  236.                                                         final T hl, final double f, final double fbar,
  237.                                                         final double akp3, final double akp24) {
  238.         final double threshold = 120000;
  239.         if (alti.getReal() < threshold) {
  240.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
  241.                                       alti, threshold);
  242.         }
  243.         final FieldComputation<T> result = new FieldComputation<>(day, alti.divide(1000), lon, lat, hl,
  244.                                                                   new double[] {
  245.                                                                       0, f, 0
  246.                                                                   }, new double[] {
  247.                                                                       0, fbar, 0
  248.                                                                   }, new double[] {
  249.                                                                       0, akp3, 0, akp24, 0
  250.                                                                   });
  251.         return result.ro.multiply(1000);
  252.     }

  253.     /** Store the DTM model elements coefficients in internal arrays.
  254.      */
  255.     private static void readcoefficients() {

  256.         final int size = NLATM + 1;
  257.         tt   = new double[size];
  258.         h    = new double[size];
  259.         he   = new double[size];
  260.         o    = new double[size];
  261.         az2  = new double[size];
  262.         o2   = new double[size];
  263.         az   = new double[size];
  264.         t0   = new double[size];
  265.         tp   = new double[size];

  266.         try (InputStream in = checkNull(DTM2000.class.getResourceAsStream(DTM2000));
  267.              BufferedReader r = new BufferedReader(new InputStreamReader(in, StandardCharsets.UTF_8))) {

  268.             r.readLine();
  269.             r.readLine();
  270.             for (String line = r.readLine(); line != null; line = r.readLine()) {
  271.                 final int num = Integer.parseInt(line.substring(0, 4).replace(' ', '0'));
  272.                 line = line.substring(4);
  273.                 tt[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  274.                 line = line.substring(13 + 9);
  275.                 h[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  276.                 line = line.substring(13 + 9);
  277.                 he[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  278.                 line = line.substring(13 + 9);
  279.                 o[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  280.                 line = line.substring(13 + 9);
  281.                 az2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  282.                 line = line.substring(13 + 9);
  283.                 o2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  284.                 line = line.substring(13 + 9);
  285.                 az[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  286.                 line = line.substring(13 + 9);
  287.                 t0[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  288.                 line = line.substring(13 + 9);
  289.                 tp[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  290.             }
  291.         } catch (IOException ioe) {
  292.             throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
  293.         }
  294.     }

  295.     /** Get the local density.
  296.      * @param date current date
  297.      * @param position current position in frame
  298.      * @param frame the frame in which is defined the position
  299.      * @return local density (kg/m³)
  300.      */
  301.     public double getDensity(final AbsoluteDate date, final Vector3D position,
  302.                              final Frame frame) {

  303.         // check if data are available :
  304.         if (date.compareTo(inputParams.getMaxDate()) > 0 ||
  305.             date.compareTo(inputParams.getMinDate()) < 0) {
  306.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  307.                                       date, inputParams.getMinDate(), inputParams.getMaxDate());
  308.         }

  309.         // compute day number in current year
  310.         final int day = date.getComponents(utc).getDate().getDayOfYear();
  311.         //position in ECEF so we only have to do the transform once
  312.         final Frame ecef = earth.getBodyFrame();
  313.         final Vector3D pEcef = frame.getStaticTransformTo(ecef, date)
  314.                 .transformPosition(position);
  315.         // compute geodetic position
  316.         final GeodeticPoint inBody = earth.transform(pEcef, ecef, date);
  317.         final double alti = inBody.getAltitude();
  318.         final double lon = inBody.getLongitude();
  319.         final double lat = inBody.getLatitude();

  320.         // compute local solar time
  321.         final Vector3D sunPos = sun.getPVCoordinates(date, ecef).getPosition();
  322.         final double hl = FastMath.PI + FastMath.atan2(
  323.                 sunPos.getX() * pEcef.getY() - sunPos.getY() * pEcef.getX(),
  324.                 sunPos.getX() * pEcef.getX() + sunPos.getY() * pEcef.getY());

  325.         // get current solar activity data and compute
  326.         return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(date),
  327.                           inputParams.getMeanFlux(date), inputParams.getThreeHourlyKP(date),
  328.                           inputParams.get24HoursKp(date));

  329.     }

  330.     /** {@inheritDoc} */
  331.     @Override
  332.     public <T extends CalculusFieldElement<T>> T
  333.         getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position,
  334.                    final Frame frame) {
  335.         // check if data are available :
  336.         final AbsoluteDate dateD = date.toAbsoluteDate();
  337.         if (dateD.compareTo(inputParams.getMaxDate()) > 0 ||
  338.             dateD.compareTo(inputParams.getMinDate()) < 0) {
  339.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  340.                                       dateD, inputParams.getMinDate(), inputParams.getMaxDate());
  341.         }

  342.         // compute day number in current year
  343.         final int day = date.getComponents(utc).getDate().getDayOfYear();
  344.         // position in ECEF so we only have to do the transform once
  345.         final Frame ecef = earth.getBodyFrame();
  346.         final FieldVector3D<T> pEcef = frame.getTransformTo(ecef, date).transformPosition(position);
  347.         // compute geodetic position
  348.         final FieldGeodeticPoint<T> inBody = earth.transform(pEcef, ecef, date);
  349.         final T alti = inBody.getAltitude();
  350.         final T lon = inBody.getLongitude();
  351.         final T lat = inBody.getLatitude();

  352.         // compute local solar time
  353.         final Vector3D sunPos = sun.getPVCoordinates(dateD, ecef).getPosition();
  354.         final T y  = pEcef.getY().multiply(sunPos.getX()).subtract(pEcef.getX().multiply(sunPos.getY()));
  355.         final T x  = pEcef.getX().multiply(sunPos.getX()).add(pEcef.getY().multiply(sunPos.getY()));
  356.         final T hl = y.atan2(x).add(y.getPi());

  357.         // get current solar activity data and compute
  358.         return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(dateD),
  359.                           inputParams.getMeanFlux(dateD), inputParams.getThreeHourlyKP(dateD),
  360.                           inputParams.get24HoursKp(dateD));
  361.     }

  362.     /**
  363.      * Helper method to check for null resources. Throws an exception if {@code
  364.      * stream} is null.
  365.      *
  366.      * @param stream loaded from the class resources.
  367.      * @return {@code stream}.
  368.      */
  369.     private static InputStream checkNull(final InputStream stream) {
  370.         if (stream == null) {
  371.             throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_RESOURCE, DTM2000);
  372.         }
  373.         return stream;
  374.     }

  375.     /** Local holder for intermediate results ensuring the model is reentrant. */
  376.     private static class Computation {

  377.         /** Number of days in current year. */
  378.         private final int day;

  379.         /** Instant solar flux. f[1] = instantaneous flux; f[2] = 0. (not used). */
  380.         private final double[] f;

  381.         /** Mean solar flux. fbar[1] = mean flux; fbar[2] = 0. (not used). */
  382.         private final double[] fbar;

  383.         /** Kp coefficients.
  384.          * <ul>
  385.          *   <li>akp[1] = 3-hourly kp</li>
  386.          *   <li>akp[2] = 0 (not used)</li>
  387.          *   <li>akp[3] = mean kp of last 24 hours</li>
  388.          *   <li>akp[4] = 0 (not used)</li>
  389.          * </ul>
  390.          */
  391.         private final double[] akp;

  392.         /** Cosine of the longitude. */
  393.         private final double clfl;

  394.         /** Sine of the longitude. */
  395.         private final double slfl;

  396.         /** Total density (g/cm3). */
  397.         private final double ro;

  398.         // CHECKSTYLE: stop JavadocVariable check

  399.         /** Legendre coefficients. */
  400.         private final double p10;
  401.         private final double p20;
  402.         private final double p30;
  403.         private final double p40;
  404.         private final double p50;
  405.         private final double p60;
  406.         private final double p11;
  407.         private final double p21;
  408.         private final double p31;
  409.         private final double p41;
  410.         private final double p51;
  411.         private final double p22;
  412.         private final double p32;
  413.         private final double p42;
  414.         private final double p52;
  415.         private final double p62;
  416.         private final double p33;
  417.         private final double p10mg;
  418.         private final double p20mg;
  419.         private final double p40mg;

  420.         /** Local time intermediate values. */
  421.         private final double hl0;
  422.         private final double ch;
  423.         private final double sh;
  424.         private final double c2h;
  425.         private final double s2h;
  426.         private final double c3h;
  427.         private final double s3h;

  428.         /** Simple constructor.
  429.          * @param day day of year
  430.          * @param altiKM altitude <em>in kilometers</em>
  431.          * @param lon local longitude (rad)
  432.          * @param lat local latitude (rad)
  433.          * @param hl local solar time in rad (O hr = 0 rad)
  434.          * @param f instantaneous solar flux (F10.7)
  435.          * @param fbar mean solar flux (F10.7)
  436.          * @param akp geomagnetic activity index
  437.          */
  438.         Computation(final int day,
  439.                     final double altiKM, final double lon, final double lat,
  440.                     final double hl, final double[] f, final double[] fbar,
  441.                     final double[] akp) {

  442.             this.day  = day;
  443.             this.f    = f;
  444.             this.fbar = fbar;
  445.             this.akp  = akp;

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

  449.             // compute Legendre polynomials wrt geographic pole
  450.             final double c = scLat.sin();
  451.             final double c2 = c * c;
  452.             final double c4 = c2 * c2;
  453.             final double s = scLat.cos();
  454.             final double s2 = s * s;
  455.             p10 = c;
  456.             p20 = 1.5 * c2 - 0.5;
  457.             p30 = c * (2.5 * c2 - 1.5);
  458.             p40 = 4.375 * c4 - 3.75 * c2 + 0.375;
  459.             p50 = c * (7.875 * c4 - 8.75 * c2 + 1.875);
  460.             p60 = (5.5 * c * p50 - 2.5 * p40) / 3.0;
  461.             p11 = s;
  462.             p21 = 3.0 * c * s;
  463.             p31 = s * (7.5 * c2 - 1.5);
  464.             p41 = c * s * (17.5 * c2 - 7.5);
  465.             p51 = s * (39.375 * c4 - 26.25 * c2 + 1.875);
  466.             p22 = 3.0 * s2;
  467.             p32 = 15.0 * c * s2;
  468.             p42 = s2 * (52.5 * c2 - 7.5);
  469.             p52 = 3.0 * c * p42 - 2.0 * p32;
  470.             p62 = 2.75 * c * p52 - 1.75 * p42;
  471.             p33 = 15.0 * s * s2;

  472.             // compute Legendre polynomials wrt magnetic pole (79N, 71W)
  473.             final double clmlmg = FastMath.cos(lon - XLMG);
  474.             final double cmg  = s * CPMG * clmlmg + c * SPMG;
  475.             final double cmg2 = cmg * cmg;
  476.             final double cmg4 = cmg2 * cmg2;
  477.             p10mg = cmg;
  478.             p20mg = 1.5 * cmg2 - 0.5;
  479.             p40mg = 4.375 * cmg4 - 3.75 * cmg2 + 0.375;

  480.             clfl = scLon.cos();
  481.             slfl = scLon.sin();

  482.             // local time
  483.             hl0 = hl;
  484.             final SinCos scHlo = FastMath.sinCos(hl0);
  485.             ch  = scHlo.cos();
  486.             sh  = scHlo.sin();
  487.             c2h = ch * ch - sh * sh;
  488.             s2h = 2.0 * ch * sh;
  489.             c3h = c2h * ch - s2h * sh;
  490.             s3h = s2h * ch + c2h * sh;

  491.             final double zlb = ZLB0; // + dzlb ??

  492.             final double[] dtt  = new double[tt.length];
  493.             final double[] dh   = new double[tt.length];
  494.             final double[] dhe  = new double[tt.length];
  495.             final double[] dox  = new double[tt.length];
  496.             final double[] daz2 = new double[tt.length];
  497.             final double[] do2  = new double[tt.length];
  498.             final double[] daz  = new double[tt.length];
  499.             final double[] dt0  = new double[tt.length];
  500.             final double[] dtp  = new double[tt.length];

  501.             Arrays.fill(dtt,  Double.NaN);
  502.             Arrays.fill(dh,   Double.NaN);
  503.             Arrays.fill(dhe,  Double.NaN);
  504.             Arrays.fill(dox,  Double.NaN);
  505.             Arrays.fill(daz2, Double.NaN);
  506.             Arrays.fill(do2,  Double.NaN);
  507.             Arrays.fill(daz,  Double.NaN);
  508.             Arrays.fill(dt0,  Double.NaN);
  509.             Arrays.fill(dtp,  Double.NaN);

  510.             //  compute function g(l) / tinf, t120, tp120
  511.             int kleq = 1;
  512.             final double gdelt = gFunction(tt, dtt, 1, kleq);
  513.             dtt[1] = 1.0 + gdelt;
  514.             final double tinf   = tt[1] * dtt[1];

  515.             kleq = 0; // equinox

  516.             if (day < 59 || day > 284) {
  517.                 kleq = -1; // north winter
  518.             }
  519.             if (day > 99 && day < 244) {
  520.                 kleq = 1; // north summer
  521.             }

  522.             final double gdelt0 =  gFunction(t0, dt0, 0, kleq);
  523.             dt0[1] = (t0[1] + gdelt0) / t0[1];
  524.             final double t120 = t0[1] + gdelt0;
  525.             final double gdeltp = gFunction(tp, dtp, 0, kleq);
  526.             dtp[1] = (tp[1] + gdeltp) / tp[1];
  527.             final double tp120 = tp[1] + gdeltp;

  528.             // compute n(z) concentrations: H, He, O, N2, O2, N
  529.             final double sigma   = tp120 / (tinf - t120);
  530.             final double dzeta   = (RE + zlb) / (RE + altiKM);
  531.             final double zeta    = (altiKM - zlb) * dzeta;
  532.             final double sigzeta = sigma * zeta;
  533.             final double expsz   = FastMath.exp(-sigzeta);
  534.             final double tz = tinf - (tinf - t120) * expsz;

  535.             final double[] dbase = new double[7];

  536.             kleq = 1;

  537.             final double gdelh = gFunction(h, dh, 0, kleq);
  538.             dh[1] = FastMath.exp(gdelh);
  539.             dbase[1] = h[1] * dh[1];

  540.             final double gdelhe = gFunction(he, dhe, 0, kleq);
  541.             dhe[1] = FastMath.exp(gdelhe);
  542.             dbase[2] = he[1] * dhe[1];

  543.             final double gdelo = gFunction(o, dox, 1, kleq);
  544.             dox[1] = FastMath.exp(gdelo);
  545.             dbase[3] = o[1] * dox[1];

  546.             final double gdelaz2 = gFunction(az2, daz2, 1, kleq);
  547.             daz2[1] = FastMath.exp(gdelaz2);
  548.             dbase[4] = az2[1] * daz2[1];

  549.             final double gdelo2 = gFunction(o2, do2, 1, kleq);
  550.             do2[1] = FastMath.exp(gdelo2);
  551.             dbase[5] = o2[1] * do2[1];

  552.             final double gdelaz = gFunction(az, daz, 1, kleq);
  553.             daz[1] = FastMath.exp(gdelaz);
  554.             dbase[6] = az[1] * daz[1];

  555.             final double zlbre  = 1.0 + zlb / RE;
  556.             final double glb    = (GSURF / (zlbre * zlbre)) / (sigma * RGAS * tinf);
  557.             final double t120tz = t120 / tz;

  558.             // Partial densities in (g/cm3).
  559.             // d(1) = hydrogen
  560.             // d(2) = helium
  561.             // d(3) = atomic oxygen
  562.             // d(4) = molecular nitrogen
  563.             // d(5) = molecular oxygen
  564.             // d(6) = atomic nitrogen
  565.             double tmpro = 0;
  566.             for (int i = 1; i <= 6; i++) {
  567.                 final double gamma = MA[i] * glb;
  568.                 final double upapg = 1.0 + ALEFA[i] + gamma;
  569.                 final double fzI = FastMath.pow(t120tz, upapg) * FastMath.exp(-sigzeta * gamma);
  570.                 // concentrations of H, He, O, N2, O2, N (particles/cm³)
  571.                 final double ccI = dbase[i] * fzI;
  572.                 // contribution of densities of H, He, O, N2, O2, N (g/cm³)
  573.                 tmpro += ccI * VMA[i];
  574.             }
  575.             this.ro = tmpro;

  576.         }

  577.         /** Computation of function G.
  578.          * @param a vector of coefficients for computation
  579.          * @param da vector of partial derivatives
  580.          * @param ff0 coefficient flag (1 for Ox, Az, He, T°; 0 for H and tp120)
  581.          * @param kle_eq season indicator flag (summer, winter, equinox)
  582.          * @return value of G
  583.          */
  584.         private double gFunction(final double[] a, final double[] da,
  585.                                  final int ff0, final int kle_eq) {

  586.             final double[] fmfb   = new double[3];
  587.             final double[] fbm150 = new double[3];

  588.             // latitude terms
  589.             da[2]  = p20;
  590.             da[3]  = p40;
  591.             da[74] = p10;
  592.             double a74 = a[74];
  593.             double a77 = a[77];
  594.             double a78 = a[78];
  595.             if (kle_eq == -1) {
  596.                 // winter
  597.                 a74 = -a74;
  598.                 a77 = -a77;
  599.                 a78 = -a78;
  600.             }
  601.             if (kle_eq == 0 ) {
  602.                 // equinox
  603.                 a74 = semestrialCorrection(a74);
  604.                 a77 = semestrialCorrection(a77);
  605.                 a78 = semestrialCorrection(a78);
  606.             }
  607.             da[77] = p30;
  608.             da[78] = p50;
  609.             da[79] = p60;

  610.             // flux terms
  611.             fmfb[1]   = f[1] - fbar[1];
  612.             fmfb[2]   = f[2] - fbar[2];
  613.             fbm150[1] = fbar[1] - 150.0;
  614.             fbm150[2] = fbar[2];
  615.             da[4]     = fmfb[1];
  616.             da[6]     = fbm150[1];
  617.             da[4]     = da[4] + a[70] * fmfb[2];
  618.             da[6]     = da[6] + a[71] * fbm150[2];
  619.             da[70]    = fmfb[2] * (a[4] + 2.0 * a[5] * da[4] + a[82] * p10 +
  620.                                    a[83] * p20 + a[84] * p30);
  621.             da[71]    = fbm150[2] * (a[6] + 2.0 * a[69] * da[6] + a[85] * p10 +
  622.                                      a[86] * p20 + a[87] * p30);
  623.             da[5]     = da[4] * da[4];
  624.             da[69]    = da[6] * da[6];
  625.             da[82]    = da[4] * p10;
  626.             da[83]    = da[4] * p20;
  627.             da[84]    = da[4] * p30;
  628.             da[85]    = da[6] * p20;
  629.             da[86]    = da[6] * p30;
  630.             da[87]    = da[6] * p40;

  631.             // Kp terms
  632.             final int ikp  = 62;
  633.             final int ikpm = 67;
  634.             final double c2fi = 1.0 - p10mg * p10mg;
  635.             final double dkp  = akp[1] + (a[ikp] + c2fi * a[ikp + 1]) * akp[2];
  636.             double dakp = a[7] + a[8] * p20mg + a[68] * p40mg +
  637.                           2.0 * dkp * (a[60] + a[61] * p20mg +
  638.                                        a[75] * 2.0 * dkp * dkp);
  639.             da[ikp] = dakp * akp[2];
  640.             da[ikp + 1] = da[ikp] * c2fi;
  641.             final double dkpm  = akp[3] + a[ikpm] * akp[4];
  642.             final double dakpm = a[64] + a[65] * p20mg + a[72] * p40mg +
  643.                                  2.0 * dkpm * (a[66] + a[73] * p20mg +
  644.                                                a[76] * 2.0 * dkpm * dkpm);
  645.             da[ikpm] = dakpm * akp[4];
  646.             da[7]    = dkp;
  647.             da[8]    = p20mg * dkp;
  648.             da[68]   = p40mg * dkp;
  649.             da[60]   = dkp * dkp;
  650.             da[61]   = p20mg * da[60];
  651.             da[75]   = da[60] * da[60];
  652.             da[64]   = dkpm;
  653.             da[65]   = p20mg * dkpm;
  654.             da[72]   = p40mg * dkpm;
  655.             da[66]   = dkpm * dkpm;
  656.             da[73]   = p20mg * da[66];
  657.             da[76]   = da[66] * da[66];

  658.             // non-periodic g(l) function
  659.             double f0 = a[4]  * da[4]  + a[5]  * da[5]  + a[6]  * da[6]  +
  660.                         a[69] * da[69] + a[82] * da[82] + a[83] * da[83] +
  661.                         a[84] * da[84] + a[85] * da[85] + a[86] * da[86] +
  662.                         a[87] * da[87];
  663.             final double f1f = 1.0 + f0 * ff0;

  664.             f0 = f0 + a[2] * da[2] + a[3] * da[3] + a74 * da[74] +
  665.                  a77 * da[77] + a[7] * da[7] + a[8] * da[8] +
  666.                  a[60] * da[60] + a[61] * da[61] + a[68] * da[68] +
  667.                  a[64] * da[64] + a[65] * da[65] + a[66] * da[66] +
  668.                  a[72] * da[72] + a[73] * da[73] + a[75] * da[75] +
  669.                  a[76] * da[76] + a78   * da[78] + a[79] * da[79];
  670. //          termes annuels symetriques en latitude
  671.             da[9]  = FastMath.cos(ROT * (day - a[11]));
  672.             da[10] = p20 * da[9];
  673. //          termes semi-annuels symetriques en latitude
  674.             da[12] = FastMath.cos(ROT2 * (day - a[14]));
  675.             da[13] = p20 * da[12];
  676. //          termes annuels non symetriques en latitude
  677.             final double coste = FastMath.cos(ROT * (day - a[18]));
  678.             da[15] = p10 * coste;
  679.             da[16] = p30 * coste;
  680.             da[17] = p50 * coste;
  681. //          terme  semi-annuel  non symetrique  en latitude
  682.             final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
  683.             da[19] = p10 * cos2te;
  684.             da[39] = p30 * cos2te;
  685.             da[59] = p50 * cos2te;
  686. //          termes diurnes [et couples annuel]
  687.             da[21] = p11 * ch;
  688.             da[22] = p31 * ch;
  689.             da[23] = p51 * ch;
  690.             da[24] = da[21] * coste;
  691.             da[25] = p21 * ch * coste;
  692.             da[26] = p11 * sh;
  693.             da[27] = p31 * sh;
  694.             da[28] = p51 * sh;
  695.             da[29] = da[26] * coste;
  696.             da[30] = p21 * sh * coste;
  697. //          termes semi-diurnes [et couples annuel]
  698.             da[31] = p22 * c2h;
  699.             da[37] = p42 * c2h;
  700.             da[32] = p32 * c2h * coste;
  701.             da[33] = p22 * s2h;
  702.             da[38] = p42 * s2h;
  703.             da[34] = p32 * s2h * coste;
  704.             da[88] = p32 * c2h;
  705.             da[89] = p32 * s2h;
  706.             da[90] = p52 * c2h;
  707.             da[91] = p52 * s2h;
  708.             double a88 = a[88];
  709.             double a89 = a[89];
  710.             double a90 = a[90];
  711.             double a91 = a[91];
  712.             if (kle_eq == -1) {            //hiver
  713.                 a88 = -a88;
  714.                 a89 = -a89;
  715.                 a90 = -a90;
  716.                 a91 = -a91;
  717.             }
  718.             if (kle_eq == 0) {             //equinox
  719.                 a88 = semestrialCorrection(a88);
  720.                 a89 = semestrialCorrection(a89);
  721.                 a90 = semestrialCorrection(a90);
  722.                 a91 = semestrialCorrection(a91);
  723.             }
  724.             da[92] = p62 * c2h;
  725.             da[93] = p62 * s2h;
  726. //          termes ter-diurnes
  727.             da[35] = p33 * c3h;
  728.             da[36] = p33 * s3h;
  729. //          fonction g[l] periodique
  730.             double fp = a[9]  * da[9]  + a[10] * da[10] + a[12] * da[12] + a[13] * da[13] +
  731.                         a[15] * da[15] + a[16] * da[16] + a[17] * da[17] + a[19] * da[19] +
  732.                         a[21] * da[21] + a[22] * da[22] + a[23] * da[23] + a[24] * da[24] +
  733.                         a[25] * da[25] + a[26] * da[26] + a[27] * da[27] + a[28] * da[28] +
  734.                         a[29] * da[29] + a[30] * da[30] + a[31] * da[31] + a[32] * da[32] +
  735.                         a[33] * da[33] + a[34] * da[34] + a[35] * da[35] + a[36] * da[36] +
  736.                         a[37] * da[37] + a[38] * da[38] + a[39] * da[39] + a[59] * da[59] +
  737.                         a88   * da[88] + a89   * da[89] + a90   * da[90] + a91   * da[91] +
  738.                         a[92] * da[92] + a[93] * da[93];
  739. //          termes d'activite magnetique
  740.             da[40] = p10 * coste * dkp;
  741.             da[41] = p30 * coste * dkp;
  742.             da[42] = p50 * coste * dkp;
  743.             da[43] = p11 * ch * dkp;
  744.             da[44] = p31 * ch * dkp;
  745.             da[45] = p51 * ch * dkp;
  746.             da[46] = p11 * sh * dkp;
  747.             da[47] = p31 * sh * dkp;
  748.             da[48] = p51 * sh * dkp;

  749. //          fonction g[l] periodique supplementaire
  750.             fp += a[40] * da[40] + a[41] * da[41] + a[42] * da[42] + a[43] * da[43] +
  751.                   a[44] * da[44] + a[45] * da[45] + a[46] * da[46] + a[47] * da[47] +
  752.                   a[48] * da[48];

  753.             dakp = (a[40] * p10 + a[41] * p30 + a[42] * p50) * coste +
  754.                    (a[43] * p11 + a[44] * p31 + a[45] * p51) * ch +
  755.                    (a[46] * p11 + a[47] * p31 + a[48] * p51) * sh;
  756.             da[ikp] += dakp * akp[2];
  757.             da[ikp + 1] = da[ikp] + dakp * c2fi * akp[2];
  758. //          termes de longitude
  759.             da[49] = p11 * clfl;
  760.             da[50] = p21 * clfl;
  761.             da[51] = p31 * clfl;
  762.             da[52] = p41 * clfl;
  763.             da[53] = p51 * clfl;
  764.             da[54] = p11 * slfl;
  765.             da[55] = p21 * slfl;
  766.             da[56] = p31 * slfl;
  767.             da[57] = p41 * slfl;
  768.             da[58] = p51 * slfl;

  769. //          fonction g[l] periodique supplementaire
  770.             fp += a[49] * da[49] + a[50] * da[50] + a[51] * da[51] + a[52] * da[52] +
  771.                   a[53] * da[53] + a[54] * da[54] + a[55] * da[55] + a[56] * da[56] +
  772.                   a[57] * da[57] + a[58] * da[58];

  773. //          fonction g(l) totale (couplage avec le flux)
  774.             return f0 + fp * f1f;

  775.         }


  776.         /** Apply a correction coefficient to the given parameter.
  777.          * @param param the parameter to correct
  778.          * @return the corrected parameter
  779.          */
  780.         private double semestrialCorrection(final double param) {
  781.             final int debeq_pr = 59;
  782.             final int debeq_au = 244;
  783.             final double result;
  784.             if (day >= 100) {
  785.                 final double xmult  = (day - debeq_au) / 40.0;
  786.                 result = param - 2.0 * param * xmult;
  787.             } else {
  788.                 final double xmult  = (day - debeq_pr) / 40.0;
  789.                 result = 2.0 * param * xmult - param;
  790.             }
  791.             return result;
  792.         }


  793.     }

  794.     /** Local holder for intermediate results ensuring the model is reentrant.
  795.      * @param <T> type of the field elements
  796.      */
  797.     private static class FieldComputation<T extends CalculusFieldElement<T>> {

  798.         /** Number of days in current year. */
  799.         private int day;

  800.         /** Instant solar flux. f[1] = instantaneous flux; f[2] = 0. (not used). */
  801.         private double[] f = new double[3];

  802.         /** Mean solar flux. fbar[1] = mean flux; fbar[2] = 0. (not used). */
  803.         private double[] fbar = new double[3];

  804.         /** Kp coefficients.
  805.          * <ul>
  806.          *   <li>akp[1] = 3-hourly kp</li>
  807.          *   <li>akp[2] = 0 (not used)</li>
  808.          *   <li>akp[3] = mean kp of last 24 hours</li>
  809.          *   <li>akp[4] = 0 (not used)</li>
  810.          * </ul>
  811.          */
  812.         private double[] akp = new double[5];

  813.         /** Cosine of the longitude. */
  814.         private final T clfl;

  815.         /** Sine of the longitude. */
  816.         private final T slfl;

  817.         /** Total density (g/cm3). */
  818.         private final T ro;

  819.         // CHECKSTYLE: stop JavadocVariable check

  820.         /** Legendre coefficients. */
  821.         private final T p10;
  822.         private final T p20;
  823.         private final T p30;
  824.         private final T p40;
  825.         private final T p50;
  826.         private final T p60;
  827.         private final T p11;
  828.         private final T p21;
  829.         private final T p31;
  830.         private final T p41;
  831.         private final T p51;
  832.         private final T p22;
  833.         private final T p32;
  834.         private final T p42;
  835.         private final T p52;
  836.         private final T p62;
  837.         private final T p33;
  838.         private final T p10mg;
  839.         private final T p20mg;
  840.         private final T p40mg;

  841.         /** Local time intermediate values. */
  842.         private final T hl0;
  843.         private final T ch;
  844.         private final T sh;
  845.         private final T c2h;
  846.         private final T s2h;
  847.         private final T c3h;
  848.         private final T s3h;

  849.         /** Simple constructor.
  850.          * @param day day of year
  851.          * @param altiKM altitude <em>in kilometers</em>
  852.          * @param lon local longitude (rad)
  853.          * @param lat local latitude (rad)
  854.          * @param hl local solar time in rad (O hr = 0 rad)
  855.          * @param f instantaneous solar flux (F10.7)
  856.          * @param far mean solar flux (F10.7)
  857.          * @param akp geomagnetic activity index
  858.          */
  859.         FieldComputation(final int day,
  860.                          final T altiKM, final T lon, final T lat,
  861.                          final T hl, final double[] f, final double[] far,
  862.                          final double[] akp) {

  863.             this.day  = day;
  864.             this.f    = f;
  865.             this.fbar = far;
  866.             this.akp  = akp;

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

  870.             // compute Legendre polynomials wrt geographic pole
  871.             final T c = scLat.sin();
  872.             final T c2 = c.multiply(c);
  873.             final T c4 = c2.multiply(c2);
  874.             final T s = scLat.cos();
  875.             final T s2 = s.multiply(s);
  876.             p10 = c;
  877.             p20 = c2.multiply(1.5).subtract(0.5);
  878.             p30 = c.multiply(c2.multiply(2.5).subtract(1.5));
  879.             p40 = c4.multiply(4.375).subtract(c2.multiply(3.75)).add(0.375);
  880.             p50 = c.multiply(c4.multiply(7.875).subtract(c2.multiply(8.75)).add(1.875));
  881.             p60 = (c.multiply(5.5).multiply(p50).subtract(p40.multiply(2.5))).divide(3.0);
  882.             p11 = s;
  883.             p21 = c.multiply(3.0).multiply(s);
  884.             p31 = s.multiply(c2.multiply(7.5).subtract(1.5));
  885.             p41 = c.multiply(s).multiply(c2.multiply(17.5).subtract(7.5));
  886.             p51 = s.multiply(c4.multiply(39.375).subtract(c2.multiply(26.25)).add(1.875));
  887.             p22 = s2.multiply(3.0);
  888.             p32 = c.multiply(15.0).multiply(s2);
  889.             p42 = s2.multiply(c2.multiply(52.5).subtract(7.5));
  890.             p52 = c.multiply(3.0).multiply(p42).subtract(p32.multiply(2.0));
  891.             p62 = c.multiply(2.75).multiply(p52).subtract(p42.multiply(1.75));
  892.             p33 = s.multiply(15.0).multiply(s2);

  893.             // compute Legendre polynomials wrt magnetic pole (79N, 71W)
  894.             final T clmlmg = lon.subtract(XLMG).cos();
  895.             final T cmg  = s.multiply(CPMG).multiply(clmlmg).add(c.multiply(SPMG));
  896.             final T cmg2 = cmg.multiply(cmg);
  897.             final T cmg4 = cmg2.multiply(cmg2);
  898.             p10mg = cmg;
  899.             p20mg = cmg2.multiply(1.5).subtract(0.5);
  900.             p40mg = cmg4.multiply(4.375).subtract(cmg2.multiply(3.75)).add(0.375);

  901.             clfl = scLon.cos();
  902.             slfl = scLon.sin();

  903.             // local time
  904.             hl0 = hl;
  905.             final FieldSinCos<T> scHlo = FastMath.sinCos(hl0);
  906.             ch  = scHlo.cos();
  907.             sh  = scHlo.sin();
  908.             c2h = ch.multiply(ch).subtract(sh.multiply(sh));
  909.             s2h = ch.multiply(sh).multiply(2);
  910.             c3h = c2h.multiply(ch).subtract(s2h.multiply(sh));
  911.             s3h = s2h.multiply(ch).add(c2h.multiply(sh));

  912.             final double zlb = ZLB0; // + dzlb ??

  913.             final T[] dtt  = MathArrays.buildArray(altiKM.getField(), tt.length);
  914.             final T[] dh   = MathArrays.buildArray(altiKM.getField(), tt.length);
  915.             final T[] dhe  = MathArrays.buildArray(altiKM.getField(), tt.length);
  916.             final T[] dox  = MathArrays.buildArray(altiKM.getField(), tt.length);
  917.             final T[] daz2 = MathArrays.buildArray(altiKM.getField(), tt.length);
  918.             final T[] do2  = MathArrays.buildArray(altiKM.getField(), tt.length);
  919.             final T[] daz  = MathArrays.buildArray(altiKM.getField(), tt.length);
  920.             final T[] dt0  = MathArrays.buildArray(altiKM.getField(), tt.length);
  921.             final T[] dtp  = MathArrays.buildArray(altiKM.getField(), tt.length);

  922.             //  compute function g(l) / tinf, t120, tp120
  923.             int kleq = 1;
  924.             final T gdelt = gFunction(tt, dtt, 1, kleq);
  925.             dtt[1] = gdelt.add(1);
  926.             final T tinf   = dtt[1].multiply(tt[1]);

  927.             kleq = 0; // equinox

  928.             if (day < 59 || day > 284) {
  929.                 kleq = -1; // north winter
  930.             }
  931.             if (day > 99 && day < 244) {
  932.                 kleq = 1; // north summer
  933.             }

  934.             final T gdelt0 =  gFunction(t0, dt0, 0, kleq);
  935.             dt0[1] = gdelt0.add(t0[1]).divide(t0[1]);
  936.             final T t120 = gdelt0.add(t0[1]);
  937.             final T gdeltp = gFunction(tp, dtp, 0, kleq);
  938.             dtp[1] = gdeltp.add(tp[1]).divide(tp[1]);
  939.             final T tp120 = gdeltp.add(tp[1]);

  940.             // compute n(z) concentrations: H, He, O, N2, O2, N
  941.             final T sigma   = tp120.divide(tinf.subtract(t120));
  942.             final T dzeta   = altiKM.add(RE).reciprocal().multiply(zlb + RE);
  943.             final T zeta    = altiKM.subtract(zlb).multiply(dzeta);
  944.             final T sigzeta = sigma.multiply(zeta);
  945.             final T expsz   = sigzeta.negate().exp();
  946.             final T tz = tinf.subtract(tinf.subtract(t120).multiply(expsz));

  947.             final T[] dbase = MathArrays.buildArray(altiKM.getField(), 7);

  948.             kleq = 1;

  949.             final T gdelh = gFunction(h, dh, 0, kleq);
  950.             dh[1] = gdelh.exp();
  951.             dbase[1] = dh[1].multiply(h[1]);

  952.             final T gdelhe = gFunction(he, dhe, 0, kleq);
  953.             dhe[1] = gdelhe.exp();
  954.             dbase[2] = dhe[1].multiply(he[1]);

  955.             final T gdelo = gFunction(o, dox, 1, kleq);
  956.             dox[1] = gdelo.exp();
  957.             dbase[3] = dox[1].multiply(o[1]);

  958.             final T gdelaz2 = gFunction(az2, daz2, 1, kleq);
  959.             daz2[1] = gdelaz2.exp();
  960.             dbase[4] = daz2[1].multiply(az2[1]);

  961.             final T gdelo2 = gFunction(o2, do2, 1, kleq);
  962.             do2[1] = gdelo2.exp();
  963.             dbase[5] = do2[1].multiply(o2[1]);

  964.             final T gdelaz = gFunction(az, daz, 1, kleq);
  965.             daz[1] = gdelaz.exp();
  966.             dbase[6] = daz[1].multiply(az[1]);

  967.             final double zlbre  = 1.0 + zlb / RE;
  968.             final T glb    = sigma.multiply(RGAS).multiply(tinf).reciprocal().multiply(GSURF / (zlbre * zlbre));
  969.             final T t120tz = t120.divide(tz);

  970.             // Partial densities in (g/cm3).
  971.             // d(1) = hydrogen
  972.             // d(2) = helium
  973.             // d(3) = atomic oxygen
  974.             // d(4) = molecular nitrogen
  975.             // d(5) = molecular oxygen
  976.             // d(6) = atomic nitrogen
  977.             T tmpro = altiKM.getField().getZero();
  978.             for (int i = 1; i <= 6; i++) {
  979.                 final T gamma = glb.multiply(MA[i]);
  980.                 final T upapg = gamma.add(1.0 + ALEFA[i]);
  981.                 final T fzI = t120tz.pow(upapg).multiply(sigzeta.negate().multiply(gamma).exp());
  982.                 // concentrations of H, He, O, N2, O2, N (particles/cm³)
  983.                 final T ccI = dbase[i].multiply(fzI);
  984.                 // contribution of densities of H, He, O, N2, O2, N (g/cm³)
  985.                 tmpro = tmpro.add(ccI.multiply(VMA[i]));
  986.             }
  987.             this.ro = tmpro;

  988.         }

  989.         /** Computation of function G.
  990.          * @param a vector of coefficients for computation
  991.          * @param da vector of partial derivatives
  992.          * @param ff0 coefficient flag (1 for Ox, Az, He, T°; 0 for H and tp120)
  993.          * @param kle_eq season indicator flag (summer, winter, equinox)
  994.          * @return value of G
  995.          */
  996.         private T gFunction(final double[] a, final T[] da,
  997.                             final int ff0, final int kle_eq) {

  998.             final T zero = da[0].getField().getZero();
  999.             final double[] fmfb   = new double[3];
  1000.             final double[] fbm150 = new double[3];

  1001.             // latitude terms
  1002.             da[2]  = p20;
  1003.             da[3]  = p40;
  1004.             da[74] = p10;
  1005.             double a74 = a[74];
  1006.             double a77 = a[77];
  1007.             double a78 = a[78];
  1008.             if (kle_eq == -1) {
  1009.                 // winter
  1010.                 a74 = -a74;
  1011.                 a77 = -a77;
  1012.                 a78 = -a78;
  1013.             }
  1014.             if (kle_eq == 0 ) {
  1015.                 // equinox
  1016.                 a74 = semestrialCorrection(a74);
  1017.                 a77 = semestrialCorrection(a77);
  1018.                 a78 = semestrialCorrection(a78);
  1019.             }
  1020.             da[77] = p30;
  1021.             da[78] = p50;
  1022.             da[79] = p60;

  1023.             // flux terms
  1024.             fmfb[1]   = f[1] - fbar[1];
  1025.             fmfb[2]   = f[2] - fbar[2];
  1026.             fbm150[1] = fbar[1] - 150.0;
  1027.             fbm150[2] = fbar[2];
  1028.             da[4]     = zero.add(fmfb[1]);
  1029.             da[6]     = zero.add(fbm150[1]);
  1030.             da[4]     = da[4].add(a[70] * fmfb[2]);
  1031.             da[6]     = da[6].add(a[71] * fbm150[2]);
  1032.             da[70]    = da[4].multiply(a[ 5]).multiply(2).
  1033.                             add(p10.multiply(a[82])).
  1034.                             add(p20.multiply(a[83])).
  1035.                             add(p30.multiply(a[84])).
  1036.                             add(a[4]).
  1037.                         multiply(fmfb[2]);
  1038.             da[71]    = da[6].multiply(a[69]).multiply(2).
  1039.                             add(p10.multiply(a[85])).
  1040.                             add(p20.multiply(a[86])).
  1041.                             add(p30.multiply(a[87])).
  1042.                             add(a[6]).
  1043.                         multiply(fbm150[2]);
  1044.             da[5]     = da[4].multiply(da[4]);
  1045.             da[69]    = da[6].multiply(da[6]);
  1046.             da[82]    = da[4].multiply(p10);
  1047.             da[83]    = da[4].multiply(p20);
  1048.             da[84]    = da[4].multiply(p30);
  1049.             da[85]    = da[6].multiply(p20);
  1050.             da[86]    = da[6].multiply(p30);
  1051.             da[87]    = da[6].multiply(p40);

  1052.             // Kp terms
  1053.             final int ikp  = 62;
  1054.             final int ikpm = 67;
  1055.             final T c2fi = p10mg.multiply(p10mg).negate().add(1);
  1056.             final T dkp  = c2fi.multiply(a[ikp + 1]).add(a[ikp]).multiply(akp[2]).add(akp[1]);
  1057.             T dakp = p20mg.multiply(a[8]).add(p40mg.multiply(a[68])).
  1058.                      add(p20mg.multiply(a[61]).add(dkp.multiply(dkp).multiply(2 * a[75]).add(a[60])).multiply(dkp.multiply(2))).
  1059.                      add(a[7]);
  1060.             da[ikp]     = dakp.multiply(akp[2]);
  1061.             da[ikp + 1] = da[ikp].multiply(c2fi);
  1062.             final double dkpm  = akp[3] + a[ikpm] * akp[4];
  1063.             final T dakpm = p20mg.multiply(a[65]).add(p40mg.multiply(a[72])).
  1064.                             add(p20mg.multiply(a[73]).add(a[66] + a[76] * 2.0 * dkpm * dkpm).multiply( 2.0 * dkpm)).
  1065.                             add(a[64]);
  1066.             da[ikpm] = dakpm.multiply(akp[4]);
  1067.             da[7]    = dkp;
  1068.             da[8]    = p20mg.multiply(dkp);
  1069.             da[68]   = p40mg.multiply(dkp);
  1070.             da[60]   = dkp.multiply(dkp);
  1071.             da[61]   = p20mg.multiply(da[60]);
  1072.             da[75]   = da[60].multiply(da[60]);
  1073.             da[64]   = zero.add(dkpm);
  1074.             da[65]   = p20mg.multiply(dkpm);
  1075.             da[72]   = p40mg.multiply(dkpm);
  1076.             da[66]   = zero.add(dkpm * dkpm);
  1077.             da[73]   = p20mg.multiply(da[66]);
  1078.             da[76]   = da[66].multiply(da[66]);

  1079.             // non-periodic g(l) function
  1080.             T f0 = da[4].multiply(a[4]).
  1081.                    add(da[5].multiply(a[5])).
  1082.                    add(da[6].multiply(a[6])).
  1083.                    add(da[69].multiply(a[69])).
  1084.                    add(da[82].multiply(a[82])).
  1085.                    add(da[83].multiply(a[83])).
  1086.                    add(da[84].multiply(a[84])).
  1087.                    add(da[85].multiply(a[85])).
  1088.                    add(da[86].multiply(a[86])).
  1089.                    add(da[87].multiply(a[87]));
  1090.             final T f1f = f0.multiply(ff0).add(1);

  1091.             f0 = f0.
  1092.                  add(da[2].multiply(a[2])).
  1093.                  add(da[3].multiply(a[3])).
  1094.                  add(da[7].multiply(a[7])).
  1095.                  add(da[8].multiply(a[8])).
  1096.                  add(da[60].multiply(a[60])).
  1097.                  add(da[61].multiply(a[61])).
  1098.                  add(da[68].multiply(a[68])).
  1099.                  add(da[64].multiply(a[64])).
  1100.                  add(da[65].multiply(a[65])).
  1101.                  add(da[66].multiply(a[66])).
  1102.                  add(da[72].multiply(a[72])).
  1103.                  add(da[73].multiply(a[73])).
  1104.                  add(da[74].multiply(a74)).
  1105.                  add(da[75].multiply(a[75])).
  1106.                  add(da[76].multiply(a[76])).
  1107.                  add(da[77].multiply(a77)).
  1108.                  add(da[78].multiply(a78)).
  1109.                  add(da[79].multiply(a[79]));
  1110. //          termes annuels symetriques en latitude
  1111.             da[9]  = zero.add(FastMath.cos(ROT * (day - a[11])));
  1112.             da[10] = p20.multiply(da[9]);
  1113. //          termes semi-annuels symetriques en latitude
  1114.             da[12] = zero.add(FastMath.cos(ROT2 * (day - a[14])));
  1115.             da[13] = p20.multiply(da[12]);
  1116. //          termes annuels non symetriques en latitude
  1117.             final double coste = FastMath.cos(ROT * (day - a[18]));
  1118.             da[15] = p10.multiply(coste);
  1119.             da[16] = p30.multiply(coste);
  1120.             da[17] = p50.multiply(coste);
  1121. //          terme  semi-annuel  non symetrique  en latitude
  1122.             final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
  1123.             da[19] = p10.multiply(cos2te);
  1124.             da[39] = p30.multiply(cos2te);
  1125.             da[59] = p50.multiply(cos2te);
  1126. //          termes diurnes [et couples annuel]
  1127.             da[21] = p11.multiply(ch);
  1128.             da[22] = p31.multiply(ch);
  1129.             da[23] = p51.multiply(ch);
  1130.             da[24] = da[21].multiply(coste);
  1131.             da[25] = p21.multiply(ch).multiply(coste);
  1132.             da[26] = p11.multiply(sh);
  1133.             da[27] = p31.multiply(sh);
  1134.             da[28] = p51.multiply(sh);
  1135.             da[29] = da[26].multiply(coste);
  1136.             da[30] = p21.multiply(sh).multiply(coste);
  1137. //          termes semi-diurnes [et couples annuel]
  1138.             da[31] = p22.multiply(c2h);
  1139.             da[37] = p42.multiply(c2h);
  1140.             da[32] = p32.multiply(c2h).multiply(coste);
  1141.             da[33] = p22.multiply(s2h);
  1142.             da[38] = p42.multiply(s2h);
  1143.             da[34] = p32.multiply(s2h).multiply(coste);
  1144.             da[88] = p32.multiply(c2h);
  1145.             da[89] = p32.multiply(s2h);
  1146.             da[90] = p52.multiply(c2h);
  1147.             da[91] = p52.multiply(s2h);
  1148.             double a88 = a[88];
  1149.             double a89 = a[89];
  1150.             double a90 = a[90];
  1151.             double a91 = a[91];
  1152.             if (kle_eq == -1) {            //hiver
  1153.                 a88 = -a88;
  1154.                 a89 = -a89;
  1155.                 a90 = -a90;
  1156.                 a91 = -a91;
  1157.             }
  1158.             if (kle_eq == 0) {             //equinox
  1159.                 a88 = semestrialCorrection(a88);
  1160.                 a89 = semestrialCorrection(a89);
  1161.                 a90 = semestrialCorrection(a90);
  1162.                 a91 = semestrialCorrection(a91);
  1163.             }
  1164.             da[92] = p62.multiply(c2h);
  1165.             da[93] = p62.multiply(s2h);
  1166. //          termes ter-diurnes
  1167.             da[35] = p33.multiply(c3h);
  1168.             da[36] = p33.multiply(s3h);
  1169. //          fonction g[l] periodique
  1170.             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])).
  1171.                    add(da[15].multiply(a[15])).add(da[16].multiply(a[16])).add(da[17].multiply(a[17])).add(da[19].multiply(a[19])).
  1172.                    add(da[21].multiply(a[21])).add(da[22].multiply(a[22])).add(da[23].multiply(a[23])).add(da[24].multiply(a[24])).
  1173.                    add(da[25].multiply(a[25])).add(da[26].multiply(a[26])).add(da[27].multiply(a[27])).add(da[28].multiply(a[28])).
  1174.                    add(da[29].multiply(a[29])).add(da[30].multiply(a[30])).add(da[31].multiply(a[31])).add(da[32].multiply(a[32])).
  1175.                    add(da[33].multiply(a[33])).add(da[34].multiply(a[34])).add(da[35].multiply(a[35])).add(da[36].multiply(a[36])).
  1176.                    add(da[37].multiply(a[37])).add(da[38].multiply(a[38])).add(da[39].multiply(a[39])).add(da[59].multiply(a[59])).
  1177.                    add(da[88].multiply(a88))  .add(da[89].multiply(a89)  ).add(da[90].multiply(a90)  ).add(da[91].multiply(a91)  ).
  1178.                    add(da[92].multiply(a[92])).add(da[93].multiply(a[93]));
  1179. //          termes d'activite magnetique
  1180.             da[40] = p10.multiply(coste).multiply(dkp);
  1181.             da[41] = p30.multiply(coste).multiply(dkp);
  1182.             da[42] = p50.multiply(coste).multiply(dkp);
  1183.             da[43] = p11.multiply(ch).multiply(dkp);
  1184.             da[44] = p31.multiply(ch).multiply(dkp);
  1185.             da[45] = p51.multiply(ch).multiply(dkp);
  1186.             da[46] = p11.multiply(sh).multiply(dkp);
  1187.             da[47] = p31.multiply(sh).multiply(dkp);
  1188.             da[48] = p51.multiply(sh).multiply(dkp);

  1189. //          fonction g[l] periodique supplementaire
  1190.             fp = fp.
  1191.                   add(da[40].multiply(a[40])).
  1192.                   add(da[41].multiply(a[41])).
  1193.                   add(da[42].multiply(a[42])).
  1194.                   add(da[43].multiply(a[43])).
  1195.                   add(da[44].multiply(a[44])).
  1196.                   add(da[45].multiply(a[45])).
  1197.                   add(da[46].multiply(a[46])).
  1198.                   add(da[47].multiply(a[47])).
  1199.                   add(da[48].multiply(a[48]));

  1200.             dakp =     p10.multiply(a[40]).add(p30.multiply(a[41])).add(p50.multiply(a[42])).multiply(coste).
  1201.                    add(p11.multiply(a[40]).add(p31.multiply(a[44])).add(p51.multiply(a[45])).multiply(ch)).
  1202.                    add(p11.multiply(a[46]).add(p31.multiply(a[47])).add(p51.multiply(a[48])).multiply(sh));
  1203.             da[ikp] = da[ikp].add(dakp.multiply(akp[2]));
  1204.             da[ikp + 1] = da[ikp].add(dakp.multiply(c2fi).multiply(akp[2]));
  1205. //          termes de longitude
  1206.             da[49] = p11.multiply(clfl);
  1207.             da[50] = p21.multiply(clfl);
  1208.             da[51] = p31.multiply(clfl);
  1209.             da[52] = p41.multiply(clfl);
  1210.             da[53] = p51.multiply(clfl);
  1211.             da[54] = p11.multiply(slfl);
  1212.             da[55] = p21.multiply(slfl);
  1213.             da[56] = p31.multiply(slfl);
  1214.             da[57] = p41.multiply(slfl);
  1215.             da[58] = p51.multiply(slfl);

  1216. //          fonction g[l] periodique supplementaire
  1217.             fp = fp.
  1218.                  add(da[49].multiply(a[49])).
  1219.                  add(da[50].multiply(a[50])).
  1220.                  add(da[51].multiply(a[51])).
  1221.                  add(da[52].multiply(a[52])).
  1222.                  add(da[53].multiply(a[53])).
  1223.                  add(da[54].multiply(a[54])).
  1224.                  add(da[55].multiply(a[55])).
  1225.                  add(da[56].multiply(a[56])).
  1226.                  add(da[57].multiply(a[57])).
  1227.                  add(da[58].multiply(a[58]));

  1228. //          fonction g(l) totale (couplage avec le flux)
  1229.             return f0.add(fp.multiply(f1f));

  1230.         }


  1231.         /** Apply a correction coefficient to the given parameter.
  1232.          * @param param the parameter to correct
  1233.          * @return the corrected parameter
  1234.          */
  1235.         private double semestrialCorrection(final double param) {
  1236.             final int debeq_pr = 59;
  1237.             final int debeq_au = 244;
  1238.             final double result;
  1239.             if (day >= 100) {
  1240.                 final double xmult  = (day - debeq_au) / 40.0;
  1241.                 result = param - 2.0 * param * xmult;
  1242.             } else {
  1243.                 final double xmult  = (day - debeq_pr) / 40.0;
  1244.                 result = 2.0 * param * xmult - param;
  1245.             }
  1246.             return result;
  1247.         }

  1248.     }

  1249. }