JB2008.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.hipparchus.util.SinCos;
  24. import org.orekit.annotation.DefaultDataContext;
  25. import org.orekit.bodies.BodyShape;
  26. import org.orekit.data.DataContext;
  27. import org.orekit.time.AbsoluteDate;
  28. import org.orekit.time.FieldAbsoluteDate;
  29. import org.orekit.time.TimeScale;

  31. import org.orekit.utils.Constants;
  32. import org.orekit.utils.ExtendedPositionProvider;

  34. /** This is the realization of the Jacchia-Bowman 2008 atmospheric model.
  35.  * <p>
  36.  * It is described in the paper:<br>
  37.  * <a href="https://www.researchgate.net/publication/228621668_A_New_Empirical_Thermospheric_Density_Model_JB2008_Using_New_Solar_and_Geomagnetic_Indices">A
  38.  * New Empirical Thermospheric Density Model JB2008 Using New Solar Indices</a><br>
  39.  * <i>Bruce R. Bowman &amp; al.</i><br>
  40.  * AIAA 2008-6438<br>
  41.  * </p>
  42.  * <p>
  43.  * Two computation methods are proposed to the user:
  44.  * <ul>
  45.  * <li>one OREKIT independent and compliant with initial FORTRAN routine entry values:
  46.  *     {@link #getDensity(double, double, double, double, double, double, double, double, double, double, double, double, double, double, double)}. </li>
  47.  * <li>one compliant with OREKIT Atmosphere interface, necessary to the
  48.  *     {@link org.orekit.forces.drag.DragForce drag force model} computation.</li>
  49.  * </ul>
  50.  * <p>
  51.  * This model provides dense output for all altitudes and positions. Output data are :
  52.  * <ul>
  53.  * <li>Exospheric Temperature above Input Position (deg K)</li>
  54.  * <li>Temperature at Input Position (deg K)</li>
  55.  * <li>Total Mass-Density at Input Position (kg/m³)</li>
  56.  * </ul>
  57.  * <p>
  58.  * The model needs geographical and time information to compute general values,
  59.  * but also needs space weather data : mean and daily solar flux, retrieved through
  60.  * different indices, and planetary geomagnetic indices.<br>
  61.  * More information on these indices can be found on the <a
  62.  * href="http://sol.spacenvironment.net/~JB2008/indices.html">
  63.  * official JB2008 website.</a>
  64.  * </p>
  65.  *
  66.  * @author Bruce R Bowman (HQ AFSPC, Space Analysis Division), 2008: FORTRAN routine
  67.  * @author Pascal Parraud (java translation)
  68.  */
  69. public class JB2008 extends AbstractJacchiaBowmanModel {

  71.     /** FZ global model values (1997-2006 fit).  */
  72.     private static final double[] FZM = {
  73.         0.2689e+00, -0.1176e-01, 0.2782e-01, -0.2782e-01, 0.3470e-03
  74.     };

  76.     /** GT global model values (1997-2006 fit). */
  77.     private static final double[] GTM = {
  78.         -0.3633e+00, 0.8506e-01,  0.2401e+00, -0.1897e+00, -0.2554e+00,
  79.         -0.1790e-01, 0.5650e-03, -0.6407e-03, -0.3418e-02, -0.1252e-02
  80.     };

  82.     /** External data container. */
  83.     private final JB2008InputParameters inputParams;

  85.     /** Constructor with space environment information for internal computation.
  86.      *
  87.      * <p>This method uses the {@link DataContext#getDefault() default data context}.
  88.      *
  89.      * @param parameters the solar and magnetic activity data
  90.      * @param sun the sun position
  91.      * @param earth the earth body shape
  92.      * @see #JB2008(JB2008InputParameters, ExtendedPositionProvider, BodyShape, TimeScale)
  93.      */
  94.     @DefaultDataContext
  95.     public JB2008(final JB2008InputParameters parameters,
  96.                   final ExtendedPositionProvider sun, final BodyShape earth) {
  97.         this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
  98.     }

  100.     /**
  101.      * Constructor with space environment information for internal computation.
  102.      *
  103.      * @param parameters the solar and magnetic activity data
  104.      * @param sun        the sun position
  105.      * @param earth      the earth body shape
  106.      * @param utc        UTC time scale. Used to computed the day fraction.
  107.      * @since 10.1
  108.      */
  109.     public JB2008(final JB2008InputParameters parameters,
  110.                   final ExtendedPositionProvider sun,
  111.                   final BodyShape earth,
  112.                   final TimeScale utc) {
  113.         super(sun, utc, earth,
  114.               parameters == null ? AbsoluteDate.PAST_INFINITY : parameters.getMinDate(),
  115.               parameters == null ? AbsoluteDate.FUTURE_INFINITY : parameters.getMaxDate());
  116.         this.inputParams = parameters;
  117.     }

  119.     /** Get the local density with initial entries.
  120.      * <p>
  121.      * The method creates a new instance of {@link JB2008} model and set the input
  122.      * solar activity data equal to the provided one. These data are then
  123.      * available form {@link AbsoluteDate#PAST_INFINITY} to {@link AbsoluteDate#FUTURE_INFINITY}.
  124.      * </p>
  125.      * @param dateMJD date and time, in modified julian days and fraction
  126.      * @param sunRA Right Ascension of Sun (radians).
  127.      * @param sunDecli Declination of Sun (radians).
  128.      * @param satLon Right Ascension of position (radians)
  129.      * @param satLat Geocentric latitude of position (radians)
  130.      * @param satAlt Height of position (m)
  131.      * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
  132.      *        (Tabular time 1.0 day earlier)
  133.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
  134.      *        (Tabular time 1.0 day earlier)
  135.      * @param s10 EUV index (26-34 nm) scaled to F10<br>
  136.      *        (Tabular time 1 day earlier)
  137.      * @param s10B UV 81-day averaged centered index
  138.      *        (Tabular time 1 day earlier)
  139.      * @param xm10 MG2 index scaled to F10<br>
  140.      *        (Tabular time 2.0 days earlier)
  141.      * @param xm10B MG2 81-day ave. centered index<br>
  142.      *        (Tabular time 2.0 days earlier)
  143.      * @param y10 Solar X-Ray &amp; Lya index scaled to F10<br>
  144.      *        (Tabular time 5.0 days earlier)
  145.      * @param y10B Solar X-Ray &amp; Lya 81-day ave. centered index<br>
  146.      *        (Tabular time 5.0 days earlier)
  147.      * @param dstdtc Temperature change computed from Dst index
  148.      * @return total mass-Density at input position (kg/m³)
  149.      */
  150.     public double getDensity(final double dateMJD, final double sunRA, final double sunDecli,
  151.                              final double satLon, final double satLat, final double satAlt,
  152.                              final double f10, final double f10B, final double s10,
  153.                              final double s10B, final double xm10, final double xm10B,
  154.                              final double y10, final double y10B, final double dstdtc) {
  155.         final LocalProvider provider = new LocalProvider(f10, f10B, s10, s10B, xm10, xm10B, y10, y10B, dstdtc);
  156.         final JB2008 modelWithLocalData = new JB2008(provider, getSun(), getEarth(), getUtc());
  157.         final double mjd = FastMath.floor(dateMJD);
  158.         final AbsoluteDate computationEpoch = AbsoluteDate.createMJDDate((int) mjd, (dateMJD - mjd) * Constants.JULIAN_DAY, getUtc());
  159.         return modelWithLocalData.computeDensity(computationEpoch, sunRA, sunDecli, satLon, satLat, satAlt);
  160.     }

  162.     /** Get the local density with initial entries.
  163.      * <p>
  164.      * The method creates a new instance of {@link JB2008} model and set the input
  165.      * solar activity data equal to the provided one. These data are then
  166.      * available form {@link AbsoluteDate#PAST_INFINITY} to {@link AbsoluteDate#FUTURE_INFINITY}.
  167.      * </p>
  168.      * @param dateMJD date and time, in modified julian days and fraction
  169.      * @param sunRA Right Ascension of Sun (radians).
  170.      * @param sunDecli Declination of Sun (radians).
  171.      * @param satLon Right Ascension of position (radians)
  172.      * @param satLat Geocentric latitude of position (radians)
  173.      * @param satAlt Height of position (m)
  174.      * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
  175.      *        (Tabular time 1.0 day earlier)
  176.      * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
  177.      *        (Tabular time 1.0 day earlier)
  178.      * @param s10 EUV index (26-34 nm) scaled to F10<br>
  179.      *        (Tabular time 1 day earlier)
  180.      * @param s10B UV 81-day averaged centered index
  181.      *        (Tabular time 1 day earlier)
  182.      * @param xm10 MG2 index scaled to F10<br>
  183.      *        (Tabular time 2.0 days earlier)
  184.      * @param xm10B MG2 81-day ave. centered index<br>
  185.      *        (Tabular time 2.0 days earlier)
  186.      * @param y10 Solar X-Ray &amp; Lya index scaled to F10<br>
  187.      *        (Tabular time 5.0 days earlier)
  188.      * @param y10B Solar X-Ray &amp; Lya 81-day ave. centered index<br>
  189.      *        (Tabular time 5.0 days earlier)
  190.      * @param dstdtc Temperature change computed from Dst index
  191.      * @param <T> type of the elements
  192.      * @return total mass-Density at input position (kg/m³)
  193.      */
  194.     public <T extends CalculusFieldElement<T>> T getDensity(final T dateMJD, final T sunRA, final T sunDecli,
  195.                                                             final T satLon, final T satLat, final T satAlt,
  196.                                                             final double f10, final double f10B, final double s10,
  197.                                                             final double s10B, final double xm10, final double xm10B,
  198.                                                             final double y10, final double y10B, final double dstdtc) {
  199.         final LocalProvider provider = new LocalProvider(f10, f10B, s10, s10B, xm10, xm10B, y10, y10B, dstdtc);
  200.         final JB2008 modelWithLocalData = new JB2008(provider, getSun(), getEarth(), getUtc());
  201.         final T mjd = FastMath.floor(dateMJD);
  202.         final FieldAbsoluteDate<T> computationEpoch = FieldAbsoluteDate.createMJDDate((int) mjd.getReal(), dateMJD.subtract(mjd).multiply(Constants.JULIAN_DAY), getUtc());
  203.         return modelWithLocalData.computeDensity(computationEpoch, sunRA, sunDecli, satLon, satLat, satAlt);
  204.     }

  206.         /** {@inheritDoc} */
  207.     @Override
  208.     protected double computeTInf(final AbsoluteDate date, final double tsubl, final double dtclst) {
  209.         return tsubl + inputParams.getDSTDTC(date) + dtclst;
  210.     }

  212.     /** {@inheritDoc} */
  213.     @Override
  214.     protected double computeTc(final AbsoluteDate date) {
  215.         final double f10   = inputParams.getF10(date);
  216.         final double f10B  = inputParams.getF10B(date);
  217.         final double s10   = inputParams.getS10(date);
  218.         final double s10B  = inputParams.getS10B(date);
  219.         final double xm10  = inputParams.getXM10(date);
  220.         final double xm10B = inputParams.getXM10B(date);
  221.         final double y10   = inputParams.getY10(date);
  222.         final double y10B  = inputParams.getY10B(date);
  223.         final double fn    = FastMath.min(1.0, FastMath.pow(f10B / 240., 0.25));
  224.         final double fsb   = f10B * fn + s10B * (1. - fn);
  225.         return 392.4 + 3.227 * fsb + 0.298 * (f10 - f10B) + 2.259 * (s10 - s10B) + 0.312 * (xm10 - xm10B) + 0.178 * (y10 - y10B);
  226.     }

  228.     /** {@inheritDoc} */
  229.     @Override
  230.     protected double getF10(final AbsoluteDate date) {
  231.         return inputParams.getF10(date);
  232.     }

  234.     /** {@inheritDoc} */
  235.     @Override
  236.     protected double getF10B(final AbsoluteDate date) {
  237.         return inputParams.getF10B(date);
  238.     }

  240.     /** {@inheritDoc} */
  241.     @Override
  242.     protected <T extends CalculusFieldElement<T>> T computeTInf(final AbsoluteDate date, final T tsubl, final T dtclst) {
  243.         return tsubl.add(inputParams.getDSTDTC(date)).add(dtclst);
  244.     }

  246.     /** {@inheritDoc} */
  247.     @Override
  248.     protected double semian(final AbsoluteDate date, final double day, final double altKm) {

  250.         final double f10B = inputParams.getF10B(date);
  251.         final double s10B = inputParams.getS10B(date);
  252.         final double xm10B = inputParams.getXM10B(date);

  254.         final double htz = altKm / 1000.0;

  256.         // COMPUTE NEW 81-DAY CENTERED SOLAR INDEX FOR FZ
  257.         double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;

  259.         // SEMIANNUAL AMPLITUDE
  260.         final double fzz = FZM[0] + fsmb * (FZM[1] + htz * (FZM[2] + FZM[3] * htz + FZM[4] * fsmb));

  262.         // COMPUTE DAILY 81-DAY CENTERED SOLAR INDEX FOR GT
  263.         fsmb  = f10B - 0.75 * s10B - 0.37 * xm10B;

  265.         // SEMIANNUAL PHASE FUNCTION
  266.         final double tau   = MathUtils.TWO_PI * (day - 1.0) / 365;
  267.         final SinCos sc1P = FastMath.sinCos(tau);
  268.         final SinCos sc2P = SinCos.sum(sc1P, sc1P);
  269.         final double gtz = GTM[0] + GTM[1] * sc1P.sin() + GTM[2] * sc1P.cos() + GTM[3] * sc2P.sin() + GTM[4] * sc2P.cos() +
  270.                    fsmb * (GTM[5] + GTM[6] * sc1P.sin() + GTM[7] * sc1P.cos() + GTM[8] * sc2P.sin() + GTM[9] * sc2P.cos());

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

  274.     }

  276.     /** {@inheritDoc} */
  277.     @Override
  278.     protected <T extends CalculusFieldElement<T>>  T semian(final AbsoluteDate date, final T day, final T altKm) {

  280.         final double f10B = inputParams.getF10B(date);
  281.         final double s10B = inputParams.getS10B(date);
  282.         final double xm10B = inputParams.getXM10B(date);

  284.         final T htz = altKm.divide(1000.0);

  286.         // COMPUTE NEW 81-DAY CENTERED SOLAR INDEX FOR FZ
  287.         double fsmb = f10B - 0.70 * s10B - 0.04 * xm10B;

  289.         // SEMIANNUAL AMPLITUDE
  290.         final T fzz = htz.multiply(FZM[3]).add(FZM[2] + FZM[4] * fsmb).multiply(htz).add(FZM[1]).multiply(fsmb).add(FZM[0]);

  292.         // COMPUTE DAILY 81-DAY CENTERED SOLAR INDEX FOR GT
  293.         fsmb  = f10B - 0.75 * s10B - 0.37 * xm10B;

  295.         // SEMIANNUAL PHASE FUNCTION
  296.         final T tau   = day.subtract(1).divide(365).multiply(MathUtils.TWO_PI);
  297.         final FieldSinCos<T> sc1P = FastMath.sinCos(tau);
  298.         final FieldSinCos<T> sc2P = FieldSinCos.sum(sc1P, sc1P);
  299.         final T gtz =           sc2P.cos().multiply(GTM[9]).
  300.                             add(sc2P.sin().multiply(GTM[8])).
  301.                             add(sc1P.cos().multiply(GTM[7])).
  302.                             add(sc1P.sin().multiply(GTM[6])).
  303.                             add(GTM[5]).multiply(fsmb).
  304.                         add(    sc2P.cos().multiply(GTM[4]).
  305.                             add(sc2P.sin().multiply(GTM[3])).
  306.                             add(sc1P.cos().multiply(GTM[2])).
  307.                             add(sc1P.sin().multiply(GTM[1])).
  308.                             add(GTM[0]));

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

  312.     }

  314.     /** Local provider for solar activity data. */
  315.     private static class LocalProvider implements JB2008InputParameters {

  317.         /** 10.7-cm Solar flux. */
  318.         private double f10;

  320.         /** 10.7-cm Solar Flux, averaged 81-day centered on the input time. */
  321.         private double f10B;

  323.         /** EUV index (26-34 nm) scaled to F10. */
  324.         private double s10;

  326.         /** UV 81-day averaged centered index. */
  327.         private double s10B;

  329.         /** MG2 index scaled to F10. */
  330.         private double xm10;

  332.         /** MG2 81-day ave. centered index. */
  333.         private double xm10B;

  335.         /** Solar X-Ray &amp; Lya index scaled to F10. */
  336.         private double y10;

  338.         /** Solar X-Ray &amp; Lya 81-day ave. centered index. */
  339.         private double y10B;

  341.         /** Temperature change computed from Dst index. */
  342.         private double dstdtc;

  344.         /** Constructor.
  345.          * @param f10 10.7-cm Solar flux (1e<sup>-22</sup>*Watt/(m²*Hertz))<br>
  346.          *        (Tabular time 1.0 day earlier)
  347.          * @param f10B 10.7-cm Solar Flux, averaged 81-day centered on the input time<br>
  348.          *        (Tabular time 1.0 day earlier)
  349.          * @param s10 EUV index (26-34 nm) scaled to F10<br>
  350.          *        (Tabular time 1 day earlier)
  351.          * @param s10B UV 81-day averaged centered index
  352.          *        (Tabular time 1 day earlier)
  353.          * @param xm10 MG2 index scaled to F10<br>
  354.          *        (Tabular time 2.0 days earlier)
  355.          * @param xm10B MG2 81-day ave. centered index<br>
  356.          *        (Tabular time 2.0 days earlier)
  357.          * @param y10 Solar X-Ray &amp; Lya index scaled to F10<br>
  358.          *        (Tabular time 5.0 days earlier)
  359.          * @param y10B Solar X-Ray &amp; Lya 81-day ave. centered index<br>
  360.          *        (Tabular time 5.0 days earlier)
  361.          * @param dstdtc Temperature change computed from Dst index
  362.          */
  363.         LocalProvider(final double f10, final double f10B, final double s10,
  364.                       final double s10B, final double xm10, final double xm10B,
  365.                       final double y10, final double y10B, final double dstdtc) {
  366.             this.f10    = f10;
  367.             this.f10B   = f10B;
  368.             this.s10    = s10;
  369.             this.s10B   = s10B;
  370.             this.xm10   = xm10;
  371.             this.xm10B  = xm10B;
  372.             this.y10    = y10;
  373.             this.y10B   = y10B;
  374.             this.dstdtc = dstdtc;
  375.         }

  377.         /** {@inheritDoc} */
  378.         @Override
  379.         public AbsoluteDate getMinDate() {
  380.             return AbsoluteDate.PAST_INFINITY;
  381.         }

  383.         /** {@inheritDoc} */
  384.         @Override
  385.         public AbsoluteDate getMaxDate() {
  386.             return AbsoluteDate.FUTURE_INFINITY;
  387.         }

  389.         /** {@inheritDoc} */
  390.         @Override
  391.         public double getF10(final AbsoluteDate date) {
  392.             return f10;
  393.         }

  395.         /** {@inheritDoc} */
  396.         @Override
  397.         public double getF10B(final AbsoluteDate date) {
  398.             return f10B;
  399.         }

  401.         /** {@inheritDoc} */
  402.         @Override
  403.         public double getS10(final AbsoluteDate date) {
  404.             return s10;
  405.         }

  407.         /** {@inheritDoc} */
  408.         @Override
  409.         public double getS10B(final AbsoluteDate date) {
  410.             return s10B;
  411.         }

  413.         /** {@inheritDoc} */
  414.         @Override
  415.         public double getXM10(final AbsoluteDate date) {
  416.             return xm10;
  417.         }

  419.         /** {@inheritDoc} */
  420.         @Override
  421.         public double getXM10B(final AbsoluteDate date) {
  422.             return xm10B;
  423.         }

  425.         /** {@inheritDoc} */
  426.         @Override
  427.         public double getY10(final AbsoluteDate date) {
  428.             return y10;
  429.         }

  431.         /** {@inheritDoc} */
  432.         @Override
  433.         public double getY10B(final AbsoluteDate date) {
  434.             return y10B;
  435.         }

  437.         /** {@inheritDoc} */
  438.         @Override
  439.         public double getDSTDTC(final AbsoluteDate date) {
  440.             return dstdtc;
  441.         }
  442.     }
  443. }
Created with JaCoCo 0.8.12.202403310830