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