JB2006.java

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

  19. import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
  20. import org.apache.commons.math3.util.FastMath;
  21. import org.orekit.bodies.BodyShape;
  22. import org.orekit.bodies.GeodeticPoint;
  23. import org.orekit.errors.OrekitException;
  24. import org.orekit.errors.OrekitMessages;
  25. import org.orekit.frames.Frame;
  26. import org.orekit.frames.Transform;
  27. import org.orekit.time.AbsoluteDate;
  28. import org.orekit.utils.Constants;
  29. import org.orekit.utils.PVCoordinates;
  30. import org.orekit.utils.PVCoordinatesProvider;

  32. /** This is the realization of the Jacchia-Bowman 2006 atmospheric model.
  33.  * <p>
  34.  * It is described in the paper: <br>
  35.  *
  36.  * <a href="http://sol.spacenvironment.net/~JB2006/pubs/JB2006_AIAA-6166_model.pdf">A
  37.  * New Empirical Thermospheric Density Model JB2006 Using New Solar Indices</a><br>
  38.  *
  39.  * <i>Bruce R. Bowman, W. Kent Tobiska and Frank A. Marcos</i> <br>
  40.  *
  41.  * AIAA 2006-6166<br>
  42.  *</p>
  43.  * <p>
  44.  * Two computation methods are proposed to the user:
  45.  * <ul>
  46.  * <li> one OREKIT independent and compliant with initial FORTRAN routine entry values:
  47.  *        {@link #getDensity(double, double, double, double, double, double, double, double, double, double, double, double, double)}. </li>
  48.  * <li> one compliant with OREKIT Atmosphere interface, necessary to the
  49.  *        {@link org.orekit.forces.drag.DragForce
  50.  *        drag force model} computation.</li>
  51.  * </ul>
  52.  * </p>
  53.  * <p>
  54.  * This model provides dense output for all altitudes and positions. Output data are :
  55.  * <ul>
  56.  * <li>Exospheric Temperature above Input Position (deg K)</li>
  57.  * <li>Temperature at Input Position (deg K)</li>
  58.  * <li>Total Mass-Density at Input Position (kg/m<sup>3</sup>)</li>
  59.  * </ul>
  60.  * </p>
  61.  * <p>
  62.  * The model needs geographical and time information to compute general values,
  63.  * but also needs space weather data : mean and daily solar flux, retrieved threw
  64.  * different indices, and planetary geomagnetic indices. <br>
  65.  * More information on these indices can be found on the  <a
  66.  * href="http://sol.spacenvironment.net/~JB2006/JB2006_index.html">
  67.  * official JB2006 website.</a>
  68.  *</p>
  69.  *
  70.  * @author Bruce R Bowman (HQ AFSPC, Space Analysis Division), Feb 2006: FORTRAN routine
  71.  * @author Fabien Maussion (java translation)
  72.  */
  73. public class JB2006 implements Atmosphere {

  75.     /** Serializable UID. */
  76.     private static final long serialVersionUID = -4201270765122160831L;

  78.     /** The alpha are the thermal diffusion coefficients in equation (6). */
  79.     private static final double[] ALPHA = {
  80.         0, 0, 0, 0, 0, -0.38
  81.     };

  83.     /** Natural logarithm of 10.0. */
  84.     private static final double AL10  = 2.3025851;

  86.     /** Molecular weights in order: N2, O2, O, Ar, He and H. */
  87.     private static final double[] AMW = {
  88.         0,
  89.         28.0134, 31.9988, 15.9994, 39.9480, 4.0026, 1.00797
  90.     };

  92.     /** Avogadro's number in mks units (molecules/kmol). */
  93.     private static final double AVOGAD = 6.02257e26;

  95.     /** Approximate value for 2 &pi;. */
  96.     private static final double TWOPI   = 6.2831853;

  98.     /** Approximate value for &pi;. */
  99.     private static final double PI      = 3.1415927;

  101.     /** Approximate value for &pi; / 2. */
  102.     private static final double PIOV2   = 1.5707963;

  104.     /** The FRAC are the assumed sea-level volume fractions in order: N2, O2, Ar, and He. */
  105.     private static final double[] FRAC = {
  106.         0,
  107.         0.78110, 0.20955, 9.3400e-3, 1.2890e-5
  108.     };

  110.     /** Universal gas-constant in mks units (joules/K/kmol). */
  111.     private static final double RSTAR = 8314.32;

  113.     /** Value used to establish height step sizes in the regime 90km to 105km. */
  114.     private static final double R1 = 0.010;

  116.     /** Value used to establish height step sizes in the regime 105km to 500km. */
  117.     private static final double R2 = 0.025;

  119.     /** Value used to establish height step sizes in the regime above 500km. */
  120.     private static final double R3 = 0.075;

  122.     /** Weights for the Newton-Cotes five-points quadrature formula. */
  123.     private static final double[] WT = {
  124.         0,
  125.         0.311111111111111, 1.422222222222222,
  126.         0.533333333333333, 1.422222222222222,
  127.         0.311111111111111
  128.     };

  130.     /** Coefficients for high altitude density correction. */
  131.     private static final double[] CHT = {
  132.         0, 0.22, -0.20e-02, 0.115e-02, -0.211e-05
  133.     };

  135.     /** FZ global model values (1978-2004 fit).  */
  136.     private static final double[] FZM = {
  137.         0,
  138.         0.111613e+00, -0.159000e-02, 0.126190e-01,
  139.         -0.100064e-01, -0.237509e-04, 0.260759e-04
  140.     };

  142.     /** GT global model values (1978-2004 fit). */
  143.     private static final double[] GTM = {
  144.         0,
  145.         -0.833646e+00, -0.265450e+00, 0.467603e+00, -0.299906e+00,
  146.         -0.105451e+00, -0.165537e-01, -0.380037e-01, -0.150991e-01,
  147.         -0.541280e-01, 0.119554e-01, 0.437544e-02, -0.369016e-02,
  148.         0.206763e-02, -0.142888e-02, -0.867124e-05, 0.189032e-04,
  149.         0.156988e-03,  0.491286e-03, -0.391484e-04, -0.126854e-04,
  150.         0.134078e-04, -0.614176e-05, 0.343423e-05
  151.     };

  153.     /** XAMBAR relative data. */
  154.     private static final double[] CXAMB = {
  155.         0,
  156.         28.15204, -8.5586e-2, +1.2840e-4, -1.0056e-5,
  157.         -1.0210e-5, +1.5044e-6, +9.9826e-8
  158.     };

  160.     /** DTSUB relative data. */
  161.     private static final double[] BDT_SUB = {
  162.         0,
  163.         -0.457512297e+01, -0.512114909e+01, -0.693003609e+02,
  164.         0.203716701e+03,  0.703316291e+03, -0.194349234e+04,
  165.         0.110651308e+04, -0.174378996e+03,  0.188594601e+04,
  166.         -0.709371517e+04,  0.922454523e+04, -0.384508073e+04,
  167.         -0.645841789e+01,  0.409703319e+02, -0.482006560e+03,
  168.         0.181870931e+04, -0.237389204e+04,  0.996703815e+03,
  169.         0.361416936e+02
  170.     };

  172.     /** DTSUB relative data. */
  173.     private static final double[] CDT_SUB = {
  174.         0,
  175.         -0.155986211e+02, -0.512114909e+01, -0.693003609e+02,
  176.         0.203716701e+03,  0.703316291e+03, -0.194349234e+04,
  177.         0.110651308e+04, -0.220835117e+03,  0.143256989e+04,
  178.         -0.318481844e+04,  0.328981513e+04, -0.135332119e+04,
  179.         0.199956489e+02, -0.127093998e+02,  0.212825156e+02,
  180.         -0.275555432e+01,  0.110234982e+02,  0.148881951e+03,
  181.         -0.751640284e+03,  0.637876542e+03,  0.127093998e+02,
  182.         -0.212825156e+02,  0.275555432e+01
  183.     };

  185.     /** Temperatures.
  186.      *  <p><ul>
  187.      *  <li>TEMP(1): Exospheric Temperature above Input Position (deg K)</li>
  188.      *  <li>TEMP(2): Temperature at Input Position (deg K)</li>
  189.      *  </ul></p>
  190.      */
  191.     private double[] temp = new double[3];

  193.     /** Total Mass-Density at Input Position (kg/m<sup>3</sup>). */
  194.     private double rho;

  196.     /** Sun position. */
  197.     private PVCoordinatesProvider sun;

  199.     /** External data container. */
  200.     private JB2006InputParameters inputParams;

  202.     /** Earth body shape. */
  203.     private BodyShape earth;

  205.     /** Constructor with space environment information for internal computation.
  206.      * @param parameters the solar and magnetic activity data
  207.      * @param sun the sun position
  208.      * @param earth the earth body shape
  209.      */
  210.     public JB2006(final JB2006InputParameters parameters,
  211.                   final PVCoordinatesProvider sun, final BodyShape earth) {
  212.         this.earth = earth;
  213.         this.sun = sun;
  214.         this.inputParams = parameters;
  215.     }

  217.     /** {@inheritDoc} */
  218.     public Frame getFrame() {
  219.         return earth.getBodyFrame();
  220.     }

  222.     /** Get the local density with initial entries.
  223.      * @param dateMJD date and time, in modified julian days and fraction
  224.      * @param sunRA Right Ascension of Sun (radians)
  225.      * @param sunDecli Declination of Sun (radians)
  226.      * @param satLon Right Ascension of position (radians)
  227.      * @param satLat Geocentric latitude of position (radians)
  228.      * @param satAlt Height of position (m)
  229.      * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m<sup>2</sup>*Hertz)).
  230.      *            Tabular time 1.0 day earlier
  231.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time
  232.      * @param ap Geomagnetic planetary 3-hour index A<sub>p</sub>
  233.      *            for a tabular time 6.7 hours earlier
  234.      * @param s10 EUV index (26-34 nm) scaled to F10. Tabular time 1 day earlier.
  235.      * @param s10B UV 81-day averaged centered index
  236.      * @param xm10 MG2 index scaled to F10
  237.      * @param xm10B MG2 81-day ave. centered index. Tabular time 5.0 days earlier.
  238.      * @return total mass-Density at input position (kg/m<sup>3</sup>)
  239.      */
  240.     public double getDensity(final double dateMJD, final double sunRA, final double sunDecli,
  241.                              final double satLon, final double satLat, final double satAlt,
  242.                              final double f10, final double f10B, final double ap,
  243.                              final double s10, final double s10B, final double xm10, final double xm10B) {

  245.         final double scaledSatAlt = satAlt / 1000.0;

  247.         // Equation (14)
  248.         final double tc = 379 + 3.353 * f10B + 0.358 * (f10 - f10B) +
  249.                           2.094 * (s10 - s10B) + 0.343 * (xm10 - xm10B);

  251.         // Equation (15)
  252.         final double eta =   0.5 * FastMath.abs(satLat - sunDecli);
  253.         final double theta = 0.5 * FastMath.abs(satLat + sunDecli);

  255.         // Equation (16)
  256.         final double h     = satLon - sunRA;
  257.         final double tau   = h - 0.64577182 + 0.10471976 * FastMath.sin(h + 0.75049158);
  258.         double solTimeHour = FastMath.toDegrees(h + PI) / 15.0;
  259.         if (solTimeHour >= 24) {
  260.             solTimeHour = solTimeHour - 24.;
  261.         }
  262.         if (solTimeHour < 0) {
  263.             solTimeHour = solTimeHour + 24.;
  264.         }

  266.         // Equation (17)
  267.         final double C = FastMath.pow(FastMath.cos(eta), 2.5);
  268.         final double S = FastMath.pow(FastMath.sin(theta), 2.5);
  269.         final double tmp = FastMath.abs(FastMath.cos(0.5 * tau));
  270.         final double DF = S + (C - S) * tmp * tmp * tmp;
  271.         final double TSUBL = tc * (1. + 0.31 * DF);

  273.         // Equation (18)
  274.         final double EXPAP = FastMath.exp(-0.08 * ap);
  275.         final double DTG = ap + 100. * (1. - EXPAP);

  277.         // Compute correction to dTc for local solar time and lat correction
  278.         final double DTCLST = dTc(f10, solTimeHour, satLat, scaledSatAlt);

  280.         // Compute the local exospheric temperature.
  281.         final double TINF = TSUBL + DTG + DTCLST;
  282.         temp[1] = TINF;

  284.         // Equation (9)
  285.         final double TSUBX = 444.3807 + 0.02385 * TINF - 392.8292 * FastMath.exp(-0.0021357 * TINF);

  287.         // Equation (11)
  288.         final double GSUBX = 0.054285714 * (TSUBX - 183.);

  290.         // The TC array will be an argument in the call to
  291.         // XLOCAL, which evaluates Equation (10) or Equation (13)
  292.         final double[] TC = new double[4];
  293.         TC[0] = TSUBX;
  294.         TC[1] = GSUBX;

  296.         //   A AND GSUBX/A OF Equation (13)
  297.         TC[2] = (TINF - TSUBX) / PIOV2;
  298.         TC[3] = GSUBX / TC[2];

  300.         // Equation (5)
  301.         final double Z1 = 90.;
  302.         final double Z2 = FastMath.min(scaledSatAlt, 105.0);
  303.         double AL = FastMath.log(Z2 / Z1);
  304.         int N = (int) FastMath.floor(AL / R1) + 1;
  305.         double ZR = FastMath.exp(AL / N);
  306.         final double AMBAR1 = xAmbar(Z1);
  307.         final double TLOC1 = xLocal(Z1, TC);
  308.         double ZEND   = Z1;
  309.         double SUM2   = 0.;
  310.         double AIN    = AMBAR1 * xGrav(Z1) / TLOC1;
  311.         double AMBAR2 = 0;
  312.         double TLOC2  = 0;
  313.         double Z      = 0;
  314.         double GRAVL  = 0;

  316.         for (int i = 1; i <= N; ++i) {
  317.             Z = ZEND;
  318.             ZEND = ZR * Z;
  319.             final double DZ = 0.25 * (ZEND - Z);
  320.             double SUM1 = WT[1] * AIN;
  321.             for (int j = 2; j <= 5; ++j) {
  322.                 Z += DZ;
  323.                 AMBAR2 = xAmbar(Z);
  324.                 TLOC2  = xLocal(Z, TC);
  325.                 GRAVL  = xGrav(Z);
  326.                 AIN    = AMBAR2 * GRAVL / TLOC2;
  327.                 SUM1  += WT[j] * AIN;
  328.             }
  329.             SUM2 = SUM2 + DZ * SUM1;
  330.         }
  331.         final double FACT1 = 1000.0 / RSTAR;
  332.         rho = 3.46e-6 * AMBAR2 * TLOC1 * FastMath.exp(-FACT1 * SUM2) / (AMBAR1 * TLOC2);

  334.         // Equation (2)
  335.         final double ANM = AVOGAD * rho;
  336.         double AN  = ANM / AMBAR2;

  338.         // Equation (3)
  339.         double FACT2  = ANM / 28.960;
  340.         final double[] ALN = new double[7];
  341.         ALN[1] = FastMath.log(FRAC[1] * FACT2);
  342.         ALN[4] = FastMath.log(FRAC[3] * FACT2);
  343.         ALN[5] = FastMath.log(FRAC[4] * FACT2);

  345.         // Equation (4)
  346.         ALN[2] = FastMath.log(FACT2 * (1. + FRAC[2]) - AN);
  347.         ALN[3] = FastMath.log(2. * (AN - FACT2));

  349.         if (scaledSatAlt <= 105.0) {
  350.             temp[2] = TLOC2;
  351.             // Put in negligible hydrogen for use in DO-LOOP 13
  352.             ALN[6] = ALN[5] - 25.0;
  353.         } else {
  354.             // Equation (6)
  355.             final double Z3 = FastMath.min(scaledSatAlt, 500.0);
  356.             AL   = FastMath.log(Z3 / Z);
  357.             N    = (int) FastMath.floor(AL / R2) + 1;
  358.             ZR   = FastMath.exp(AL / N);
  359.             SUM2 = 0.;
  360.             AIN  = GRAVL / TLOC2;

  362.             double TLOC3 = 0;
  363.             for (int I = 1; I <= N; ++I) {
  364.                 Z = ZEND;
  365.                 ZEND = ZR * Z;
  366.                 final double DZ = 0.25 * (ZEND - Z);
  367.                 double SUM1 = WT[1] * AIN;
  368.                 for (int J = 2; J <= 5; ++J) {
  369.                     Z    += DZ;
  370.                     TLOC3 = xLocal(Z, TC);
  371.                     GRAVL = xGrav(Z);
  372.                     AIN   = GRAVL / TLOC3;
  373.                     SUM1  = SUM1 + WT[J] * AIN;
  374.                 }
  375.                 SUM2 = SUM2 + DZ * SUM1;
  376.             }

  378.             final double Z4 = FastMath.max(scaledSatAlt, 500.0);
  379.             AL = FastMath.log(Z4 / Z);
  380.             double R = R2;
  381.             if (scaledSatAlt > 500.0) {
  382.                 R = R3;
  383.             }
  384.             N = (int) FastMath.floor(AL / R) + 1;
  385.             ZR = FastMath.exp(AL / N);
  386.             double SUM3 = 0.;
  387.             double TLOC4 = 0;
  388.             for (int I = 1; I <= N; ++I) {
  389.                 Z = ZEND;
  390.                 ZEND = ZR * Z;
  391.                 final double DZ = 0.25 * (ZEND - Z);
  392.                 double SUM1 = WT[1] * AIN;
  393.                 for (int J = 2; J <= 5; ++J) {
  394.                     Z    += DZ;
  395.                     TLOC4 = xLocal(Z, TC);
  396.                     GRAVL = xGrav(Z);
  397.                     AIN   = GRAVL / TLOC4;
  398.                     SUM1  = SUM1 + WT[J] * AIN;
  399.                 }
  400.                 SUM3 = SUM3 + DZ * SUM1;
  401.             }
  402.             double ALTR;
  403.             double HSIGN;
  404.             if (scaledSatAlt <= 500.) {
  405.                 temp[2] = TLOC3;
  406.                 ALTR = FastMath.log(TLOC3 / TLOC2);
  407.                 FACT2 = FACT1 * SUM2;
  408.                 HSIGN = 1.0;

  410.             } else {
  411.                 temp[2] = TLOC4;
  412.                 ALTR = FastMath.log(TLOC4 / TLOC2);
  413.                 FACT2 = FACT1 * (SUM2 + SUM3);
  414.                 HSIGN = -1.0;
  415.             }
  416.             for (int I = 1; I <= 5; ++I) {
  417.                 ALN[I] = ALN[I] - (1.0 + ALPHA[I]) * ALTR - FACT2 * AMW[I];
  418.             }

  420.             // Equation (7) - Note that in CIRA72, AL10T5 = DLOG10(T500)
  421.             final double AL10T5 = FastMath.log(TINF) / FastMath.log(10);
  422.             final double ALNH5 = (5.5 * AL10T5 - 39.40) * AL10T5 + 73.13;
  423.             ALN[6] = AL10 * (ALNH5 + 6.) + HSIGN * (FastMath.log(TLOC4 / TLOC3) + FACT1 * SUM3 * AMW[6]);

  425.         }

  427.         // Equation (24)  - J70 Seasonal-Latitudinal Variation
  428.         final double CAPPHI = ((dateMJD - 36204.0) / 365.2422) % 1;
  429.         final int signum = (satLat >= 0) ? 1 : -1;
  430.         final double sinLat = FastMath.sin(satLat);
  431.         final double DLRSL = 0.02 * (scaledSatAlt - 90.) * FastMath.exp(-0.045 * (scaledSatAlt - 90.)) *
  432.                              signum * FastMath.sin(TWOPI * CAPPHI + 1.72) * sinLat * sinLat;

  434.         // Equation (23) - Computes the semiannual variation
  435.         double DLRSA = 0;
  436.         if (Z < 2000.0) {
  437.             final double D1950 = dateMJD - 33281.0;
  438.             // Use new semiannual model DELTA LOG RHO
  439.             DLRSA = semian(dayOfYear(D1950), scaledSatAlt, f10B);
  440.         }

  442.         // Sum the delta-log-rhos and apply to the number densities.
  443.         // In CIRA72 the following equation contains an actual sum,
  444.         // namely DLR = AL10 * (DLRGM + DLRSA + DLRSL)
  445.         // However, for Jacchia 70, there is no DLRGM or DLRSA.
  446.         final double DLR = AL10 * (DLRSL + DLRSA);
  447.         for (int i = 1; i <= 6; ++i) {
  448.             ALN[i] += DLR;
  449.         }

  451.         // Compute mass-density and mean-molecular-weight and
  452.         // convert number density logs from natural to common.

  454.         double SUMNM = 0.0;

  456.         for (int I = 1; I <= 6; ++I) {
  457.             AN = FastMath.exp(ALN[I]);
  458.             SUMNM += AN * AMW[I];
  459.         }

  461.         rho = SUMNM / AVOGAD;

  463.         // Compute the high altitude exospheric density correction factor
  464.         double FEX = 1.0;
  465.         if ((scaledSatAlt >= 1000.0) && (scaledSatAlt < 1500.0)) {
  466.             final double ZETA   = (scaledSatAlt - 1000.) * 0.002;
  467.             final double ZETA2  =  ZETA * ZETA;
  468.             final double ZETA3  =  ZETA * ZETA2;
  469.             final double F15C   = CHT[1] + CHT[2] * f10B + CHT[3] * 1500.0 + CHT[4] * f10B * 1500.0;
  470.             final double F15C_ZETA = (CHT[3] + CHT[4] * f10B) * 500.0;
  471.             final double FEX2   = 3.0 * F15C - F15C_ZETA - 3.0;
  472.             final double FEX3   = F15C_ZETA - 2.0 * F15C + 2.0;
  473.             FEX    = 1.0 + FEX2 * ZETA2 + FEX3 * ZETA3;
  474.         }
  475.         if (scaledSatAlt >= 1500.0) {
  476.             FEX    = CHT[1] + CHT[2] * f10B + CHT[3] * scaledSatAlt + CHT[4] * f10B * scaledSatAlt;
  477.         }

  479.         // Apply the exospheric density correction factor.
  480.         rho  *= FEX;

  482.         return rho;

  484.     }

  486.     /** Compute daily temperature correction for Jacchia-Bowman model.
  487.      * @param f10 solar flux index
  488.      * @param solTimeHour local solar time (hours 0-23.999)
  489.      * @param satLat sat lat (radians)
  490.      * @param satAlt height (km)
  491.      * @return dTc correction
  492.      */
  493.     private static double dTc(final double f10, final double solTimeHour,
  494.                               final double satLat, final double satAlt) {

  496.         double dTc = 0;
  497.         final double tx  = solTimeHour / 24.0;
  498.         final double tx2 = tx * tx;
  499.         final double tx3 = tx2 * tx;
  500.         final double tx4 = tx3 * tx;
  501.         final double tx5 = tx4 * tx;
  502.         final double ycs = FastMath.cos(satLat);
  503.         final double f   = (f10 - 100.0) / 100.0;
  504.         double h;
  505.         double sum;

  507.         // Calculates dTc
  508.         if ((satAlt >= 120) && (satAlt <= 200)) {
  509.             final double DTC200 = CDT_SUB[17] + CDT_SUB[18] * tx * ycs + CDT_SUB[19] * tx2 * ycs +
  510.                                   CDT_SUB[20] * tx3 * ycs + CDT_SUB[21] * f * ycs + CDT_SUB[22] * tx * f * ycs +
  511.                                   CDT_SUB[23] * tx2 * f * ycs;
  512.             sum = CDT_SUB[1] + BDT_SUB[2] * f + CDT_SUB[3] * tx * f + CDT_SUB[4] * tx2 * f +
  513.                   CDT_SUB[5] * tx3 * f + CDT_SUB[6] * tx4 * f + CDT_SUB[7] * tx5 * f +
  514.                   CDT_SUB[8] * tx * ycs + CDT_SUB[9] * tx2 * ycs + CDT_SUB[10] * tx3 * ycs +
  515.                   CDT_SUB[11] * tx4 * ycs + CDT_SUB[12] * tx5 * ycs + CDT_SUB[13] * ycs +
  516.                   CDT_SUB[14] * f * ycs + CDT_SUB[15] * tx * f * ycs  + CDT_SUB[16] * tx2 * f * ycs;
  517.             final double DTC200DZ = sum;
  518.             final double CC  = 3.0 * DTC200 - DTC200DZ;
  519.             final double DD  = DTC200 - CC;
  520.             final double ZP  = (satAlt - 120.0) / 80.0;
  521.             dTc = CC * ZP * ZP + DD * ZP * ZP * ZP;
  522.         }

  524.         if ((satAlt > 200.0) && (satAlt <= 240.0)) {
  525.             h = (satAlt - 200.0) / 50.0;
  526.             sum = CDT_SUB[1] * h + BDT_SUB[2] * f * h + CDT_SUB[3] * tx * f * h + CDT_SUB[4] * tx2 * f * h +
  527.                   CDT_SUB[5] * tx3 * f * h + CDT_SUB[6] * tx4 * f * h + CDT_SUB[7] * tx5 * f * h +
  528.                   CDT_SUB[8] * tx * ycs * h + CDT_SUB[9] * tx2 * ycs * h + CDT_SUB[10] * tx3 * ycs * h +
  529.                   CDT_SUB[11] * tx4 * ycs * h + CDT_SUB[12] * tx5 * ycs * h + CDT_SUB[13] * ycs * h +
  530.                   CDT_SUB[14] * f * ycs * h + CDT_SUB[15] * tx * f * ycs * h  + CDT_SUB[16] * tx2 * f * ycs * h +
  531.                   CDT_SUB[17] + CDT_SUB[18] * tx * ycs + CDT_SUB[19] * tx2 * ycs +
  532.                   CDT_SUB[20] * tx3 * ycs + CDT_SUB[21] * f * ycs + CDT_SUB[22] * tx * f * ycs +
  533.                   CDT_SUB[23] * tx2 * f * ycs;
  534.             dTc = sum;
  535.         }

  537.         if ((satAlt > 240.0) && (satAlt <= 300.0)) {
  538.             h = 40.0 / 50.0;
  539.             sum = CDT_SUB[1] * h + BDT_SUB[2] * f * h + CDT_SUB[3] * tx * f * h + CDT_SUB[4] * tx2 * f * h +
  540.                   CDT_SUB[5] * tx3 * f * h + CDT_SUB[6] * tx4 * f * h + CDT_SUB[7] * tx5 * f * h +
  541.                   CDT_SUB[8] * tx * ycs * h + CDT_SUB[9] * tx2 * ycs * h + CDT_SUB[10] * tx3 * ycs * h +
  542.                   CDT_SUB[11] * tx4 * ycs * h + CDT_SUB[12] * tx5 * ycs * h + CDT_SUB[13] * ycs * h +
  543.                   CDT_SUB[14] * f * ycs * h + CDT_SUB[15] * tx * f * ycs * h  + CDT_SUB[16] * tx2 * f * ycs * h +
  544.                   CDT_SUB[17] + CDT_SUB[18] * tx * ycs + CDT_SUB[19] * tx2 * ycs +
  545.                   CDT_SUB[20] * tx3 * ycs + CDT_SUB[21] * f * ycs + CDT_SUB[22] * tx * f * ycs +
  546.                   CDT_SUB[23] * tx2 * f * ycs;
  547.             final double AA = sum;
  548.             final double BB = CDT_SUB[1] + BDT_SUB[2] * f + CDT_SUB[3] * tx * f + CDT_SUB[4] * tx2 * f +
  549.                         CDT_SUB[5] * tx3 * f + CDT_SUB[6] * tx4 * f + CDT_SUB[7] * tx5 * f +
  550.                         CDT_SUB[8] * tx * ycs + CDT_SUB[9] * tx2 * ycs + CDT_SUB[10] * tx3 * ycs +
  551.                         CDT_SUB[11] * tx4 * ycs + CDT_SUB[12] * tx5 * ycs + CDT_SUB[13] * ycs +
  552.                         CDT_SUB[14] * f * ycs + CDT_SUB[15] * tx * f * ycs + CDT_SUB[16] * tx2 * f * ycs;
  553.             h   = 300.0 / 100.0;
  554.             sum = BDT_SUB[1] + BDT_SUB[2] * f  + BDT_SUB[3] * tx * f         + BDT_SUB[4] * tx2 * f +
  555.                   BDT_SUB[5] * tx3 * f      + BDT_SUB[6] * tx4 * f      + BDT_SUB[7] * tx5 * f +
  556.                   BDT_SUB[8] * tx * ycs       + BDT_SUB[9] * tx2 * ycs    + BDT_SUB[10] * tx3 * ycs +
  557.                   BDT_SUB[11] * tx4 * ycs   + BDT_SUB[12] * tx5 * ycs   + BDT_SUB[13] * h * ycs +
  558.                   BDT_SUB[14] * tx * h * ycs    + BDT_SUB[15] * tx2 * h * ycs + BDT_SUB[16] * tx3 * h * ycs +
  559.                   BDT_SUB[17] * tx4 * h * ycs + BDT_SUB[18] * tx5 * h * ycs + BDT_SUB[19] * ycs;
  560.             final double DTC300 = sum;
  561.             sum = BDT_SUB[13] * ycs +
  562.                   BDT_SUB[14] * tx * ycs + BDT_SUB[15] * tx2 * ycs + BDT_SUB[16] * tx3 * ycs +
  563.                   BDT_SUB[17] * tx4 * ycs + BDT_SUB[18] * tx5 * ycs;
  564.             final double DTC300DZ = sum;
  565.             final double CC = 3.0 * DTC300 - DTC300DZ - 3.0 * AA - 2.0 * BB;
  566.             final double  DD = DTC300 - AA - BB - CC;
  567.             final double ZP  = (satAlt - 240.0) / 60.0;
  568.             dTc = AA + BB * ZP + CC * ZP * ZP + DD * ZP * ZP * ZP;
  569.         }

  571.         if ((satAlt > 300.0) && (satAlt <= 600.0)) {
  572.             h   = satAlt / 100.0;
  573.             sum = BDT_SUB[1]    + BDT_SUB[2] * f  + BDT_SUB[3] * tx * f         + BDT_SUB[4] * tx2 * f +
  574.                   BDT_SUB[5] * tx3 * f      + BDT_SUB[6] * tx4 * f      + BDT_SUB[7] * tx5 * f +
  575.                   BDT_SUB[8] * tx * ycs       + BDT_SUB[9] * tx2 * ycs    + BDT_SUB[10] * tx3 * ycs +
  576.                   BDT_SUB[11] * tx4 * ycs   + BDT_SUB[12] * tx5 * ycs   + BDT_SUB[13] * h * ycs +
  577.                   BDT_SUB[14] * tx * h * ycs    + BDT_SUB[15] * tx2 * h * ycs + BDT_SUB[16] * tx3 * h * ycs +
  578.                   BDT_SUB[17] * tx4 * h * ycs + BDT_SUB[18] * tx5 * h * ycs + BDT_SUB[19] * ycs;
  579.             dTc = sum;
  580.         }

  582.         if ((satAlt > 600.0) && (satAlt <= 800.0)) {
  583.             final double ZP = (satAlt - 600.0) / 100.0;
  584.             final double HP = 600.0 / 100.0;
  585.             final double AA  = BDT_SUB[1]    + BDT_SUB[2] * f  + BDT_SUB[3] * tx * f         + BDT_SUB[4] * tx2 * f +
  586.                                BDT_SUB[5] * tx3 * f      + BDT_SUB[6] * tx4 * f      + BDT_SUB[7] * tx5 * f +
  587.                                BDT_SUB[8] * tx * ycs       + BDT_SUB[9] * tx2 * ycs    + BDT_SUB[10] * tx3 * ycs +
  588.                                BDT_SUB[11] * tx4 * ycs   + BDT_SUB[12] * tx5 * ycs   + BDT_SUB[13] * HP * ycs +
  589.                                BDT_SUB[14] * tx * HP * ycs   + BDT_SUB[15] * tx2 * HP * ycs + BDT_SUB[16] * tx3 * HP * ycs +
  590.                                BDT_SUB[17] * tx4 * HP * ycs + BDT_SUB[18] * tx5 * HP * ycs + BDT_SUB[19] * ycs;
  591.             final double BB  = BDT_SUB[13] * ycs +
  592.                                BDT_SUB[14] * tx * ycs    + BDT_SUB[15] * tx2 * ycs + BDT_SUB[16] * tx3 * ycs +
  593.                                BDT_SUB[17] * tx4 * ycs + BDT_SUB[18] * tx5 * ycs;
  594.             final double CC  = -(3.0 * AA + 4.0 * BB) / 4.0;
  595.             final double DD  = (AA + BB) / 4.0;
  596.             dTc = AA + BB * ZP + CC * ZP * ZP + DD * ZP * ZP * ZP;
  597.         }

  599.         return dTc;
  600.     }

  602.     /** Evaluates Equation (1).
  603.      * @param z altitude
  604.      * @return equation (1) value
  605.      */
  606.     private static double xAmbar(final double z) {
  607.         final double dz = z - 100.;
  608.         double amb = CXAMB[7];
  609.         for (int i = 6; i >= 1; --i) {
  610.             amb = dz * amb + CXAMB[i];
  611.         }
  612.         return amb;
  613.     }

  615.     /**  Evaluates Equation (10) or Equation (13), depending on Z.
  616.      * @param z altitude
  617.      * @param TC tc array ???
  618.      * @return equation (10) value
  619.      */
  620.     private static double xLocal(final double z, final double[] TC) {
  621.         final double dz = z - 125;
  622.         if (dz <= 0) {
  623.             return ((-9.8204695e-6 * dz - 7.3039742e-4) * dz * dz + 1.0) * dz * TC[1] + TC[0];
  624.         } else {
  625.             return TC[0] + TC[2] * FastMath.atan(TC[3] * dz * (1 + 4.5e-6 * FastMath.pow(dz, 2.5)));
  626.         }
  627.     }

  629.     /** Evaluates Equation (8) of gravity field.
  630.      * @param z altitude
  631.      * @return the gravity field
  632.      */
  633.     private static double xGrav(final double z) {
  634.         final double temp = 1.0 + z / 6356.766;
  635.         return 9.80665 / (temp * temp);
  636.     }

  638.     /** Compute semi-annual variation (delta log(rho)).
  639.      * @param day day of year
  640.      * @param height height (km)
  641.      * @param f10Bar average 81-day centered f10
  642.      * @return semi-annual variation
  643.      */
  644.     private static double semian (final double day, final double height, final double f10Bar) {

  646.         final double f10Bar2 = f10Bar * f10Bar;
  647.         final double htz = height / 1000.0;

  649.         // SEMIANNUAL AMPLITUDE
  650.         final double fzz = FZM[1] + FZM[2] * f10Bar  + FZM[3] * f10Bar * htz +
  651.                            FZM[4] * f10Bar * htz * htz + FZM[5] * f10Bar * f10Bar * htz +
  652.                            FZM[6] * f10Bar * f10Bar * htz * htz;

  654.         // SEMIANNUAL PHASE FUNCTION
  655.         final double tau   = TWOPI * (day - 1.0) / 365;
  656.         final double sin1P = FastMath.sin(tau);
  657.         final double cos1P = FastMath.cos(tau);
  658.         final double sin2P = FastMath.sin(2.0 * tau);
  659.         final double cos2P = FastMath.cos(2.0 * tau);
  660.         final double sin3P = FastMath.sin(3.0 * tau);
  661.         final double cos3P = FastMath.cos(3.0 * tau);
  662.         final double sin4P = FastMath.sin(4.0 * tau);
  663.         final double cos4P = FastMath.cos(4.0 * tau);
  664.         final double gtz = GTM[1] + GTM[2] * sin1P + GTM[3] * cos1P +
  665.                            GTM[4] * sin2P + GTM[5] * cos2P +
  666.                            GTM[6] * sin3P + GTM[7] * cos3P +
  667.                            GTM[8] * sin4P + GTM[9] * cos4P +
  668.                            GTM[10] * f10Bar + GTM[11] * f10Bar * sin1P + GTM[12] * f10Bar * cos1P +
  669.                            GTM[13] * f10Bar * sin2P + GTM[14] * f10Bar * cos2P +
  670.                            GTM[15] * f10Bar * sin3P + GTM[16] * f10Bar * cos3P +
  671.                            GTM[17] * f10Bar * sin4P + GTM[18] * f10Bar * cos4P +
  672.                            GTM[19] * f10Bar2  + GTM[20] * f10Bar2  * sin1P + GTM[21] * f10Bar2  * cos1P +
  673.                            GTM[22] * f10Bar2  * sin2P + GTM[23] * f10Bar2  * cos2P;

  675.         return FastMath.max(1.0e-6, fzz) * gtz;

  677.     }

  679.     /** Compute day of year.
  680.      * @param d1950 (days since 1950)
  681.      * @return the number days in year
  682.      */
  683.     private static double dayOfYear(final double d1950) {

  685.         int iyday = (int) d1950;
  686.         final double frac = d1950 - iyday;
  687.         iyday = iyday + 364;

  689.         int itemp = iyday / 1461;

  691.         iyday = iyday - itemp * 1461;
  692.         itemp = iyday / 365;
  693.         if (itemp >= 3) {
  694.             itemp = 3;
  695.         }
  696.         iyday = iyday - 365 * itemp + 1;
  697.         return iyday + frac;
  698.     }

  700.     // OUTPUT:

  702.     /** Get the exospheric temperature above input position.
  703.      * {@link #getDensity(double, double, double, double, double, double, double, double, double, double, double, double, double)}
  704.      * <b> must </b> must be called before calling this function.
  705.      * @return the exospheric temperature (deg K)
  706.      */
  707.     public double getExosphericTemp() {
  708.         return temp[1];
  709.     }

  711.     /** Get the temperature at input position.
  712.      * {@link #getDensity(double, double, double, double, double, double, double, double, double, double, double, double, double)}
  713.      * <b> must </b> must be called before calling this function.
  714.      * @return the local temperature (deg K)
  715.      */
  716.     public double getLocalTemp() {
  717.         return temp[2];
  718.     }

  720.     /** Get the local density.
  721.      * @param date current date
  722.      * @param position current position in frame
  723.      * @param frame the frame in which is defined the position
  724.      * @return local density (kg/m<sup>3</sup>)
  725.      * @exception OrekitException if date is out of range of solar activity
  726.      */
  727.     public double getDensity(final AbsoluteDate date, final Vector3D position,
  728.                              final Frame frame)
  729.         throws OrekitException {
  730.         // check if data are available :
  731.         if (date.compareTo(inputParams.getMaxDate()) > 0 ||
  732.             date.compareTo(inputParams.getMinDate()) < 0) {
  733.             throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE,
  734.                                       date, inputParams.getMinDate(), inputParams.getMaxDate());
  735.         }

  737.         // compute modified julian days date
  738.         final double dateMJD = date.durationFrom(AbsoluteDate.MODIFIED_JULIAN_EPOCH) / Constants.JULIAN_DAY;

  740.         // compute geodetic position
  741.         final GeodeticPoint inBody = earth.transform(position, frame, date);

  743.         // compute sun position
  744.         final GeodeticPoint sunInBody =
  745.             earth.transform(sun.getPVCoordinates(date, frame).getPosition(), frame, date);
  746.         return getDensity(dateMJD,
  747.                           sunInBody.getLongitude(), sunInBody.getLatitude(),
  748.                           inBody.getLongitude(), inBody.getLatitude(),
  749.                           inBody.getAltitude(), inputParams.getF10(date),
  750.                           inputParams.getF10B(date),
  751.                           inputParams.getAp(date), inputParams.getS10(date),
  752.                           inputParams.getS10B(date), inputParams.getXM10(date),
  753.                           inputParams.getXM10B(date));
  754.     }

  756.     /** Get the inertial velocity of atmosphere molecules.
  757.      * Here the case is simplified : atmosphere is supposed to have a null velocity
  758.      * in earth frame.
  759.      * @param date current date
  760.      * @param position current position in frame
  761.      * @param frame the frame in which is defined the position
  762.      * @return velocity (m/s) (defined in the same frame as the position)
  763.      * @exception OrekitException if some frame conversion cannot be performed
  764.      */
  765.     public Vector3D getVelocity(final AbsoluteDate date, final Vector3D position,
  766.                                 final Frame frame)
  767.         throws OrekitException {
  768.         final Transform bodyToFrame = earth.getBodyFrame().getTransformTo(frame, date);
  769.         final Vector3D posInBody = bodyToFrame.getInverse().transformPosition(position);
  770.         final PVCoordinates pvBody = new PVCoordinates(posInBody, new Vector3D(0, 0, 0));
  771.         final PVCoordinates pvFrame = bodyToFrame.transformPVCoordinates(pvBody);
  772.         return pvFrame.getVelocity();
  773.     }

  775. }
Created with JaCoCo 0.6.3.201306030806