DTM2000.java

  1. /* Copyright 2002-2018 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.forces.drag.atmosphere;

  19. import java.io.BufferedReader;
  20. import java.io.IOException;
  21. import java.io.InputStream;
  22. import java.io.InputStreamReader;
  23. import java.util.Arrays;
  24. import java.util.Calendar;
  25. import java.util.GregorianCalendar;

  27. import org.hipparchus.RealFieldElement;
  28. import org.hipparchus.exception.DummyLocalizable;
  29. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  30. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  31. import org.hipparchus.util.FastMath;
  32. import org.hipparchus.util.MathArrays;
  33. import org.orekit.bodies.BodyShape;
  34. import org.orekit.bodies.FieldGeodeticPoint;
  35. import org.orekit.bodies.GeodeticPoint;
  36. import org.orekit.errors.OrekitException;
  37. import org.orekit.errors.OrekitMessages;
  38. import org.orekit.frames.Frame;
  39. import org.orekit.time.AbsoluteDate;
  40. import org.orekit.time.FieldAbsoluteDate;
  41. import org.orekit.time.TimeScalesFactory;
  42. 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="http://sec.noaa.gov/Data/index.html">
  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 {

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

  85.     /** Identifier for helium. */
  86.     public static final int HELIUM = 2;

  88.     /** Identifier for atomic oxygen. */
  89.     public static final int ATOMIC_OXYGEN = 3;

  91.     /** Identifier for molecular nitrogen. */
  92.     public static final int MOLECULAR_NITROGEN = 4;

  94.     /** Identifier for molecular oxygen. */
  95.     public static final int MOLECULAR_OXYGEN = 5;

  97.     /** Identifier for atomic nitrogen. */
  98.     public static final int ATOMIC_NITROGEN = 6;

  100.     /** Serializable UID. */
  101.     private static final long serialVersionUID = 20170705L;

  103.     // Constants :

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

  108.     /** Thermal diffusion coefficient. */
  109.     private static final double[] ALEFA = {
  110.         0, -0.40, -0.38, 0, 0, 0, 0
  111.     };

  113.     /** Atomic mass  H, He, O, N2, O2, N. */
  114.     private static final double[] MA = {
  115.         0, 1, 4, 16, 28, 32, 14
  116.     };

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

  123.     /** Polar Earth radius. */
  124.     private static final double RE = 6356.77;

  126.     /** Reference altitude. */
  127.     private static final double ZLB0 = 120.0;

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

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

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

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

  141.     /** Universal gas constant. */
  142.     private static final double RGAS = 831.4;

  144.     /** 2 * π / 365. */
  145.     private static final double ROT = 0.017214206;

  147.     /** 2 * rot. */
  148.     private static final double ROT2 = 0.034428412;

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

  153.     // CHECKSTYLE: stop JavadocVariable check

  155.     /** Elements coefficients. */
  156.     private static double[] tt   = null;
  157.     private static double[] h    = null;
  158.     private static double[] he   = null;
  159.     private static double[] o    = null;
  160.     private static double[] az2  = null;
  161.     private static double[] o2   = null;
  162.     private static double[] az   = null;
  163.     private static double[] t0   = null;
  164.     private static double[] tp   = null;

  166.     /** Sun position. */
  167.     private PVCoordinatesProvider sun;

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

  172.     /** Earth body shape. */
  173.     private BodyShape earth;

  175.     /** Simple constructor for independent computation.
  176.      * @param parameters the solar and magnetic activity data
  177.      * @param sun the sun position
  178.      * @param earth the earth body shape
  179.      * @exception OrekitException if some resource file reading error occurs
  180.      */
  181.     public DTM2000(final DTM2000InputParameters parameters,
  182.                    final PVCoordinatesProvider sun, final BodyShape earth)
  183.         throws OrekitException {

  185.         synchronized (DTM2000.class) {
  186.             // lazy reading of model coefficients
  187.             if (tt == null) {
  188.                 readcoefficients();
  189.             }
  190.         }

  192.         this.earth = earth;
  193.         this.sun = sun;
  194.         this.inputParams = parameters;

  196.     }

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

  203.     /** Get the local density with initial entries.
  204.      * @param day day of year
  205.      * @param alti altitude in meters
  206.      * @param lon local longitude (rad)
  207.      * @param lat local latitude (rad)
  208.      * @param hl local solar time in rad (O hr = 0 rad)
  209.      * @param f instantaneous solar flux (F10.7)
  210.      * @param fbar mean solar flux (F10.7)
  211.      * @param akp3 3 hrs geomagnetic activity index (1-9)
  212.      * @param akp24 Mean of last 24 hrs geomagnetic activity index (1-9)
  213.      * @return the local density (kg/m³)
  214.      * @exception OrekitException if altitude is outside of supported range
  215.      */
  216.     public double getDensity(final int day,
  217.                              final double alti, final double lon, final double lat,
  218.                              final double hl, final double f, final double fbar,
  219.                              final double akp3, final double akp24)
  220.         throws OrekitException {
  221.         final double threshold = 120000;
  222.         if (alti < threshold) {
  223.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
  224.                                       alti, threshold);
  225.         }
  226.         final Computation result = new Computation(day, alti / 1000, lon, lat, hl,
  227.                                                    new double[] {
  228.                                                        0, f, 0
  229.                                                    }, new double[] {
  230.                                                        0, fbar, 0
  231.                                                    }, new double[] {
  232.                                                        0, akp3, 0, akp24, 0
  233.                                                    });
  234.         return result.ro * 1000;
  235.     }

  237.     /** Get the local density with initial entries.
  238.      * @param day day of year
  239.      * @param alti altitude in meters
  240.      * @param lon local longitude (rad)
  241.      * @param lat local latitude (rad)
  242.      * @param hl local solar time in rad (O hr = 0 rad)
  243.      * @param f instantaneous solar flux (F10.7)
  244.      * @param fbar mean solar flux (F10.7)
  245.      * @param akp3 3 hrs geomagnetic activity index (1-9)
  246.      * @param akp24 Mean of last 24 hrs geomagnetic activity index (1-9)
  247.      * @param <T> type of the field elements
  248.      * @return the local density (kg/m³)
  249.      * @exception OrekitException if altitude is outside of supported range
  250.      * @since 9.0
  251.      */
  252.     public <T extends RealFieldElement<T>> T getDensity(final int day,
  253.                                                         final T alti, final T lon, final T lat,
  254.                                                         final T hl, final double f, final double fbar,
  255.                                                         final double akp3, final double akp24)
  256.         throws OrekitException {
  257.         final double threshold = 120000;
  258.         if (alti.getReal() < threshold) {
  259.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
  260.                                       alti, threshold);
  261.         }
  262.         final FieldComputation<T> result = new FieldComputation<>(day, alti.divide(1000), lon, lat, hl,
  263.                                                                   new double[] {
  264.                                                                       0, f, 0
  265.                                                                   }, new double[] {
  266.                                                                       0, fbar, 0
  267.                                                                   }, new double[] {
  268.                                                                       0, akp3, 0, akp24, 0
  269.                                                                   });
  270.         return result.ro.multiply(1000);
  271.     }

  273.     /** Store the DTM model elements coefficients in internal arrays.
  274.      * @exception OrekitException if some resource file reading error occurs
  275.      */
  276.     private static void readcoefficients() throws OrekitException {

  278.         final int size = NLATM + 1;
  279.         tt   = new double[size];
  280.         h    = new double[size];
  281.         he   = new double[size];
  282.         o    = new double[size];
  283.         az2  = new double[size];
  284.         o2   = new double[size];
  285.         az   = new double[size];
  286.         t0   = new double[size];
  287.         tp   = new double[size];

  289.         final InputStream in = DTM2000.class.getResourceAsStream(DTM2000);
  290.         if (in == null) {
  291.             throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_RESOURCE, DTM2000);
  292.         }

  294.         BufferedReader r = null;
  295.         try {

  297.             r = new BufferedReader(new InputStreamReader(in, "UTF-8"));
  298.             r.readLine();
  299.             r.readLine();
  300.             for (String line = r.readLine(); line != null; line = r.readLine()) {
  301.                 final int num = Integer.parseInt(line.substring(0, 4).replace(' ', '0'));
  302.                 line = line.substring(4);
  303.                 tt[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  304.                 line = line.substring(13 + 9);
  305.                 h[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  306.                 line = line.substring(13 + 9);
  307.                 he[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  308.                 line = line.substring(13 + 9);
  309.                 o[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  310.                 line = line.substring(13 + 9);
  311.                 az2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  312.                 line = line.substring(13 + 9);
  313.                 o2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  314.                 line = line.substring(13 + 9);
  315.                 az[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  316.                 line = line.substring(13 + 9);
  317.                 t0[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  318.                 line = line.substring(13 + 9);
  319.                 tp[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  320.             }
  321.         } catch (IOException ioe) {
  322.             throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
  323.         } finally {
  324.             if (r != null) {
  325.                 try {
  326.                     r.close();
  327.                 } catch (IOException ioe) {
  328.                     throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
  329.                 }
  330.             }
  331.         }
  332.     }

  334.     /** Get the local density.
  335.      * @param date current date
  336.      * @param position current position in frame
  337.      * @param frame the frame in which is defined the position
  338.      * @return local density (kg/m³)
  339.      * @exception OrekitException if date is out of range of solar activity model
  340.      * or if some frame conversion cannot be performed
  341.      */
  342.     public double getDensity(final AbsoluteDate date, final Vector3D position,
  343.                              final Frame frame)
  344.         throws OrekitException {

  346.         // check if data are available :
  347.         if ((date.compareTo(inputParams.getMaxDate()) > 0) ||
  348.             (date.compareTo(inputParams.getMinDate()) < 0)) {
  349.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  350.                                       date, inputParams.getMinDate(), inputParams.getMaxDate());
  351.         }

  353.         // compute day number in current year
  354.         final Calendar cal = new GregorianCalendar();
  355.         cal.setTime(date.toDate(TimeScalesFactory.getUTC()));
  356.         final int day = cal.get(Calendar.DAY_OF_YEAR);
  357.         //position in ECEF so we only have to do the transform once
  358.         final Frame ecef = earth.getBodyFrame();
  359.         final Vector3D pEcef = frame.getTransformTo(ecef, date)
  360.                 .transformPosition(position);
  361.         // compute geodetic position
  362.         final GeodeticPoint inBody = earth.transform(pEcef, ecef, date);
  363.         final double alti = inBody.getAltitude();
  364.         final double lon = inBody.getLongitude();
  365.         final double lat = inBody.getLatitude();

  367.         // compute local solar time
  368.         final Vector3D sunPos = sun.getPVCoordinates(date, ecef).getPosition();
  369.         final double hl = FastMath.PI + FastMath.atan2(
  370.                 sunPos.getX() * pEcef.getY() - sunPos.getY() * pEcef.getX(),
  371.                 sunPos.getX() * pEcef.getX() + sunPos.getY() * pEcef.getY());

  373.         // get current solar activity data and compute
  374.         return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(date),
  375.                           inputParams.getMeanFlux(date), inputParams.getThreeHourlyKP(date),
  376.                           inputParams.get24HoursKp(date));

  378.     }

  380.     /** {@inheritDoc} */
  381.     @Override
  382.     public <T extends RealFieldElement<T>> T
  383.         getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position,
  384.                    final Frame frame)
  385.             throws OrekitException {
  386.         // check if data are available :
  387.         final AbsoluteDate dateD = date.toAbsoluteDate();
  388.         if ((dateD.compareTo(inputParams.getMaxDate()) > 0) ||
  389.             (dateD.compareTo(inputParams.getMinDate()) < 0)) {
  390.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  391.                                       dateD, inputParams.getMinDate(), inputParams.getMaxDate());
  392.         }

  394.         // compute day number in current year
  395.         final Calendar cal = new GregorianCalendar();
  396.         cal.setTime(date.toDate(TimeScalesFactory.getUTC()));
  397.         final int day = cal.get(Calendar.DAY_OF_YEAR);
  398.         //position in ECEF so we only have to do the transform once
  399.         final Frame ecef = earth.getBodyFrame();
  400.         final FieldVector3D<T> pEcef = frame.getTransformTo(ecef, date).transformPosition(position);
  401.         // compute geodetic position
  402.         final FieldGeodeticPoint<T> inBody = earth.transform(pEcef, ecef, date);
  403.         final T alti = inBody.getAltitude();
  404.         final T lon = inBody.getLongitude();
  405.         final T lat = inBody.getLatitude();

  407.         // compute local solar time
  408.         final Vector3D sunPos = sun.getPVCoordinates(dateD, ecef).getPosition();
  409.         final T y  = pEcef.getY().multiply(sunPos.getX()).subtract(pEcef.getX().multiply(sunPos.getY()));
  410.         final T x  = pEcef.getX().multiply(sunPos.getX()).add(pEcef.getY().multiply(sunPos.getY()));
  411.         final T hl = y.atan2(x).add(FastMath.PI);

  413.         // get current solar activity data and compute
  414.         return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(dateD),
  415.                           inputParams.getMeanFlux(dateD), inputParams.getThreeHourlyKP(dateD),
  416.                           inputParams.get24HoursKp(dateD));
  417.     }

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

  422.         /** Number of days in current year. */
  423.         private final int day;

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

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

  431.         /** Kp coefficients.
  432.          * <ul>
  433.          *   <li>akp[1] = 3-hourly kp</li>
  434.          *   <li>akp[2] = 0 (not used)</li>
  435.          *   <li>akp[3] = mean kp of last 24 hours</li>
  436.          *   <li>akp[4] = 0 (not used)</li>
  437.          * </ul>
  438.          */
  439.         private final double[] akp;

  441.         /** Cosine of the longitude. */
  442.         private final double clfl;

  444.         /** Sine of the longitude. */
  445.         private final double slfl;

  447.         /** Total density (g/cm3). */
  448.         private final double ro;

  450.         // CHECKSTYLE: stop JavadocVariable check

  452.         /** Legendre coefficients. */
  453.         private final double p10;
  454.         private final double p20;
  455.         private final double p30;
  456.         private final double p40;
  457.         private final double p50;
  458.         private final double p60;
  459.         private final double p11;
  460.         private final double p21;
  461.         private final double p31;
  462.         private final double p41;
  463.         private final double p51;
  464.         private final double p22;
  465.         private final double p32;
  466.         private final double p42;
  467.         private final double p52;
  468.         private final double p62;
  469.         private final double p33;
  470.         private final double p10mg;
  471.         private final double p20mg;
  472.         private final double p40mg;

  474.         /** Local time intermediate values. */
  475.         private final double hl0;
  476.         private final double ch;
  477.         private final double sh;
  478.         private final double c2h;
  479.         private final double s2h;
  480.         private final double c3h;
  481.         private final double s3h;

  483.         /** Simple constructor.
  484.          * @param day day of year
  485.          * @param altiKM altitude <em>in kilometers</em>
  486.          * @param lon local longitude (rad)
  487.          * @param lat local latitude (rad)
  488.          * @param hl local solar time in rad (O hr = 0 rad)
  489.          * @param f instantaneous solar flux (F10.7)
  490.          * @param fbar mean solar flux (F10.7)
  491.          * @param akp geomagnetic activity index
  492.          */
  493.         Computation(final int day,
  494.                     final double altiKM, final double lon, final double lat,
  495.                     final double hl, final double[] f, final double[] fbar,
  496.                     final double[] akp) {

  498.             this.day  = day;
  499.             this.f    = f;
  500.             this.fbar = fbar;
  501.             this.akp  = akp;

  503.             // compute Legendre polynomials wrt geographic pole
  504.             final double c = FastMath.sin(lat);
  505.             final double c2 = c * c;
  506.             final double c4 = c2 * c2;
  507.             final double s = FastMath.cos(lat);
  508.             final double s2 = s * s;
  509.             p10 = c;
  510.             p20 = 1.5 * c2 - 0.5;
  511.             p30 = c * (2.5 * c2 - 1.5);
  512.             p40 = 4.375 * c4 - 3.75 * c2 + 0.375;
  513.             p50 = c * (7.875 * c4 - 8.75 * c2 + 1.875);
  514.             p60 = (5.5 * c * p50 - 2.5 * p40) / 3.0;
  515.             p11 = s;
  516.             p21 = 3.0 * c * s;
  517.             p31 = s * (7.5 * c2 - 1.5);
  518.             p41 = c * s * (17.5 * c2 - 7.5);
  519.             p51 = s * (39.375 * c4 - 26.25 * c2 + 1.875);
  520.             p22 = 3.0 * s2;
  521.             p32 = 15.0 * c * s2;
  522.             p42 = s2 * (52.5 * c2 - 7.5);
  523.             p52 = 3.0 * c * p42 - 2.0 * p32;
  524.             p62 = 2.75 * c * p52 - 1.75 * p42;
  525.             p33 = 15.0 * s * s2;

  527.             // compute Legendre polynomials wrt magnetic pole (79N, 71W)
  528.             final double clmlmg = FastMath.cos(lon - XLMG);
  529.             final double cmg  = s * CPMG * clmlmg + c * SPMG;
  530.             final double cmg2 = cmg * cmg;
  531.             final double cmg4 = cmg2 * cmg2;
  532.             p10mg = cmg;
  533.             p20mg = 1.5 * cmg2 - 0.5;
  534.             p40mg = 4.375 * cmg4 - 3.75 * cmg2 + 0.375;

  536.             clfl = FastMath.cos(lon);
  537.             slfl = FastMath.sin(lon);

  539.             // local time
  540.             hl0 = hl;
  541.             ch  = FastMath.cos(hl0);
  542.             sh  = FastMath.sin(hl0);
  543.             c2h = ch * ch - sh * sh;
  544.             s2h = 2.0 * ch * sh;
  545.             c3h = c2h * ch - s2h * sh;
  546.             s3h = s2h * ch + c2h * sh;

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

  550.             final double[] dtt  = new double[tt.length];
  551.             final double[] dh   = new double[tt.length];
  552.             final double[] dhe  = new double[tt.length];
  553.             final double[] dox  = new double[tt.length];
  554.             final double[] daz2 = new double[tt.length];
  555.             final double[] do2  = new double[tt.length];
  556.             final double[] daz  = new double[tt.length];
  557.             final double[] dt0  = new double[tt.length];
  558.             final double[] dtp  = new double[tt.length];

  560.             Arrays.fill(dtt,  Double.NaN);
  561.             Arrays.fill(dh,   Double.NaN);
  562.             Arrays.fill(dhe,  Double.NaN);
  563.             Arrays.fill(dox,  Double.NaN);
  564.             Arrays.fill(daz2, Double.NaN);
  565.             Arrays.fill(do2,  Double.NaN);
  566.             Arrays.fill(daz,  Double.NaN);
  567.             Arrays.fill(dt0,  Double.NaN);
  568.             Arrays.fill(dtp,  Double.NaN);

  570.             //  compute function g(l) / tinf, t120, tp120
  571.             int kleq = 1;
  572.             final double gdelt = gFunction(tt, dtt, 1, kleq);
  573.             dtt[1] = 1.0 + gdelt;
  574.             final double tinf   = tt[1] * dtt[1];

  576.             kleq = 0; // equinox

  578.             if ((day < 59) || (day > 284)) {
  579.                 kleq = -1; // north winter
  580.             }
  581.             if ((day > 99) && (day < 244)) {
  582.                 kleq = 1; // north summer
  583.             }

  585.             final double gdelt0 =  gFunction(t0, dt0, 0, kleq);
  586.             dt0[1] = (t0[1] + gdelt0) / t0[1];
  587.             final double t120 = t0[1] + gdelt0;
  588.             final double gdeltp = gFunction(tp, dtp, 0, kleq);
  589.             dtp[1] = (tp[1] + gdeltp) / tp[1];
  590.             final double tp120 = tp[1] + gdeltp;

  592.             // compute n(z) concentrations: H, He, O, N2, O2, N
  593.             final double sigma   = tp120 / (tinf - t120);
  594.             final double dzeta   = (RE + zlb) / (RE + altiKM);
  595.             final double zeta    = (altiKM - zlb) * dzeta;
  596.             final double sigzeta = sigma * zeta;
  597.             final double expsz   = FastMath.exp(-sigzeta);
  598.             final double tz = tinf - (tinf - t120) * expsz;

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

  602.             kleq = 1;

  604.             final double gdelh = gFunction(h, dh, 0, kleq);
  605.             dh[1] = FastMath.exp(gdelh);
  606.             dbase[1] = h[1] * dh[1];

  608.             final double gdelhe = gFunction(he, dhe, 0, kleq);
  609.             dhe[1] = FastMath.exp(gdelhe);
  610.             dbase[2] = he[1] * dhe[1];

  612.             final double gdelo = gFunction(o, dox, 1, kleq);
  613.             dox[1] = FastMath.exp(gdelo);
  614.             dbase[3] = o[1] * dox[1];

  616.             final double gdelaz2 = gFunction(az2, daz2, 1, kleq);
  617.             daz2[1] = FastMath.exp(gdelaz2);
  618.             dbase[4] = az2[1] * daz2[1];

  620.             final double gdelo2 = gFunction(o2, do2, 1, kleq);
  621.             do2[1] = FastMath.exp(gdelo2);
  622.             dbase[5] = o2[1] * do2[1];

  624.             final double gdelaz = gFunction(az, daz, 1, kleq);
  625.             daz[1] = FastMath.exp(gdelaz);
  626.             dbase[6] = az[1] * daz[1];

  628.             final double zlbre  = 1.0 + zlb / RE;
  629.             final double glb    = (GSURF / (zlbre * zlbre)) / (sigma * RGAS * tinf);
  630.             final double t120tz = t120 / tz;

  632.             // Partial densities in (g/cm3).
  633.             // d(1) = hydrogen
  634.             // d(2) = helium
  635.             // d(3) = atomic oxygen
  636.             // d(4) = molecular nitrogen
  637.             // d(5) = molecular oxygen
  638.             // d(6) = atomic nitrogen
  639.             double tmpro = 0;
  640.             for (int i = 1; i <= 6; i++) {
  641.                 final double gamma = MA[i] * glb;
  642.                 final double upapg = 1.0 + ALEFA[i] + gamma;
  643.                 final double fzI = FastMath.pow(t120tz, upapg) * FastMath.exp(-sigzeta * gamma);
  644.                 // concentrations of H, He, O, N2, O2, N (particles/cm³)
  645.                 final double ccI = dbase[i] * fzI;
  646.                 // contribution of densities of H, He, O, N2, O2, N (g/cm³)
  647.                 tmpro += ccI * VMA[i];
  648.             }
  649.             this.ro = tmpro;

  651.         }

  653.         /** Computation of function G.
  654.          * @param a vector of coefficients for computation
  655.          * @param da vector of partial derivatives
  656.          * @param ff0 coefficient flag (1 for Ox, Az, He, T°; 0 for H and tp120)
  657.          * @param kle_eq season indicator flag (summer, winter, equinox)
  658.          * @return value of G
  659.          */
  660.         private double gFunction(final double[] a, final double[] da,
  661.                                  final int ff0, final int kle_eq) {

  663.             final double[] fmfb   = new double[3];
  664.             final double[] fbm150 = new double[3];

  666.             // latitude terms
  667.             da[2]  = p20;
  668.             da[3]  = p40;
  669.             da[74] = p10;
  670.             double a74 = a[74];
  671.             double a77 = a[77];
  672.             double a78 = a[78];
  673.             if (kle_eq == -1) {
  674.                 // winter
  675.                 a74 = -a74;
  676.                 a77 = -a77;
  677.                 a78 = -a78;
  678.             }
  679.             if (kle_eq == 0 ) {
  680.                 // equinox
  681.                 a74 = semestrialCorrection(a74);
  682.                 a77 = semestrialCorrection(a77);
  683.                 a78 = semestrialCorrection(a78);
  684.             }
  685.             da[77] = p30;
  686.             da[78] = p50;
  687.             da[79] = p60;

  689.             // flux terms
  690.             fmfb[1]   = f[1] - fbar[1];
  691.             fmfb[2]   = f[2] - fbar[2];
  692.             fbm150[1] = fbar[1] - 150.0;
  693.             fbm150[2] = fbar[2];
  694.             da[4]     = fmfb[1];
  695.             da[6]     = fbm150[1];
  696.             da[4]     = da[4] + a[70] * fmfb[2];
  697.             da[6]     = da[6] + a[71] * fbm150[2];
  698.             da[70]    = fmfb[2] * (a[4] + 2.0 * a[5] * da[4] + a[82] * p10 +
  699.                                    a[83] * p20 + a[84] * p30);
  700.             da[71]    = fbm150[2] * (a[6] + 2.0 * a[69] * da[6] + a[85] * p10 +
  701.                                      a[86] * p20 + a[87] * p30);
  702.             da[5]     = da[4] * da[4];
  703.             da[69]    = da[6] * da[6];
  704.             da[82]    = da[4] * p10;
  705.             da[83]    = da[4] * p20;
  706.             da[84]    = da[4] * p30;
  707.             da[85]    = da[6] * p20;
  708.             da[86]    = da[6] * p30;
  709.             da[87]    = da[6] * p40;

  711.             // Kp terms
  712.             final int ikp  = 62;
  713.             final int ikpm = 67;
  714.             final double c2fi = 1.0 - p10mg * p10mg;
  715.             final double dkp  = akp[1] + (a[ikp] + c2fi * a[ikp + 1]) * akp[2];
  716.             double dakp = a[7] + a[8] * p20mg + a[68] * p40mg +
  717.                           2.0 * dkp * (a[60] + a[61] * p20mg +
  718.                                        a[75] * 2.0 * dkp * dkp);
  719.             da[ikp] = dakp * akp[2];
  720.             da[ikp + 1] = da[ikp] * c2fi;
  721.             final double dkpm  = akp[3] + a[ikpm] * akp[4];
  722.             final double dakpm = a[64] + a[65] * p20mg + a[72] * p40mg +
  723.                                  2.0 * dkpm * (a[66] + a[73] * p20mg +
  724.                                                a[76] * 2.0 * dkpm * dkpm);
  725.             da[ikpm] = dakpm * akp[4];
  726.             da[7]    = dkp;
  727.             da[8]    = p20mg * dkp;
  728.             da[68]   = p40mg * dkp;
  729.             da[60]   = dkp * dkp;
  730.             da[61]   = p20mg * da[60];
  731.             da[75]   = da[60] * da[60];
  732.             da[64]   = dkpm;
  733.             da[65]   = p20mg * dkpm;
  734.             da[72]   = p40mg * dkpm;
  735.             da[66]   = dkpm * dkpm;
  736.             da[73]   = p20mg * da[66];
  737.             da[76]   = da[66] * da[66];

  739.             // non-periodic g(l) function
  740.             double f0 = a[4]  * da[4]  + a[5]  * da[5]  + a[6]  * da[6]  +
  741.                         a[69] * da[69] + a[82] * da[82] + a[83] * da[83] +
  742.                         a[84] * da[84] + a[85] * da[85] + a[86] * da[86] +
  743.                         a[87] * da[87];
  744.             final double f1f = 1.0 + f0 * ff0;

  746.             f0 = f0 + a[2] * da[2] + a[3] * da[3] + a74 * da[74] +
  747.                  a77 * da[77] + a[7] * da[7] + a[8] * da[8] +
  748.                  a[60] * da[60] + a[61] * da[61] + a[68] * da[68] +
  749.                  a[64] * da[64] + a[65] * da[65] + a[66] * da[66] +
  750.                  a[72] * da[72] + a[73] * da[73] + a[75] * da[75] +
  751.                  a[76] * da[76] + a78   * da[78] + a[79] * da[79];
  752. //          termes annuels symetriques en latitude
  753.             da[9]  = FastMath.cos(ROT * (day - a[11]));
  754.             da[10] = p20 * da[9];
  755. //          termes semi-annuels symetriques en latitude
  756.             da[12] = FastMath.cos(ROT2 * (day - a[14]));
  757.             da[13] = p20 * da[12];
  758. //          termes annuels non symetriques en latitude
  759.             final double coste = FastMath.cos(ROT * (day - a[18]));
  760.             da[15] = p10 * coste;
  761.             da[16] = p30 * coste;
  762.             da[17] = p50 * coste;
  763. //          terme  semi-annuel  non symetrique  en latitude
  764.             final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
  765.             da[19] = p10 * cos2te;
  766.             da[39] = p30 * cos2te;
  767.             da[59] = p50 * cos2te;
  768. //          termes diurnes [et couples annuel]
  769.             da[21] = p11 * ch;
  770.             da[22] = p31 * ch;
  771.             da[23] = p51 * ch;
  772.             da[24] = da[21] * coste;
  773.             da[25] = p21 * ch * coste;
  774.             da[26] = p11 * sh;
  775.             da[27] = p31 * sh;
  776.             da[28] = p51 * sh;
  777.             da[29] = da[26] * coste;
  778.             da[30] = p21 * sh * coste;
  779. //          termes semi-diurnes [et couples annuel]
  780.             da[31] = p22 * c2h;
  781.             da[37] = p42 * c2h;
  782.             da[32] = p32 * c2h * coste;
  783.             da[33] = p22 * s2h;
  784.             da[38] = p42 * s2h;
  785.             da[34] = p32 * s2h * coste;
  786.             da[88] = p32 * c2h;
  787.             da[89] = p32 * s2h;
  788.             da[90] = p52 * c2h;
  789.             da[91] = p52 * s2h;
  790.             double a88 = a[88];
  791.             double a89 = a[89];
  792.             double a90 = a[90];
  793.             double a91 = a[91];
  794.             if (kle_eq == -1) {            //hiver
  795.                 a88 = -a88;
  796.                 a89 = -a89;
  797.                 a90 = -a90;
  798.                 a91 = -a91;
  799.             }
  800.             if (kle_eq == 0) {             //equinox
  801.                 a88 = semestrialCorrection(a88);
  802.                 a89 = semestrialCorrection(a89);
  803.                 a90 = semestrialCorrection(a90);
  804.                 a91 = semestrialCorrection(a91);
  805.             }
  806.             da[92] = p62 * c2h;
  807.             da[93] = p62 * s2h;
  808. //          termes ter-diurnes
  809.             da[35] = p33 * c3h;
  810.             da[36] = p33 * s3h;
  811. //          fonction g[l] periodique
  812.             double fp = a[9]  * da[9]  + a[10] * da[10] + a[12] * da[12] + a[13] * da[13] +
  813.                         a[15] * da[15] + a[16] * da[16] + a[17] * da[17] + a[19] * da[19] +
  814.                         a[21] * da[21] + a[22] * da[22] + a[23] * da[23] + a[24] * da[24] +
  815.                         a[25] * da[25] + a[26] * da[26] + a[27] * da[27] + a[28] * da[28] +
  816.                         a[29] * da[29] + a[30] * da[30] + a[31] * da[31] + a[32] * da[32] +
  817.                         a[33] * da[33] + a[34] * da[34] + a[35] * da[35] + a[36] * da[36] +
  818.                         a[37] * da[37] + a[38] * da[38] + a[39] * da[39] + a[59] * da[59] +
  819.                         a88   * da[88] + a89   * da[89] + a90   * da[90] + a91   * da[91] +
  820.                         a[92] * da[92] + a[93] * da[93];
  821. //          termes d'activite magnetique
  822.             da[40] = p10 * coste * dkp;
  823.             da[41] = p30 * coste * dkp;
  824.             da[42] = p50 * coste * dkp;
  825.             da[43] = p11 * ch * dkp;
  826.             da[44] = p31 * ch * dkp;
  827.             da[45] = p51 * ch * dkp;
  828.             da[46] = p11 * sh * dkp;
  829.             da[47] = p31 * sh * dkp;
  830.             da[48] = p51 * sh * dkp;

  832. //          fonction g[l] periodique supplementaire
  833.             fp += a[40] * da[40] + a[41] * da[41] + a[42] * da[42] + a[43] * da[43] +
  834.                   a[44] * da[44] + a[45] * da[45] + a[46] * da[46] + a[47] * da[47] +
  835.                   a[48] * da[48];

  837.             dakp = (a[40] * p10 + a[41] * p30 + a[42] * p50) * coste +
  838.                    (a[43] * p11 + a[44] * p31 + a[45] * p51) * ch +
  839.                    (a[46] * p11 + a[47] * p31 + a[48] * p51) * sh;
  840.             da[ikp] += dakp * akp[2];
  841.             da[ikp + 1] = da[ikp] + dakp * c2fi * akp[2];
  842. //          termes de longitude
  843.             da[49] = p11 * clfl;
  844.             da[50] = p21 * clfl;
  845.             da[51] = p31 * clfl;
  846.             da[52] = p41 * clfl;
  847.             da[53] = p51 * clfl;
  848.             da[54] = p11 * slfl;
  849.             da[55] = p21 * slfl;
  850.             da[56] = p31 * slfl;
  851.             da[57] = p41 * slfl;
  852.             da[58] = p51 * slfl;

  854. //          fonction g[l] periodique supplementaire
  855.             fp += a[49] * da[49] + a[50] * da[50] + a[51] * da[51] + a[52] * da[52] +
  856.                   a[53] * da[53] + a[54] * da[54] + a[55] * da[55] + a[56] * da[56] +
  857.                   a[57] * da[57] + a[58] * da[58];

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

  862.         }


  865.         /** Apply a correction coefficient to the given parameter.
  866.          * @param param the parameter to correct
  867.          * @return the corrected parameter
  868.          */
  869.         private double semestrialCorrection(final double param) {
  870.             final int debeq_pr = 59;
  871.             final int debeq_au = 244;
  872.             final double result;
  873.             if (day >= 100) {
  874.                 final double xmult  = (day - debeq_au) / 40.0;
  875.                 result = param - 2.0 * param * xmult;
  876.             } else {
  877.                 final double xmult  = (day - debeq_pr) / 40.0;
  878.                 result = 2.0 * param * xmult - param;
  879.             }
  880.             return result;
  881.         }


  884.     }

  886.     /** Local holder for intermediate results ensuring the model is reentrant.
  887.      * @param <T> type of the field elements
  888.      */
  889.     private static class FieldComputation<T extends RealFieldElement<T>> {

  891.         /** Number of days in current year. */
  892.         private int day;

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

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

  900.         /** Kp coefficients.
  901.          * <ul>
  902.          *   <li>akp[1] = 3-hourly kp</li>
  903.          *   <li>akp[2] = 0 (not used)</li>
  904.          *   <li>akp[3] = mean kp of last 24 hours</li>
  905.          *   <li>akp[4] = 0 (not used)</li>
  906.          * </ul>
  907.          */
  908.         private double[] akp = new double[5];

  910.         /** Cosine of the longitude. */
  911.         private final T clfl;

  913.         /** Sine of the longitude. */
  914.         private final T slfl;

  916.         /** Total density (g/cm3). */
  917.         private final T ro;

  919.         // CHECKSTYLE: stop JavadocVariable check

  921.         /** Legendre coefficients. */
  922.         private final T p10;
  923.         private final T p20;
  924.         private final T p30;
  925.         private final T p40;
  926.         private final T p50;
  927.         private final T p60;
  928.         private final T p11;
  929.         private final T p21;
  930.         private final T p31;
  931.         private final T p41;
  932.         private final T p51;
  933.         private final T p22;
  934.         private final T p32;
  935.         private final T p42;
  936.         private final T p52;
  937.         private final T p62;
  938.         private final T p33;
  939.         private final T p10mg;
  940.         private final T p20mg;
  941.         private final T p40mg;

  943.         /** Local time intermediate values. */
  944.         private final T hl0;
  945.         private final T ch;
  946.         private final T sh;
  947.         private final T c2h;
  948.         private final T s2h;
  949.         private final T c3h;
  950.         private final T s3h;

  952.         /** Simple constructor.
  953.          * @param day day of year
  954.          * @param altiKM altitude <em>in kilometers</em>
  955.          * @param lon local longitude (rad)
  956.          * @param lat local latitude (rad)
  957.          * @param hl local solar time in rad (O hr = 0 rad)
  958.          * @param f instantaneous solar flux (F10.7)
  959.          * @param far mean solar flux (F10.7)
  960.          * @param akp geomagnetic activity index
  961.          */
  962.         FieldComputation(final int day,
  963.                          final T altiKM, final T lon, final T lat,
  964.                          final T hl, final double[] f, final double[] far,
  965.                          final double[] akp) {

  967.             this.day  = day;
  968.             this.f    = f;
  969.             this.fbar = far;
  970.             this.akp  = akp;

  972.             // compute Legendre polynomials wrt geographic pole
  973.             final T c = lat.sin();
  974.             final T c2 = c.multiply(c);
  975.             final T c4 = c2.multiply(c2);
  976.             final T s = lat.cos();
  977.             final T s2 = s.multiply(s);
  978.             p10 = c;
  979.             p20 = c2.multiply(1.5).subtract(0.5);
  980.             p30 = c.multiply(c2.multiply(2.5).subtract(1.5));
  981.             p40 = c4.multiply(4.375).subtract(c2.multiply(3.75)).add(0.375);
  982.             p50 = c.multiply(c4.multiply(7.875).subtract(c2.multiply(8.75)).add(1.875));
  983.             p60 = (c.multiply(5.5).multiply(p50).subtract(p40.multiply(2.5))).divide(3.0);
  984.             p11 = s;
  985.             p21 = c.multiply(3.0).multiply(s);
  986.             p31 = s.multiply(c2.multiply(7.5).subtract(1.5));
  987.             p41 = c.multiply(s).multiply(c2.multiply(17.5).subtract(7.5));
  988.             p51 = s.multiply(c4.multiply(39.375).subtract(c2.multiply(26.25)).add(1.875));
  989.             p22 = s2.multiply(3.0);
  990.             p32 = c.multiply(15.0).multiply(s2);
  991.             p42 = s2.multiply(c2.multiply(52.5).subtract(7.5));
  992.             p52 = c.multiply(3.0).multiply(p42).subtract(p32.multiply(2.0));
  993.             p62 = c.multiply(2.75).multiply(p52).subtract(p42.multiply(1.75));
  994.             p33 = s.multiply(15.0).multiply(s2);

  996.             // compute Legendre polynomials wrt magnetic pole (79N, 71W)
  997.             final T clmlmg = lon.subtract(XLMG).cos();
  998.             final T cmg  = s.multiply(CPMG).multiply(clmlmg).add(c.multiply(SPMG));
  999.             final T cmg2 = cmg.multiply(cmg);
  1000.             final T cmg4 = cmg2.multiply(cmg2);
  1001.             p10mg = cmg;
  1002.             p20mg = cmg2.multiply(1.5).subtract(0.5);
  1003.             p40mg = cmg4.multiply(4.375).subtract(cmg2.multiply(3.75)).add(0.375);

  1005.             clfl = lon.cos();
  1006.             slfl = lon.sin();

  1008.             // local time
  1009.             hl0 = hl;
  1010.             ch  = hl0.cos();
  1011.             sh  = hl0.sin();
  1012.             c2h = ch.multiply(ch).subtract(sh.multiply(sh));
  1013.             s2h = ch.multiply(sh).multiply(2);
  1014.             c3h = c2h.multiply(ch).subtract(s2h.multiply(sh));
  1015.             s3h = s2h.multiply(ch).add(c2h.multiply(sh));

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

  1019.             final T[] dtt  = MathArrays.buildArray(altiKM.getField(), tt.length);
  1020.             final T[] dh   = MathArrays.buildArray(altiKM.getField(), tt.length);
  1021.             final T[] dhe  = MathArrays.buildArray(altiKM.getField(), tt.length);
  1022.             final T[] dox  = MathArrays.buildArray(altiKM.getField(), tt.length);
  1023.             final T[] daz2 = MathArrays.buildArray(altiKM.getField(), tt.length);
  1024.             final T[] do2  = MathArrays.buildArray(altiKM.getField(), tt.length);
  1025.             final T[] daz  = MathArrays.buildArray(altiKM.getField(), tt.length);
  1026.             final T[] dt0  = MathArrays.buildArray(altiKM.getField(), tt.length);
  1027.             final T[] dtp  = MathArrays.buildArray(altiKM.getField(), tt.length);

  1029.             //  compute function g(l) / tinf, t120, tp120
  1030.             int kleq = 1;
  1031.             final T gdelt = gFunction(tt, dtt, 1, kleq);
  1032.             dtt[1] = gdelt.add(1);
  1033.             final T tinf   = dtt[1].multiply(tt[1]);

  1035.             kleq = 0; // equinox

  1037.             if ((day < 59) || (day > 284)) {
  1038.                 kleq = -1; // north winter
  1039.             }
  1040.             if ((day > 99) && (day < 244)) {
  1041.                 kleq = 1; // north summer
  1042.             }

  1044.             final T gdelt0 =  gFunction(t0, dt0, 0, kleq);
  1045.             dt0[1] = gdelt0.add(t0[1]).divide(t0[1]);
  1046.             final T t120 = gdelt0.add(t0[1]);
  1047.             final T gdeltp = gFunction(tp, dtp, 0, kleq);
  1048.             dtp[1] = gdeltp.add(tp[1]).divide(tp[1]);
  1049.             final T tp120 = gdeltp.add(tp[1]);

  1051.             // compute n(z) concentrations: H, He, O, N2, O2, N
  1052.             final T sigma   = tp120.divide(tinf.subtract(t120));
  1053.             final T dzeta   = altiKM.add(RE).reciprocal().multiply(zlb + RE);
  1054.             final T zeta    = altiKM.subtract(zlb).multiply(dzeta);
  1055.             final T sigzeta = sigma.multiply(zeta);
  1056.             final T expsz   = sigzeta.negate().exp();
  1057.             final T tz = tinf.subtract(tinf.subtract(t120).multiply(expsz));

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

  1061.             kleq = 1;

  1063.             final T gdelh = gFunction(h, dh, 0, kleq);
  1064.             dh[1] = gdelh.exp();
  1065.             dbase[1] = dh[1].multiply(h[1]);

  1067.             final T gdelhe = gFunction(he, dhe, 0, kleq);
  1068.             dhe[1] = gdelhe.exp();
  1069.             dbase[2] = dhe[1].multiply(he[1]);

  1071.             final T gdelo = gFunction(o, dox, 1, kleq);
  1072.             dox[1] = gdelo.exp();
  1073.             dbase[3] = dox[1].multiply(o[1]);

  1075.             final T gdelaz2 = gFunction(az2, daz2, 1, kleq);
  1076.             daz2[1] = gdelaz2.exp();
  1077.             dbase[4] = daz2[1].multiply(az2[1]);

  1079.             final T gdelo2 = gFunction(o2, do2, 1, kleq);
  1080.             do2[1] = gdelo2.exp();
  1081.             dbase[5] = do2[1].multiply(o2[1]);

  1083.             final T gdelaz = gFunction(az, daz, 1, kleq);
  1084.             daz[1] = gdelaz.exp();
  1085.             dbase[6] = daz[1].multiply(az[1]);

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

  1091.             // Partial densities in (g/cm3).
  1092.             // d(1) = hydrogen
  1093.             // d(2) = helium
  1094.             // d(3) = atomic oxygen
  1095.             // d(4) = molecular nitrogen
  1096.             // d(5) = molecular oxygen
  1097.             // d(6) = atomic nitrogen
  1098.             T tmpro = altiKM.getField().getZero();
  1099.             for (int i = 1; i <= 6; i++) {
  1100.                 final T gamma = glb.multiply(MA[i]);
  1101.                 final T upapg = gamma.add(1.0 + ALEFA[i]);
  1102.                 final T fzI = t120tz.pow(upapg).multiply(sigzeta.negate().multiply(gamma).exp());
  1103.                 // concentrations of H, He, O, N2, O2, N (particles/cm³)
  1104.                 final T ccI = dbase[i].multiply(fzI);
  1105.                 // contribution of densities of H, He, O, N2, O2, N (g/cm³)
  1106.                 tmpro = tmpro.add(ccI.multiply(VMA[i]));
  1107.             }
  1108.             this.ro = tmpro;

  1110.         }

  1112.         /** Computation of function G.
  1113.          * @param a vector of coefficients for computation
  1114.          * @param da vector of partial derivatives
  1115.          * @param ff0 coefficient flag (1 for Ox, Az, He, T°; 0 for H and tp120)
  1116.          * @param kle_eq season indicator flag (summer, winter, equinox)
  1117.          * @return value of G
  1118.          */
  1119.         private T gFunction(final double[] a, final T[] da,
  1120.                             final int ff0, final int kle_eq) {

  1122.             final T zero = da[0].getField().getZero();
  1123.             final double[] fmfb   = new double[3];
  1124.             final double[] fbm150 = new double[3];

  1126.             // latitude terms
  1127.             da[2]  = p20;
  1128.             da[3]  = p40;
  1129.             da[74] = p10;
  1130.             double a74 = a[74];
  1131.             double a77 = a[77];
  1132.             double a78 = a[78];
  1133.             if (kle_eq == -1) {
  1134.                 // winter
  1135.                 a74 = -a74;
  1136.                 a77 = -a77;
  1137.                 a78 = -a78;
  1138.             }
  1139.             if (kle_eq == 0 ) {
  1140.                 // equinox
  1141.                 a74 = semestrialCorrection(a74);
  1142.                 a77 = semestrialCorrection(a77);
  1143.                 a78 = semestrialCorrection(a78);
  1144.             }
  1145.             da[77] = p30;
  1146.             da[78] = p50;
  1147.             da[79] = p60;

  1149.             // flux terms
  1150.             fmfb[1]   = f[1] - fbar[1];
  1151.             fmfb[2]   = f[2] - fbar[2];
  1152.             fbm150[1] = fbar[1] - 150.0;
  1153.             fbm150[2] = fbar[2];
  1154.             da[4]     = zero.add(fmfb[1]);
  1155.             da[6]     = zero.add(fbm150[1]);
  1156.             da[4]     = da[4].add(a[70] * fmfb[2]);
  1157.             da[6]     = da[6].add(a[71] * fbm150[2]);
  1158.             da[70]    = da[4].multiply(a[ 5]).multiply(2).
  1159.                             add(p10.multiply(a[82])).
  1160.                             add(p20.multiply(a[83])).
  1161.                             add(p30.multiply(a[84])).
  1162.                             add(a[4]).
  1163.                         multiply(fmfb[2]);
  1164.             da[71]    = da[6].multiply(a[69]).multiply(2).
  1165.                             add(p10.multiply(a[85])).
  1166.                             add(p20.multiply(a[86])).
  1167.                             add(p30.multiply(a[87])).
  1168.                             add(a[6]).
  1169.                         multiply(fbm150[2]);
  1170.             da[5]     = da[4].multiply(da[4]);
  1171.             da[69]    = da[6].multiply(da[6]);
  1172.             da[82]    = da[4].multiply(p10);
  1173.             da[83]    = da[4].multiply(p20);
  1174.             da[84]    = da[4].multiply(p30);
  1175.             da[85]    = da[6].multiply(p20);
  1176.             da[86]    = da[6].multiply(p30);
  1177.             da[87]    = da[6].multiply(p40);

  1179.             // Kp terms
  1180.             final int ikp  = 62;
  1181.             final int ikpm = 67;
  1182.             final T c2fi = p10mg.multiply(p10mg).negate().add(1);
  1183.             final T dkp  = c2fi.multiply(a[ikp + 1]).add(a[ikp]).multiply(akp[2]).add(akp[1]);
  1184.             T dakp = p20mg.multiply(a[8]).add(p40mg.multiply(a[68])).
  1185.                      add(p20mg.multiply(a[61]).add(dkp.multiply(dkp).multiply(2 * a[75]).add(a[60])).multiply(dkp.multiply(2))).
  1186.                      add(a[7]);
  1187.             da[ikp]     = dakp.multiply(akp[2]);
  1188.             da[ikp + 1] = da[ikp].multiply(c2fi);
  1189.             final double dkpm  = akp[3] + a[ikpm] * akp[4];
  1190.             final T dakpm = p20mg.multiply(a[65]).add(p40mg.multiply(a[72])).
  1191.                             add(p20mg.multiply(a[73]).add(a[66] + a[76] * 2.0 * dkpm * dkpm).multiply( 2.0 * dkpm)).
  1192.                             add(a[64]);
  1193.             da[ikpm] = dakpm.multiply(akp[4]);
  1194.             da[7]    = dkp;
  1195.             da[8]    = p20mg.multiply(dkp);
  1196.             da[68]   = p40mg.multiply(dkp);
  1197.             da[60]   = dkp.multiply(dkp);
  1198.             da[61]   = p20mg.multiply(da[60]);
  1199.             da[75]   = da[60].multiply(da[60]);
  1200.             da[64]   = zero.add(dkpm);
  1201.             da[65]   = p20mg.multiply(dkpm);
  1202.             da[72]   = p40mg.multiply(dkpm);
  1203.             da[66]   = zero.add(dkpm * dkpm);
  1204.             da[73]   = p20mg.multiply(da[66]);
  1205.             da[76]   = da[66].multiply(da[66]);

  1207.             // non-periodic g(l) function
  1208.             T f0 = da[4].multiply(a[4]).
  1209.                    add(da[5].multiply(a[5])).
  1210.                    add(da[6].multiply(a[6])).
  1211.                    add(da[69].multiply(a[69])).
  1212.                    add(da[82].multiply(a[82])).
  1213.                    add(da[83].multiply(a[83])).
  1214.                    add(da[84].multiply(a[84])).
  1215.                    add(da[85].multiply(a[85])).
  1216.                    add(da[86].multiply(a[86])).
  1217.                    add(da[87].multiply(a[87]));
  1218.             final T f1f = f0.multiply(ff0).add(1);

  1220.             f0 = f0.
  1221.                  add(da[2].multiply(a[2])).
  1222.                  add(da[3].multiply(a[3])).
  1223.                  add(da[7].multiply(a[7])).
  1224.                  add(da[8].multiply(a[8])).
  1225.                  add(da[60].multiply(a[60])).
  1226.                  add(da[61].multiply(a[61])).
  1227.                  add(da[68].multiply(a[68])).
  1228.                  add(da[64].multiply(a[64])).
  1229.                  add(da[65].multiply(a[65])).
  1230.                  add(da[66].multiply(a[66])).
  1231.                  add(da[72].multiply(a[72])).
  1232.                  add(da[73].multiply(a[73])).
  1233.                  add(da[74].multiply(a74)).
  1234.                  add(da[75].multiply(a[75])).
  1235.                  add(da[76].multiply(a[76])).
  1236.                  add(da[77].multiply(a77)).
  1237.                  add(da[78].multiply(a78)).
  1238.                  add(da[79].multiply(a[79]));
  1239. //          termes annuels symetriques en latitude
  1240.             da[9]  = zero.add(FastMath.cos(ROT * (day - a[11])));
  1241.             da[10] = p20.multiply(da[9]);
  1242. //          termes semi-annuels symetriques en latitude
  1243.             da[12] = zero.add(FastMath.cos(ROT2 * (day - a[14])));
  1244.             da[13] = p20.multiply(da[12]);
  1245. //          termes annuels non symetriques en latitude
  1246.             final double coste = FastMath.cos(ROT * (day - a[18]));
  1247.             da[15] = p10.multiply(coste);
  1248.             da[16] = p30.multiply(coste);
  1249.             da[17] = p50.multiply(coste);
  1250. //          terme  semi-annuel  non symetrique  en latitude
  1251.             final double cos2te = FastMath.cos(ROT2 * (day - a[20]));
  1252.             da[19] = p10.multiply(cos2te);
  1253.             da[39] = p30.multiply(cos2te);
  1254.             da[59] = p50.multiply(cos2te);
  1255. //          termes diurnes [et couples annuel]
  1256.             da[21] = p11.multiply(ch);
  1257.             da[22] = p31.multiply(ch);
  1258.             da[23] = p51.multiply(ch);
  1259.             da[24] = da[21].multiply(coste);
  1260.             da[25] = p21.multiply(ch).multiply(coste);
  1261.             da[26] = p11.multiply(sh);
  1262.             da[27] = p31.multiply(sh);
  1263.             da[28] = p51.multiply(sh);
  1264.             da[29] = da[26].multiply(coste);
  1265.             da[30] = p21.multiply(sh).multiply(coste);
  1266. //          termes semi-diurnes [et couples annuel]
  1267.             da[31] = p22.multiply(c2h);
  1268.             da[37] = p42.multiply(c2h);
  1269.             da[32] = p32.multiply(c2h).multiply(coste);
  1270.             da[33] = p22.multiply(s2h);
  1271.             da[38] = p42.multiply(s2h);
  1272.             da[34] = p32.multiply(s2h).multiply(coste);
  1273.             da[88] = p32.multiply(c2h);
  1274.             da[89] = p32.multiply(s2h);
  1275.             da[90] = p52.multiply(c2h);
  1276.             da[91] = p52.multiply(s2h);
  1277.             double a88 = a[88];
  1278.             double a89 = a[89];
  1279.             double a90 = a[90];
  1280.             double a91 = a[91];
  1281.             if (kle_eq == -1) {            //hiver
  1282.                 a88 = -a88;
  1283.                 a89 = -a89;
  1284.                 a90 = -a90;
  1285.                 a91 = -a91;
  1286.             }
  1287.             if (kle_eq == 0) {             //equinox
  1288.                 a88 = semestrialCorrection(a88);
  1289.                 a89 = semestrialCorrection(a89);
  1290.                 a90 = semestrialCorrection(a90);
  1291.                 a91 = semestrialCorrection(a91);
  1292.             }
  1293.             da[92] = p62.multiply(c2h);
  1294.             da[93] = p62.multiply(s2h);
  1295. //          termes ter-diurnes
  1296.             da[35] = p33.multiply(c3h);
  1297.             da[36] = p33.multiply(s3h);
  1298. //          fonction g[l] periodique
  1299.             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])).
  1300.                    add(da[15].multiply(a[15])).add(da[16].multiply(a[16])).add(da[17].multiply(a[17])).add(da[19].multiply(a[19])).
  1301.                    add(da[21].multiply(a[21])).add(da[22].multiply(a[22])).add(da[23].multiply(a[23])).add(da[24].multiply(a[24])).
  1302.                    add(da[25].multiply(a[25])).add(da[26].multiply(a[26])).add(da[27].multiply(a[27])).add(da[28].multiply(a[28])).
  1303.                    add(da[29].multiply(a[29])).add(da[30].multiply(a[30])).add(da[31].multiply(a[31])).add(da[32].multiply(a[32])).
  1304.                    add(da[33].multiply(a[33])).add(da[34].multiply(a[34])).add(da[35].multiply(a[35])).add(da[36].multiply(a[36])).
  1305.                    add(da[37].multiply(a[37])).add(da[38].multiply(a[38])).add(da[39].multiply(a[39])).add(da[59].multiply(a[59])).
  1306.                    add(da[88].multiply(a88))  .add(da[89].multiply(a89)  ).add(da[90].multiply(a90)  ).add(da[91].multiply(a91)  ).
  1307.                    add(da[92].multiply(a[92])).add(da[93].multiply(a[93]));
  1308. //          termes d'activite magnetique
  1309.             da[40] = p10.multiply(coste).multiply(dkp);
  1310.             da[41] = p30.multiply(coste).multiply(dkp);
  1311.             da[42] = p50.multiply(coste).multiply(dkp);
  1312.             da[43] = p11.multiply(ch).multiply(dkp);
  1313.             da[44] = p31.multiply(ch).multiply(dkp);
  1314.             da[45] = p51.multiply(ch).multiply(dkp);
  1315.             da[46] = p11.multiply(sh).multiply(dkp);
  1316.             da[47] = p31.multiply(sh).multiply(dkp);
  1317.             da[48] = p51.multiply(sh).multiply(dkp);

  1319. //          fonction g[l] periodique supplementaire
  1320.             fp = fp.
  1321.                   add(da[40].multiply(a[40])).
  1322.                   add(da[41].multiply(a[41])).
  1323.                   add(da[42].multiply(a[42])).
  1324.                   add(da[43].multiply(a[43])).
  1325.                   add(da[44].multiply(a[44])).
  1326.                   add(da[45].multiply(a[45])).
  1327.                   add(da[46].multiply(a[46])).
  1328.                   add(da[47].multiply(a[47])).
  1329.                   add(da[48].multiply(a[48]));

  1331.             dakp =     p10.multiply(a[40]).add(p30.multiply(a[41])).add(p50.multiply(a[42])).multiply(coste).
  1332.                    add(p11.multiply(a[40]).add(p31.multiply(a[44])).add(p51.multiply(a[45])).multiply(ch)).
  1333.                    add(p11.multiply(a[46]).add(p31.multiply(a[47])).add(p51.multiply(a[48])).multiply(sh));
  1334.             da[ikp] = da[ikp].add(dakp.multiply(akp[2]));
  1335.             da[ikp + 1] = da[ikp].add(dakp.multiply(c2fi).multiply(akp[2]));
  1336. //          termes de longitude
  1337.             da[49] = p11.multiply(clfl);
  1338.             da[50] = p21.multiply(clfl);
  1339.             da[51] = p31.multiply(clfl);
  1340.             da[52] = p41.multiply(clfl);
  1341.             da[53] = p51.multiply(clfl);
  1342.             da[54] = p11.multiply(slfl);
  1343.             da[55] = p21.multiply(slfl);
  1344.             da[56] = p31.multiply(slfl);
  1345.             da[57] = p41.multiply(slfl);
  1346.             da[58] = p51.multiply(slfl);

  1348. //          fonction g[l] periodique supplementaire
  1349.             fp = fp.
  1350.                  add(da[49].multiply(a[49])).
  1351.                  add(da[50].multiply(a[50])).
  1352.                  add(da[51].multiply(a[51])).
  1353.                  add(da[52].multiply(a[52])).
  1354.                  add(da[53].multiply(a[53])).
  1355.                  add(da[54].multiply(a[54])).
  1356.                  add(da[55].multiply(a[55])).
  1357.                  add(da[56].multiply(a[56])).
  1358.                  add(da[57].multiply(a[57])).
  1359.                  add(da[58].multiply(a[58]));

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

  1364.         }


  1367.         /** Apply a correction coefficient to the given parameter.
  1368.          * @param param the parameter to correct
  1369.          * @return the corrected parameter
  1370.          */
  1371.         private double semestrialCorrection(final double param) {
  1372.             final int debeq_pr = 59;
  1373.             final int debeq_au = 244;
  1374.             final double result;
  1375.             if (day >= 100) {
  1376.                 final double xmult  = (day - debeq_au) / 40.0;
  1377.                 result = param - 2.0 * param * xmult;
  1378.             } else {
  1379.                 final double xmult  = (day - debeq_pr) / 40.0;
  1380.                 result = 2.0 * param * xmult - param;
  1381.             }
  1382.             return result;
  1383.         }

  1385.     }

  1387. }
Created with JaCoCo 0.7.9.201702052155