JB2006.java

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

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.util.FastMath;
  21. import org.hipparchus.util.FieldSinCos;
  22. import org.hipparchus.util.MathUtils;
  23. import org.orekit.annotation.DefaultDataContext;
  24. import org.orekit.bodies.BodyShape;
  25. import org.orekit.data.DataContext;
  26. import org.orekit.time.AbsoluteDate;
  27. import org.orekit.time.TimeScale;
  28. import org.orekit.utils.ExtendedPositionProvider;

  30. /**
  31.  * This is the realization of the Jacchia-Bowman 2006 atmospheric model.
  32.  * <p>
  33.  * It is described in the paper: <br>
  34.  *
  35.  * <a href="http://sol.spacenvironment.net/~JB2006/pubs/JB2006_AIAA-6166_model.pdf">A
  36.  * New Empirical Thermospheric Density Model JB2006 Using New Solar Indices</a><br>
  37.  *
  38.  * <i>Bruce R. Bowman, W. Kent Tobiska and Frank A. Marcos</i> <br>
  39.  * <p>
  40.  * AIAA 2006-6166<br>
  41.  * </p>
  42.  *
  43.  * <p>
  44.  * This model provides dense output for all altitudes and positions. Output data are :
  45.  * <ul>
  46.  * <li>Exospheric Temperature above Input Position (deg K)</li>
  47.  * <li>Temperature at Input Position (deg K)</li>
  48.  * <li>Total Mass-Density at Input Position (kg/m³)</li>
  49.  * </ul>
  50.  *
  51.  * <p>
  52.  * The model needs geographical and time information to compute general values,
  53.  * but also needs space weather data : mean and daily solar flux, retrieved threw
  54.  * different indices, and planetary geomagnetic indices. <br>
  55.  * More information on these indices can be found on the  <a
  56.  * href="http://sol.spacenvironment.net/~JB2006/JB2006_index.html">
  57.  * official JB2006 website.</a>
  58.  * </p>
  59.  *
  60.  * @author Bruce R Bowman (HQ AFSPC, Space Analysis Division), Feb 2006: FORTRAN routine
  61.  * @author Fabien Maussion (java translation)
  62.  * @author Bryan Cazabonne (Orekit 13 update and field translation)
  63.  * @since 13.1
  64.  */
  65. public class JB2006 extends AbstractJacchiaBowmanModel {

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

  70.     /** GT global model values (1978-2004 fit). */
  71.     private static final double[] GTM = {-0.833646e+00, -0.265450e+00, 0.467603e+00, -0.299906e+00, -0.105451e+00,
  72.                                          -0.165537e-01, -0.380037e-01, -0.150991e-01, -0.541280e-01,  0.119554e-01,
  73.                                          0.437544e-02, -0.369016e-02, 0.206763e-02, -0.142888e-02, -0.867124e-05,
  74.                                          0.189032e-04, 0.156988e-03, 0.491286e-03, -0.391484e-04, -0.126854e-04,
  75.                                          0.134078e-04, -0.614176e-05, 0.343423e-05};

  77.     /** External data container. */
  78.     private final JB2006InputParameters inputParams;

  80.     /**
  81.      * Constructor with space environment information for internal computation.
  82.      *
  83.      * @param parameters the solar and magnetic activity data
  84.      * @param sun        the sun position
  85.      * @param earth      the earth body shape
  86.      */
  87.     @DefaultDataContext
  88.     public JB2006(final JB2006InputParameters parameters, final ExtendedPositionProvider sun, final BodyShape earth) {
  89.         this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
  90.     }

  92.     /**
  93.      * Constructor with space environment information for internal computation.
  94.      *
  95.      * @param parameters the solar and magnetic activity data
  96.      * @param sun        the sun position
  97.      * @param earth      the earth body shape
  98.      * @param utc        UTC time scale. Used to computed the day fraction.
  99.      */
  100.     public JB2006(final JB2006InputParameters parameters, final ExtendedPositionProvider sun,
  101.                   final BodyShape earth, final TimeScale utc) {
  102.         super(sun, utc, earth, parameters.getMinDate(), parameters.getMaxDate());
  103.         this.inputParams = parameters;
  104.     }

  106.     /** {@inheritDoc} */
  107.     @Override
  108.     protected double computeTInf(final AbsoluteDate date, final double tsubl, final double dtclst) {
  109.         return tsubl + getDtg(date) + dtclst;
  110.     }

  112.     /** {@inheritDoc} */
  113.     @Override
  114.     protected <T extends CalculusFieldElement<T>> T computeTInf(final AbsoluteDate date, final T tsubl, final T dtclst) {
  115.         return tsubl.add(getDtg(date)).add(dtclst);
  116.     }

  118.     /** {@inheritDoc} */
  119.     @Override
  120.     protected double computeTc(final AbsoluteDate date) {
  121.         final double f10   = inputParams.getF10(date);
  122.         final double f10B  = inputParams.getF10B(date);
  123.         final double s10   = inputParams.getS10(date);
  124.         final double s10B  = inputParams.getS10B(date);
  125.         final double xm10  = inputParams.getXM10(date);
  126.         final double xm10B = inputParams.getXM10B(date);
  127.         return 379 + 3.353 * f10B + 0.358 * (f10 - f10B) + 2.094 * (s10 - s10B) + 0.343 * (xm10 - xm10B);
  128.     }

  130.     /** {@inheritDoc} */
  131.     @Override
  132.     protected double getF10(final AbsoluteDate date) {
  133.         return inputParams.getF10(date);
  134.     }

  136.     /** {@inheritDoc} */
  137.     @Override
  138.     protected double getF10B(final AbsoluteDate date) {
  139.         return inputParams.getF10B(date);
  140.     }

  142.     /** {@inheritDoc} */
  143.     @Override
  144.     protected double semian(final AbsoluteDate date, final double day, final double height) {

  146.         final double f10Bar = inputParams.getF10B(date);
  147.         final double f10Bar2 = f10Bar * f10Bar;
  148.         final double htz = height / 1000.0;

  150.         // SEMIANNUAL AMPLITUDE
  151.         final double fzz = FZM[0] + FZM[1] * f10Bar + FZM[2] * f10Bar * htz + FZM[3] * f10Bar * htz * htz + FZM[4] * f10Bar * f10Bar * htz + FZM[5] * f10Bar * f10Bar * htz * htz;

  153.         // SEMIANNUAL PHASE FUNCTION
  154.         final double tau = MathUtils.TWO_PI * (day - 1.0) / 365;
  155.         final double sin1P = FastMath.sin(tau);
  156.         final double cos1P = FastMath.cos(tau);
  157.         final double sin2P = FastMath.sin(2.0 * tau);
  158.         final double cos2P = FastMath.cos(2.0 * tau);
  159.         final double sin3P = FastMath.sin(3.0 * tau);
  160.         final double cos3P = FastMath.cos(3.0 * tau);
  161.         final double sin4P = FastMath.sin(4.0 * tau);
  162.         final double cos4P = FastMath.cos(4.0 * tau);
  163.         final double gtz = GTM[0] +
  164.                            GTM[1] * sin1P +
  165.                            GTM[2] * cos1P +
  166.                            GTM[3] * sin2P +
  167.                            GTM[4] * cos2P +
  168.                            GTM[5] * sin3P +
  169.                            GTM[6] * cos3P +
  170.                            GTM[7] * sin4P +
  171.                            GTM[8] * cos4P +
  172.                            GTM[9] * f10Bar +
  173.                            GTM[10] * f10Bar * sin1P +
  174.                            GTM[11] * f10Bar * cos1P +
  175.                            GTM[12] * f10Bar * sin2P +
  176.                            GTM[13] * f10Bar * cos2P +
  177.                            GTM[14] * f10Bar * sin3P +
  178.                            GTM[15] * f10Bar * cos3P +
  179.                            GTM[16] * f10Bar * sin4P +
  180.                            GTM[17] * f10Bar * cos4P +
  181.                            GTM[18] * f10Bar2 +
  182.                            GTM[19] * f10Bar2 * sin1P +
  183.                            GTM[20] * f10Bar2 * cos1P +
  184.                            GTM[21] * f10Bar2 * sin2P +
  185.                            GTM[22] * f10Bar2 * cos2P;

  187.         return FastMath.max(1.0e-6, fzz) * gtz;
  188.     }
  189.     /** {@inheritDoc} */
  190.     @Override
  191.     protected <T extends CalculusFieldElement<T>> T semian(final AbsoluteDate date, final T doy, final T height) {

  193.         final double f10Bar = getF10B(date);
  194.         final double f10Bar2 = f10Bar * f10Bar;
  195.         final T htz = height.divide(1000.0);

  197.         // SEMIANNUAL AMPLITUDE
  198.         final T fzz = htz.multiply(FZM[2] * f10Bar).add(htz.square().multiply(FZM[3] * f10Bar)).add(htz.multiply(FZM[4] * f10Bar * f10Bar)).add(htz.square().multiply(FZM[5] * f10Bar * f10Bar)).add(FZM[0] + FZM[1] * f10Bar);

  200.         // SEMIANNUAL PHASE FUNCTION
  201.         final T tau   = doy.subtract(1).divide(365).multiply(MathUtils.TWO_PI);
  202.         final FieldSinCos<T> sc1P = FastMath.sinCos(tau);
  203.         final FieldSinCos<T> sc2P = FastMath.sinCos(tau.multiply(2.0));
  204.         final FieldSinCos<T> sc3P = FastMath.sinCos(tau.multiply(3.0));
  205.         final FieldSinCos<T> sc4P = FastMath.sinCos(tau.multiply(4.0));
  206.         final T gtz = sc1P.sin().multiply(GTM[1]).add(
  207.                       sc1P.cos().multiply(GTM[2])).add(
  208.                       sc2P.sin().multiply(GTM[3])).add(
  209.                       sc2P.cos().multiply(GTM[4])).add(
  210.                       sc3P.sin().multiply(GTM[5])).add(
  211.                       sc3P.cos().multiply(GTM[6])).add(
  212.                       sc4P.sin().multiply(GTM[7])).add(
  213.                       sc4P.cos().multiply(GTM[8])).add(
  214.                       GTM[9] * f10Bar).add(
  215.                       sc1P.sin().multiply(f10Bar).multiply(GTM[10])).add(
  216.                       sc1P.cos().multiply(f10Bar).multiply(GTM[11])).add(
  217.                       sc2P.sin().multiply(f10Bar).multiply(GTM[12])).add(
  218.                       sc2P.cos().multiply(f10Bar).multiply(GTM[13])).add(
  219.                       sc3P.sin().multiply(f10Bar).multiply(GTM[14])).add(
  220.                       sc3P.cos().multiply(f10Bar).multiply(GTM[15])).add(
  221.                       sc4P.sin().multiply(f10Bar).multiply(GTM[16])).add(
  222.                       sc4P.cos().multiply(f10Bar).multiply(GTM[17])).add(
  223.                       GTM[18] * f10Bar2).add(
  224.                       sc1P.sin().multiply(f10Bar2).multiply(GTM[19])).add(
  225.                       sc1P.cos().multiply(f10Bar2).multiply(GTM[20])).add(
  226.                       sc2P.sin().multiply(f10Bar2).multiply(GTM[21])).add(
  227.                       sc2P.cos().multiply(f10Bar2).multiply(GTM[22])).add(GTM[0]);

  229.         return fzz.getReal() > 1.0e-6 ? gtz.multiply(fzz) : gtz.multiply(1.0e-6);
  230.     }

  232.     /** Computes the temperature computed by Equation (18).
  233.      * @param date computation epoch
  234.      * @return the temperature given by Equation (18)
  235.      */
  236.     private double getDtg(final AbsoluteDate date) {
  237.         // Equation (18)
  238.         final double ap = inputParams.getAp(date);
  239.         final double expAp = FastMath.exp(-0.08 * ap);
  240.         return ap + 100. * (1. - expAp);
  241.     }
  242. }
Created with JaCoCo 0.8.12.202403310830