DTM2000.java

  1. /* Copyright 2002-2016 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;

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

  25. import org.apache.commons.math3.exception.util.DummyLocalizable;
  26. import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
  27. import org.apache.commons.math3.util.FastMath;
  28. import org.orekit.bodies.BodyShape;
  29. import org.orekit.bodies.GeodeticPoint;
  30. import org.orekit.errors.OrekitException;
  31. import org.orekit.errors.OrekitMessages;
  32. import org.orekit.frames.Frame;
  33. import org.orekit.frames.Transform;
  34. import org.orekit.time.AbsoluteDate;
  35. import org.orekit.time.TimeScalesFactory;
  36. import org.orekit.utils.PVCoordinates;
  37. import org.orekit.utils.PVCoordinatesProvider;

  38. /** This atmosphere model is the realization of the DTM-2000 model.
  39.  * <p>
  40.  * It is described in the paper: <br>
  41.  *
  42.  * <b>The DTM-2000 empirical thermosphere model with new data assimilation
  43.  *  and constraints at lower boundary: accuracy and properties</b><br>
  44.  *
  45.  * <i>S. Bruinsma, G. Thuillier and F. Barlier</i> <br>
  46.  *
  47.  * Journal of Atmospheric and Solar-Terrestrial Physics 65 (2003) 1053–1070<br>
  48.  *
  49.  *</p>
  50.  * <p>
  51.  * Two computation methods are proposed to the user:
  52.  * <ul>
  53.  *   <li>one OREKIT independent and compliant with initial FORTRAN
  54.  *       routine entry values:
  55.  *        {@link #getDensity(int, double, double, double, double, double, double, double, double)}.</li>
  56.  *   <li> one compliant with OREKIT Atmosphere interface, necessary to the
  57.  *        {@link org.orekit.forces.drag.DragForce drag force model}
  58.  *        computation.</li>
  59.  * </ul>
  60.  *</p>
  61.  *<p>
  62.  * This model provides dense output for altitudes beyond 120 km. Computed data are:
  63.  * <ul>
  64.  *   <li>Temperature at altitude z (K)</li>
  65.  *   <li>Exospheric temperature above input position (K)</li>
  66.  *   <li>Vertical gradient of T a 120 km</li>
  67.  *   <li>Total density (kg/m³)</li>
  68.  *   <li>Mean atomic mass</li>
  69.  *   <li>Partial densities in (kg/m³) : hydrogen, helium, atomic oxygen,
  70.  *   molecular nitrogen, molecular oxygen, atomic nitrogen</li>
  71.  * </ul>
  72.  * </p>
  73.  * <p>
  74.  * The model needs geographical and time information to compute general values,
  75.  * but also needs space weather data : mean and instantaneous solar flux and
  76.  * geomagnetic indices.
  77.  * </p>
  78.  * <p>
  79.  * Mean solar flux is (for the moment) represented by the F10.7 indices. Instantaneous
  80.  * flux can be set to the mean value if the data is not available. Geomagnetic
  81.  * activity is represented by the Kp indice, which goes from 1 (very low activity) to
  82.  * 9 (high activity).
  83.  * <p>
  84.  * </p>
  85.  * All these data can be found on the <a href="http://sec.noaa.gov/Data/index.html">
  86.  * NOAA (National Oceanic and Atmospheric Administration) website.</a>
  87.  * </p>
  88.  *
  89.  *
  90.  * @author R. Biancale, S. Bruinsma: original fortran routine
  91.  * @author Fabien Maussion (java translation)
  92.  */
  93. public class DTM2000 implements Atmosphere {

  94.     /** Identifier for hydrogen. */
  95.     public static final int HYDROGEN = 1;

  96.     /** Identifier for helium. */
  97.     public static final int HELIUM = 2;

  98.     /** Identifier for atomic oxygen. */
  99.     public static final int ATOMIC_OXYGEN = 3;

  100.     /** Identifier for molecular nitrogen. */
  101.     public static final int MOLECULAR_NITROGEN = 4;

  102.     /** Identifier for molecular oxygen. */
  103.     public static final int MOLECULAR_OXYGEN = 5;

  104.     /** Identifier for atomic nitrogen. */
  105.     public static final int ATOMIC_NITROGEN = 6;

  106.     /** Serializable UID. */
  107.     private static final long serialVersionUID = -8901940398967553588L;

  108.     // Constants :

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

  111.     /** Thermal diffusion coefficient. */
  112.     private static final double[] ALEFA = {
  113.         0, -0.40, -0.38, 0, 0, 0, 0
  114.     };

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

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

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

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

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

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

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

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

  135.     /** Universal gas constant. */
  136.     private static final double RGAS = 831.4;

  137.     /** 2 * π / 365. */
  138.     private static final double ROT = 0.017214206;

  139.     /** 2 * rot. */
  140.     private static final double ROT2 = 0.034428412;

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

  143.     // CHECKSTYLE: stop JavadocVariable check

  144.     /** Elements coefficients. */
  145.     private static double[] tt   = null;
  146.     private static double[] h    = null;
  147.     private static double[] he   = null;
  148.     private static double[] o    = null;
  149.     private static double[] az2  = null;
  150.     private static double[] o2   = null;
  151.     private static double[] az   = null;
  152.     private static double[] t0   = null;
  153.     private static double[] tp   = null;

  154.     /** Partial derivatives. */
  155.     private static double[] dtt  = null;
  156.     private static double[] dh   = null;
  157.     private static double[] dhe  = null;
  158.     private static double[] dox  = null;
  159.     private static double[] daz2 = null;
  160.     private static double[] do2  = null;
  161.     private static double[] daz  = null;
  162.     private static double[] dt0  = null;
  163.     private static double[] dtp  = null;

  164.     // CHECKSTYLE: resume JavadocVariable check

  165.     /** Number of days in current year. */
  166.     private int cachedDay;

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

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

  171.     /** Kp coefficients.
  172.      * <p><ul>
  173.      *   <li>akp[1] = 3-hourly kp</li>
  174.      *   <li>akp[2] = 0 (not used)</li>
  175.      *   <li>akp[3] = mean kp of last 24 hours</li>
  176.      *   <li>akp[4] = 0 (not used)</li>
  177.      * </ul></p>
  178.      */
  179.     private double[] akp = new double[5];

  180.     /** Geodetic altitude in km (minimum altitude: 120 km). */
  181.     private double cachedAlti;

  182.     /** Local solar time (rad). */
  183.     private double cachedHl;

  184.     /** Geodetic Latitude (rad). */
  185.     private double alat;

  186.     /** Geodetic longitude (rad). */
  187.     private double xlon;

  188.     /** Temperature at altitude z (K). */
  189.     private double tz;

  190.     /** Exospheric temperature. */
  191.     private double tinf;

  192.     /** Vertical gradient of T a 120 km. */
  193.     private double tp120;

  194.     /** Total density (g/cm3). */
  195.     private double ro;

  196.     /** Mean atomic mass. */
  197.     private double wmm;

  198.     /** Partial densities in (g/cm3).
  199.      * d(1) = hydrogen
  200.      * d(2) = helium
  201.      * d(3) = atomic oxygen
  202.      * d(4) = molecular nitrogen
  203.      * d(5) = molecular oxygen
  204.      * d(6) = atomic nitrogen
  205.      */
  206.     private double[] d = new double[7];

  207.     // CHECKSTYLE: stop JavadocVariable check

  208.     /** Legendre coefficients. */
  209.     private double p10;
  210.     private double p20;
  211.     private double p30;
  212.     private double p40;
  213.     private double p50;
  214.     private double p60;
  215.     private double p11;
  216.     private double p21;
  217.     private double p31;
  218.     private double p41;
  219.     private double p51;
  220.     private double p22;
  221.     private double p32;
  222.     private double p42;
  223.     private double p52;
  224.     private double p62;
  225.     private double p33;
  226.     private double p10mg;
  227.     private double p20mg;
  228.     private double p40mg;

  229.     /** Local time intermediate values. */
  230.     private double hl0;
  231.     private double ch;
  232.     private double sh;
  233.     private double c2h;
  234.     private double s2h;
  235.     private double c3h;
  236.     private double s3h;

  237.     // CHECKSTYLE: resume JavadocVariable check

  238.     /** Sun position. */
  239.     private PVCoordinatesProvider sun;

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

  242.     /** Earth body shape. */
  243.     private BodyShape earth;

  244.     /** Simple constructor for independent computation.
  245.      * @param parameters the solar and magnetic activity data
  246.      * @param sun the sun position
  247.      * @param earth the earth body shape
  248.      * @exception OrekitException if some resource file reading error occurs
  249.      */
  250.     public DTM2000(final DTM2000InputParameters parameters,
  251.                    final PVCoordinatesProvider sun, final BodyShape earth)
  252.         throws OrekitException {
  253.         this.earth = earth;
  254.         this.sun = sun;
  255.         this.inputParams = parameters;
  256.         if (tt == null) {
  257.             readcoefficients();
  258.         }
  259.     }

  260.     /** {@inheritDoc} */
  261.     public Frame getFrame() {
  262.         return earth.getBodyFrame();
  263.     }

  264.     /** Get the local density with initial entries.
  265.      * @param day day of year
  266.      * @param alti altitude in meters
  267.      * @param lon local longitude (rad)
  268.      * @param lat local latitude (rad)
  269.      * @param hl local solar time in rad (O hr = 0 rad)
  270.      * @param f instantaneous solar flux (F10.7)
  271.      * @param fbar mean solar flux (F10.7)
  272.      * @param akp3 3 hrs geomagnetic activity index (1-9)
  273.      * @param akp24 Mean of last 24 hrs geomagnetic activity index (1-9)
  274.      * @return the local density (kg/m³)
  275.      * @exception OrekitException if altitude is outside of supported range
  276.      */
  277.     public double getDensity(final int day,
  278.                              final double alti, final double lon, final double lat,
  279.                              final double hl, final double f, final double fbar,
  280.                              final double akp3, final double akp24)
  281.         throws OrekitException {
  282.         final double threshold = 120000;
  283.         if (alti < threshold) {
  284.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD,
  285.                                       alti, threshold);
  286.         }
  287.         this.cachedDay  = day;
  288.         this.cachedAlti = alti / 1000;
  289.         xlon = lon;
  290.         alat = lat;
  291.         this.cachedHl = hl;
  292.         this.cachedF[1] = f;
  293.         this.cachedFbar[1] = fbar;
  294.         akp[1] = akp3;
  295.         akp[3] = akp24;
  296.         computation();
  297.         return ro * 1000;
  298.     }


  299.     /** Computes output vales once the inputs are set.
  300.      */
  301.     private void computation() {
  302.         ro = 0.0;

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

  304.         // compute Legendre polynomials wrt geographic pole
  305.         final double c = FastMath.sin(alat);
  306.         final double c2 = c * c;
  307.         final double c4 = c2 * c2;
  308.         final double s = FastMath.cos(alat);
  309.         final double s2 = s * s;
  310.         p10 = c;
  311.         p20 = 1.5 * c2 - 0.5;
  312.         p30 = c * (2.5 * c2 - 1.5);
  313.         p40 = 4.375 * c4 - 3.75 * c2 + 0.375;
  314.         p50 = c * (7.875 * c4 - 8.75 * c2 + 1.875);
  315.         p60 = (5.5 * c * p50 - 2.5 * p40) / 3.0;
  316.         p11 = s;
  317.         p21 = 3.0 * c * s;
  318.         p31 = s * (7.5 * c2 - 1.5);
  319.         p41 = c * s * (17.5 * c2 - 7.5);
  320.         p51 = s * (39.375 * c4 - 26.25 * c2 + 1.875);
  321.         p22 = 3.0 * s2;
  322.         p32 = 15.0 * c * s2;
  323.         p42 = s2 * (52.5 * c2 - 7.5);
  324.         p52 = 3.0 * c * p42 - 2.0 * p32;
  325.         p62 = 2.75 * c * p52 - 1.75 * p42;
  326.         p33 = 15.0 * s * s2;

  327.         // compute Legendre polynomials wrt magnetic pole (79N, 71W)
  328.         final double clmlmg = FastMath.cos(xlon - XLMG);
  329.         final double cmg  = s * CPMG * clmlmg + c * SPMG;
  330.         final double cmg2 = cmg * cmg;
  331.         final double cmg4 = cmg2 * cmg2;
  332.         p10mg = cmg;
  333.         p20mg = 1.5 * cmg2 - 0.5;
  334.         p40mg = 4.375 * cmg4 - 3.75 * cmg2 + 0.375;

  335.         // local time
  336.         hl0 = cachedHl;
  337.         ch  = FastMath.cos(hl0);
  338.         sh  = FastMath.sin(hl0);
  339.         c2h = ch * ch - sh * sh;
  340.         s2h = 2.0 * ch * sh;
  341.         c3h = c2h * ch - s2h * sh;
  342.         s3h = s2h * ch + c2h * sh;

  343.         //  compute function g(l) / tinf, t120, tp120
  344.         int kleq = 1;
  345.         final double gdelt = gFunction(tt, dtt, 1, kleq);
  346.         dtt[1] = 1.0 + gdelt;
  347.         tinf   = tt[1] * dtt[1];

  348.         kleq = 0; // equinox

  349.         if ((cachedDay < 59) || (cachedDay > 284)) {
  350.             kleq = -1; // north winter
  351.         }
  352.         if ((cachedDay > 99) && (cachedDay < 244)) {
  353.             kleq = 1; // north summer
  354.         }

  355.         final double gdelt0 =  gFunction(t0, dt0, 0, kleq);
  356.         dt0[1] = (t0[1] + gdelt0) / t0[1];
  357.         final double t120 = t0[1] + gdelt0;
  358.         final double gdeltp = gFunction(tp, dtp, 0, kleq);
  359.         dtp[1] = (tp[1] + gdeltp) / tp[1];
  360.         tp120 = tp[1] + gdeltp;

  361.         // compute n(z) concentrations: H, He, O, N2, O2, N
  362.         final double sigma   = tp120 / (tinf - t120);
  363.         final double dzeta   = (RE + zlb) / (RE + cachedAlti);
  364.         final double zeta    = (cachedAlti - zlb) * dzeta;
  365.         final double sigzeta = sigma * zeta;
  366.         final double expsz   = FastMath.exp(-sigzeta);
  367.         tz = tinf - (tinf - t120) * expsz;

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

  369.         kleq = 1;

  370.         final double gdelh = gFunction(h, dh, 0, kleq);
  371.         dh[1] = FastMath.exp(gdelh);
  372.         dbase[1] = h[1] * dh[1];

  373.         final double gdelhe = gFunction(he, dhe, 0, kleq);
  374.         dhe[1] = FastMath.exp(gdelhe);
  375.         dbase[2] = he[1] * dhe[1];

  376.         final double gdelo = gFunction(o, dox, 1, kleq);
  377.         dox[1] = FastMath.exp(gdelo);
  378.         dbase[3] = o[1] * dox[1];

  379.         final double gdelaz2 = gFunction(az2, daz2, 1, kleq);
  380.         daz2[1] = FastMath.exp(gdelaz2);
  381.         dbase[4] = az2[1] * daz2[1];

  382.         final double gdelo2 = gFunction(o2, do2, 1, kleq);
  383.         do2[1] = FastMath.exp(gdelo2);
  384.         dbase[5] = o2[1] * do2[1];

  385.         final double gdelaz = gFunction(az, daz, 1, kleq);
  386.         daz[1] = FastMath.exp(gdelaz);
  387.         dbase[6] = az[1] * daz[1];

  388.         final double zlbre  = 1.0 + zlb / RE;
  389.         final double glb    = (GSURF / (zlbre * zlbre)) / (sigma * RGAS * tinf);
  390.         final double t120tz = t120 / tz;

  391.         final double[] cc = new double[7];
  392.         final double[] fz = new double[7];

  393.         for (int i = 1; i <= 6; i++) {
  394.             final double gamma = MA[i] * glb;
  395.             final double upapg = 1.0 + ALEFA[i] + gamma;
  396.             fz[i] = FastMath.pow(t120tz, upapg) * FastMath.exp(-sigzeta * gamma);
  397.             // concentrations of H, He, O, N2, O2, N (particles/cm³)
  398.             cc[i] = dbase[i] * fz[i];
  399.             // densities of H, He, O, N2, O2, N (g/cm³)
  400.             d[i]  = cc[i] * VMA[i];
  401.             // total density
  402.             ro += d[i];
  403.         }

  404.         // mean atomic mass
  405.         wmm = ro / (VMA[1] * (cc[1] + cc[2] + cc[3] + cc[4] + cc[5] + cc[6]));

  406.     }

  407.     /** Computation of function G.
  408.      * @param a vector of coefficients for computation
  409.      * @param da vector of partial derivatives
  410.      * @param ff0 coefficient flag (1 for Ox, Az, He, T°; 0 for H and tp120)
  411.      * @param kle_eq season indicator flag (summer, winter, equinox)
  412.      * @return value of G
  413.      */
  414.     private double gFunction(final double[] a, final double[] da,
  415.                              final int ff0, final int kle_eq) {

  416.         final double[] fmfb   = new double[3];
  417.         final double[] fbm150 = new double[3];

  418.         // latitude terms
  419.         da[2]  = p20;
  420.         da[3]  = p40;
  421.         da[74] = p10;
  422.         double a74 = a[74];
  423.         double a77 = a[77];
  424.         double a78 = a[78];
  425.         if (kle_eq == -1) {
  426.             // winter
  427.             a74 = -a74;
  428.             a77 = -a77;
  429.             a78 = -a78;
  430.         }
  431.         if (kle_eq == 0 ) {
  432.             // equinox
  433.             a74 = semestrialCorrection(a74);
  434.             a77 = semestrialCorrection(a77);
  435.             a78 = semestrialCorrection(a78);
  436.         }
  437.         da[77] = p30;
  438.         da[78] = p50;
  439.         da[79] = p60;

  440.         // flux terms
  441.         fmfb[1]   = cachedF[1] - cachedFbar[1];
  442.         fmfb[2]   = cachedF[2] - cachedFbar[2];
  443.         fbm150[1] = cachedFbar[1] - 150.0;
  444.         fbm150[2] = cachedFbar[2];
  445.         da[4]     = fmfb[1];
  446.         da[6]     = fbm150[1];
  447.         da[4]     = da[4] + a[70] * fmfb[2];
  448.         da[6]     = da[6] + a[71] * fbm150[2];
  449.         da[70]    = fmfb[2] * (a[4] + 2.0 * a[5] * da[4] + a[82] * p10 +
  450.                                a[83] * p20 + a[84] * p30);
  451.         da[71]    = fbm150[2] * (a[6] + 2.0 * a[69] * da[6] + a[85] * p10 +
  452.                                  a[86] * p20 + a[87] * p30);
  453.         da[5]     = da[4] * da[4];
  454.         da[69]    = da[6] * da[6];
  455.         da[82]    = da[4] * p10;
  456.         da[83]    = da[4] * p20;
  457.         da[84]    = da[4] * p30;
  458.         da[85]    = da[6] * p20;
  459.         da[86]    = da[6] * p30;
  460.         da[87]    = da[6] * p40;

  461.         // Kp terms
  462.         final int ikp  = 62;
  463.         final int ikpm = 67;
  464.         final double c2fi = 1.0 - p10mg * p10mg;
  465.         final double dkp  = akp[1] + (a[ikp] + c2fi * a[ikp + 1]) * akp[2];
  466.         double dakp = a[7] + a[8] * p20mg + a[68] * p40mg +
  467.                       2.0 * dkp * (a[60] + a[61] * p20mg +
  468.                                    a[75] * 2.0 * dkp * dkp);
  469.         da[ikp] = dakp * akp[2];
  470.         da[ikp + 1] = da[ikp] * c2fi;
  471.         final double dkpm  = akp[3] + a[ikpm] * akp[4];
  472.         final double dakpm = a[64] + a[65] * p20mg + a[72] * p40mg +
  473.                              2.0 * dkpm * (a[66] + a[73] * p20mg +
  474.                                            a[76] * 2.0 * dkpm * dkpm);
  475.         da[ikpm] = dakpm * akp[4];
  476.         da[7]    = dkp;
  477.         da[8]    = p20mg * dkp;
  478.         da[68]   = p40mg * dkp;
  479.         da[60]   = dkp * dkp;
  480.         da[61]   = p20mg * da[60];
  481.         da[75]   = da[60] * da[60];
  482.         da[64]   = dkpm;
  483.         da[65]   = p20mg * dkpm;
  484.         da[72]   = p40mg * dkpm;
  485.         da[66]   = dkpm * dkpm;
  486.         da[73]   = p20mg * da[66];
  487.         da[76]   = da[66] * da[66];

  488.         // non-periodic g(l) function
  489.         double f0 = a[4]  * da[4]  + a[5]  * da[5]  + a[6]  * da[6]  +
  490.                     a[69] * da[69] + a[82] * da[82] + a[83] * da[83] +
  491.                     a[84] * da[84] + a[85] * da[85] + a[86] * da[86] +
  492.                     a[87] * da[87];
  493.         final double f1f = 1.0 + f0 * ff0;

  494.         f0 = f0 + a[2] * da[2] + a[3] * da[3] + a74 * da[74] +
  495.              a77 * da[77] + a[7] * da[7] + a[8] * da[8] +
  496.              a[60] * da[60] + a[61] * da[61] + a[68] * da[68] +
  497.              a[64] * da[64] + a[65] * da[65] + a[66] * da[66] +
  498.              a[72] * da[72] + a[73] * da[73] + a[75] * da[75] +
  499.              a[76] * da[76] + a78   * da[78] + a[79] * da[79];
  500. //      termes annuels symetriques en latitude
  501.         da[9]  = FastMath.cos(ROT * (cachedDay - a[11]));
  502.         da[10] = p20 * da[9];
  503. //      termes semi-annuels symetriques en latitude
  504.         da[12] = FastMath.cos(ROT2 * (cachedDay - a[14]));
  505.         da[13] = p20 * da[12];
  506. //      termes annuels non symetriques en latitude
  507.         final double coste = FastMath.cos(ROT * (cachedDay - a[18]));
  508.         da[15] = p10 * coste;
  509.         da[16] = p30 * coste;
  510.         da[17] = p50 * coste;
  511. //      terme  semi-annuel  non symetrique  en latitude
  512.         final double cos2te = FastMath.cos(ROT2 * (cachedDay - a[20]));
  513.         da[19] = p10 * cos2te;
  514.         da[39] = p30 * cos2te;
  515.         da[59] = p50 * cos2te;
  516. //      termes diurnes [et couples annuel]
  517.         da[21] = p11 * ch;
  518.         da[22] = p31 * ch;
  519.         da[23] = p51 * ch;
  520.         da[24] = da[21] * coste;
  521.         da[25] = p21 * ch * coste;
  522.         da[26] = p11 * sh;
  523.         da[27] = p31 * sh;
  524.         da[28] = p51 * sh;
  525.         da[29] = da[26] * coste;
  526.         da[30] = p21 * sh * coste;
  527. //      termes semi-diurnes [et couples annuel]
  528.         da[31] = p22 * c2h;
  529.         da[37] = p42 * c2h;
  530.         da[32] = p32 * c2h * coste;
  531.         da[33] = p22 * s2h;
  532.         da[38] = p42 * s2h;
  533.         da[34] = p32 * s2h * coste;
  534.         da[88] = p32 * c2h;
  535.         da[89] = p32 * s2h;
  536.         da[90] = p52 * c2h;
  537.         da[91] = p52 * s2h;
  538.         double a88 = a[88];
  539.         double a89 = a[89];
  540.         double a90 = a[90];
  541.         double a91 = a[91];
  542.         if (kle_eq == -1) {            //hiver
  543.             a88 = -a88;
  544.             a89 = -a89;
  545.             a90 = -a90;
  546.             a91 = -a91;
  547.         }
  548.         if (kle_eq == 0) {             //equinox
  549.             a88 = semestrialCorrection(a88);
  550.             a89 = semestrialCorrection(a89);
  551.             a90 = semestrialCorrection(a90);
  552.             a91 = semestrialCorrection(a91);
  553.         }
  554.         da[92] = p62 * c2h;
  555.         da[93] = p62 * s2h;
  556. //      termes ter-diurnes
  557.         da[35] = p33 * c3h;
  558.         da[36] = p33 * s3h;
  559. //      fonction g[l] periodique
  560.         double fp = a[9]  * da[9]  + a[10] * da[10] + a[12] * da[12] + a[13] * da[13] +
  561.                     a[15] * da[15] + a[16] * da[16] + a[17] * da[17] + a[19] * da[19] +
  562.                     a[21] * da[21] + a[22] * da[22] + a[23] * da[23] + a[24] * da[24] +
  563.                     a[25] * da[25] + a[26] * da[26] + a[27] * da[27] + a[28] * da[28] +
  564.                     a[29] * da[29] + a[30] * da[30] + a[31] * da[31] + a[32] * da[32] +
  565.                     a[33] * da[33] + a[34] * da[34] + a[35] * da[35] + a[36] * da[36] +
  566.                     a[37] * da[37] + a[38] * da[38] + a[39] * da[39] + a[59] * da[59] +
  567.                     a88   * da[88] + a89   * da[89] + a90   * da[90] + a91   * da[91] +
  568.                     a[92] * da[92] + a[93] * da[93];
  569. //      termes d'activite magnetique
  570.         da[40] = p10 * coste * dkp;
  571.         da[41] = p30 * coste * dkp;
  572.         da[42] = p50 * coste * dkp;
  573.         da[43] = p11 * ch * dkp;
  574.         da[44] = p31 * ch * dkp;
  575.         da[45] = p51 * ch * dkp;
  576.         da[46] = p11 * sh * dkp;
  577.         da[47] = p31 * sh * dkp;
  578.         da[48] = p51 * sh * dkp;

  579. //      fonction g[l] periodique supplementaire
  580.         fp += a[40] * da[40] + a[41] * da[41] + a[42] * da[42] + a[43] * da[43] +
  581.               a[44] * da[44] + a[45] * da[45] + a[46] * da[46] + a[47] * da[47] +
  582.               a[48] * da[48];

  583.         dakp = (a[40] * p10 + a[41] * p30 + a[42] * p50) * coste +
  584.                (a[43] * p11 + a[44] * p31 + a[45] * p51) * ch +
  585.                (a[46] * p11 + a[47] * p31 + a[48] * p51) * sh;
  586.         da[ikp] += dakp * akp[2];
  587.         da[ikp + 1] = da[ikp] + dakp * c2fi * akp[2];
  588. //      termes de longitude
  589.         final double clfl = FastMath.cos(xlon);
  590.         da[49] = p11 * clfl;
  591.         da[50] = p21 * clfl;
  592.         da[51] = p31 * clfl;
  593.         da[52] = p41 * clfl;
  594.         da[53] = p51 * clfl;
  595.         final double slfl = FastMath.sin(xlon);
  596.         da[54] = p11 * slfl;
  597.         da[55] = p21 * slfl;
  598.         da[56] = p31 * slfl;
  599.         da[57] = p41 * slfl;
  600.         da[58] = p51 * slfl;

  601. //      fonction g[l] periodique supplementaire
  602.         fp += a[49] * da[49] + a[50] * da[50] + a[51] * da[51] + a[52] * da[52] +
  603.               a[53] * da[53] + a[54] * da[54] + a[55] * da[55] + a[56] * da[56] +
  604.               a[57] * da[57] + a[58] * da[58];

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

  607.     }


  608.     /** Apply a correction coefficient to the given parameter.
  609.      * @param param the parameter to correct
  610.      * @return the corrected parameter
  611.      */
  612.     private double semestrialCorrection(final double param) {
  613.         final int debeq_pr = 59;
  614.         final int debeq_au = 244;
  615.         double xmult;
  616.         double result;
  617.         if (cachedDay >= 100) {
  618.             xmult  = (cachedDay - debeq_au) / 40.0;
  619.             result = param - 2.0 * param * xmult;
  620.         } else {
  621.             xmult  = (cachedDay - debeq_pr) / 40.0;
  622.             result = 2.0 * param * xmult - param;
  623.         }
  624.         return result;
  625.     }


  626.     /** Store the DTM model elements coefficients in internal arrays.
  627.      * @exception OrekitException if some resource file reading error occurs
  628.      */
  629.     private static void readcoefficients() throws OrekitException {

  630.         final int size = NLATM + 1;
  631.         tt   = new double[size];
  632.         h    = new double[size];
  633.         he   = new double[size];
  634.         o    = new double[size];
  635.         az2  = new double[size];
  636.         o2   = new double[size];
  637.         az   = new double[size];
  638.         t0   = new double[size];
  639.         tp   = new double[size];
  640.         dtt  = new double[size];
  641.         dh   = new double[size];
  642.         dhe  = new double[size];
  643.         dox  = new double[size];
  644.         daz2 = new double[size];
  645.         do2  = new double[size];
  646.         daz  = new double[size];
  647.         dt0  = new double[size];
  648.         dtp  = new double[size];

  649.         Arrays.fill(dtt,  Double.NaN);
  650.         Arrays.fill(dh,   Double.NaN);
  651.         Arrays.fill(dhe,  Double.NaN);
  652.         Arrays.fill(dox,  Double.NaN);
  653.         Arrays.fill(daz2, Double.NaN);
  654.         Arrays.fill(do2,  Double.NaN);
  655.         Arrays.fill(daz,  Double.NaN);
  656.         Arrays.fill(dt0,  Double.NaN);
  657.         Arrays.fill(dtp,  Double.NaN);

  658.         final InputStream in = DTM2000.class.getResourceAsStream(DTM2000);
  659.         if (in == null) {
  660.             throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_RESOURCE, DTM2000);
  661.         }

  662.         BufferedReader r = null;
  663.         try {

  664.             r = new BufferedReader(new InputStreamReader(in, "UTF-8"));
  665.             r.readLine();
  666.             r.readLine();
  667.             for (String line = r.readLine(); line != null; line = r.readLine()) {
  668.                 final int num = Integer.parseInt(line.substring(0, 4).replace(' ', '0'));
  669.                 line = line.substring(4);
  670.                 tt[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  671.                 line = line.substring(13 + 9);
  672.                 h[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  673.                 line = line.substring(13 + 9);
  674.                 he[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  675.                 line = line.substring(13 + 9);
  676.                 o[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  677.                 line = line.substring(13 + 9);
  678.                 az2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  679.                 line = line.substring(13 + 9);
  680.                 o2[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  681.                 line = line.substring(13 + 9);
  682.                 az[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  683.                 line = line.substring(13 + 9);
  684.                 t0[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  685.                 line = line.substring(13 + 9);
  686.                 tp[num] = Double.parseDouble(line.substring(0, 13).replace(' ', '0'));
  687.             }
  688.         } catch (IOException ioe) {
  689.             throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
  690.         } finally {
  691.             if (r != null) {
  692.                 try {
  693.                     r.close();
  694.                 } catch (IOException ioe) {
  695.                     throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
  696.                 }
  697.             }
  698.         }
  699.     }



  700.     /** Get the current exospheric temperature above input position.
  701.      * {@link #getDensity(int, double, double, double, double, double, double, double, double) getDensity}
  702.      * method <b>must</b> be called before calling this function.
  703.      * @return the exospheric temperature (K)
  704.      */
  705.     public double getTinf() {
  706.         return tinf;
  707.     }

  708.     /** Get the local temperature.
  709.      * {@link #getDensity(int, double, double, double, double, double, double, double, double) getDensity}
  710.      * method <b>must</b> be called before calling this function.
  711.      * @return the temperature at altitude z (K)
  712.      */
  713.     public double getT() {
  714.         return tz;
  715.     }

  716.     /** Get the local mean atomic mass.
  717.      * {@link #getDensity(int, double, double, double, double, double, double, double, double) getDensity}
  718.      * method <b>must</b> be called before calling this function.
  719.      * @return the local mean atomic mass
  720.      */
  721.     public double getMam() {
  722.         return wmm;
  723.     }

  724.     /** Get the local partial density of the selected element.
  725.      * {@link #getDensity(int, double, double, double, double, double, double, double, double) getDensity}
  726.      * method <b>must</b> be called before calling this function.
  727.      * @param identifier one of the six elements : {@link #HYDROGEN}, {@link #HELIUM},
  728.      * {@link #ATOMIC_OXYGEN}, {@link #MOLECULAR_NITROGEN},  {@link #MOLECULAR_OXYGEN}, {@link #ATOMIC_NITROGEN}
  729.      * @return the local partial density (kg/m³)
  730.      */
  731.     public double getPartialDensities(final int identifier) {
  732.         if (identifier < 1 || identifier > 6) {
  733.             throw new IllegalArgumentException("element identifier is not correct");
  734.         }
  735.         return d[identifier] * 1000;
  736.     }

  737.     /** Get the local density.
  738.      * @param date current date
  739.      * @param position current position in frame
  740.      * @param frame the frame in which is defined the position
  741.      * @return local density (kg/m³)
  742.      * @exception OrekitException if date is out of range of solar activity model
  743.      * or if some frame conversion cannot be performed
  744.      */
  745.     public double getDensity(final AbsoluteDate date, final Vector3D position,
  746.                              final Frame frame)
  747.         throws OrekitException {

  748.         // check if data are available :
  749.         if ((date.compareTo(inputParams.getMaxDate()) > 0) ||
  750.             (date.compareTo(inputParams.getMinDate()) < 0)) {
  751.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  752.                                       date, inputParams.getMinDate(), inputParams.getMaxDate());
  753.         }

  754.         // compute day number in current year
  755.         final Calendar cal = new GregorianCalendar();
  756.         cal.setTime(date.toDate(TimeScalesFactory.getUTC()));
  757.         final int day = cal.get(Calendar.DAY_OF_YEAR);
  758.         //position in ECEF so we only have to do the transform once
  759.         final Frame ecef = earth.getBodyFrame();
  760.         final Vector3D pEcef = frame.getTransformTo(ecef, date)
  761.                 .transformPosition(position);
  762.         // compute geodetic position
  763.         final GeodeticPoint inBody = earth.transform(pEcef, ecef, date);
  764.         final double alti = inBody.getAltitude();
  765.         final double lon = inBody.getLongitude();
  766.         final double lat = inBody.getLatitude();

  767.         // compute local solar time
  768.         final Vector3D sunPos = sun.getPVCoordinates(date, ecef).getPosition();
  769.         final double hl = FastMath.PI + FastMath.atan2(
  770.                 sunPos.getX() * pEcef.getY() - sunPos.getY() * pEcef.getX(),
  771.                 sunPos.getX() * pEcef.getX() + sunPos.getY() * pEcef.getY());

  772.         // get current solar activity data and compute
  773.         return getDensity(day, alti, lon, lat, hl, inputParams.getInstantFlux(date),
  774.                           inputParams.getMeanFlux(date), inputParams.getThreeHourlyKP(date),
  775.                           inputParams.get24HoursKp(date));

  776.     }

  777.     /** Get the inertial velocity of atmosphere molecules.
  778.      * Here the case is simplified : atmosphere is supposed to have a null velocity
  779.      * in earth frame.
  780.      * @param date current date
  781.      * @param position current position in frame
  782.      * @param frame the frame in which is defined the position
  783.      * @return velocity (m/s) (defined in the same frame as the position)
  784.      * @exception OrekitException if some frame conversion cannot be performed
  785.      */
  786.     public Vector3D getVelocity(final AbsoluteDate date, final Vector3D position,
  787.                                 final Frame frame)
  788.         throws OrekitException {
  789.         final Transform bodyToFrame = earth.getBodyFrame().getTransformTo(frame, date);
  790.         final Vector3D posInBody = bodyToFrame.getInverse().transformPosition(position);
  791.         final PVCoordinates pvBody = new PVCoordinates(posInBody, new Vector3D(0, 0, 0));
  792.         final PVCoordinates pvFrame = bodyToFrame.transformPVCoordinates(pvBody);
  793.         return pvFrame.getVelocity();
  794.     }

  795. }