HarrisPriester.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.geometry.euclidean.threed.FieldVector3D;
  21. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.Precision;
  24. import org.hipparchus.util.SinCos;
  25. import org.orekit.bodies.OneAxisEllipsoid;
  26. import org.orekit.errors.OrekitException;
  27. import org.orekit.errors.OrekitMessages;
  28. import org.orekit.frames.Frame;
  29. import org.orekit.time.AbsoluteDate;
  30. import org.orekit.time.FieldAbsoluteDate;
  31. import org.orekit.utils.ExtendedPositionProvider;


  34. /** This atmosphere model is the realization of the Modified Harris-Priester model.
  35.  * <p>
  36.  * This model is a static one that takes into account the diurnal density bulge.
  37.  * It doesn't need any space weather data but a density vs. altitude table, which
  38.  * depends on solar activity.
  39.  * </p>
  40.  * <p>
  41.  * The implementation relies on the book:<br>
  42.  * <b>Satellite Orbits</b><br>
  43.  * <i>Oliver Montenbruck, Eberhard Gill</i><br>
  44.  * Springer 2005
  45.  * </p>
  46.  * @author Pascal Parraud
  47.  */
  48. public class HarrisPriester extends AbstractSunInfluencedAtmosphere {

  50.     // Constants :

  52.     /** Default cosine exponent value. */
  53.     private static final int N_DEFAULT = 4;

  55.     /** Minimal value for calculating poxer of cosine. */
  56.     private static final double MIN_COS = 1.e-12;

  58.     /** Lag angle for diurnal bulge. */
  59.     private static final double LAG = FastMath.toRadians(30.0);

  61.     /** Lag angle sine and cosine. */
  62.     private static final SinCos SCLAG = FastMath.sinCos(LAG);

  64.     // CHECKSTYLE: stop NoWhitespaceAfter check
  65.     /** Harris-Priester min-max density (kg/m3) vs. altitude (m) table.
  66.      *  These data are valid for a mean solar activity. */
  67.     private static final double[][] ALT_RHO = {
  68.         {  100000.0, 4.974e-07, 4.974e-07 },
  69.         {  120000.0, 2.490e-08, 2.490e-08 },
  70.         {  130000.0, 8.377e-09, 8.710e-09 },
  71.         {  140000.0, 3.899e-09, 4.059e-09 },
  72.         {  150000.0, 2.122e-09, 2.215e-09 },
  73.         {  160000.0, 1.263e-09, 1.344e-09 },
  74.         {  170000.0, 8.008e-10, 8.758e-10 },
  75.         {  180000.0, 5.283e-10, 6.010e-10 },
  76.         {  190000.0, 3.617e-10, 4.297e-10 },
  77.         {  200000.0, 2.557e-10, 3.162e-10 },
  78.         {  210000.0, 1.839e-10, 2.396e-10 },
  79.         {  220000.0, 1.341e-10, 1.853e-10 },
  80.         {  230000.0, 9.949e-11, 1.455e-10 },
  81.         {  240000.0, 7.488e-11, 1.157e-10 },
  82.         {  250000.0, 5.709e-11, 9.308e-11 },
  83.         {  260000.0, 4.403e-11, 7.555e-11 },
  84.         {  270000.0, 3.430e-11, 6.182e-11 },
  85.         {  280000.0, 2.697e-11, 5.095e-11 },
  86.         {  290000.0, 2.139e-11, 4.226e-11 },
  87.         {  300000.0, 1.708e-11, 3.526e-11 },
  88.         {  320000.0, 1.099e-11, 2.511e-11 },
  89.         {  340000.0, 7.214e-12, 1.819e-11 },
  90.         {  360000.0, 4.824e-12, 1.337e-11 },
  91.         {  380000.0, 3.274e-12, 9.955e-12 },
  92.         {  400000.0, 2.249e-12, 7.492e-12 },
  93.         {  420000.0, 1.558e-12, 5.684e-12 },
  94.         {  440000.0, 1.091e-12, 4.355e-12 },
  95.         {  460000.0, 7.701e-13, 3.362e-12 },
  96.         {  480000.0, 5.474e-13, 2.612e-12 },
  97.         {  500000.0, 3.916e-13, 2.042e-12 },
  98.         {  520000.0, 2.819e-13, 1.605e-12 },
  99.         {  540000.0, 2.042e-13, 1.267e-12 },
  100.         {  560000.0, 1.488e-13, 1.005e-12 },
  101.         {  580000.0, 1.092e-13, 7.997e-13 },
  102.         {  600000.0, 8.070e-14, 6.390e-13 },
  103.         {  620000.0, 6.012e-14, 5.123e-13 },
  104.         {  640000.0, 4.519e-14, 4.121e-13 },
  105.         {  660000.0, 3.430e-14, 3.325e-13 },
  106.         {  680000.0, 2.632e-14, 2.691e-13 },
  107.         {  700000.0, 2.043e-14, 2.185e-13 },
  108.         {  720000.0, 1.607e-14, 1.779e-13 },
  109.         {  740000.0, 1.281e-14, 1.452e-13 },
  110.         {  760000.0, 1.036e-14, 1.190e-13 },
  111.         {  780000.0, 8.496e-15, 9.776e-14 },
  112.         {  800000.0, 7.069e-15, 8.059e-14 },
  113.         {  840000.0, 4.680e-15, 5.741e-14 },
  114.         {  880000.0, 3.200e-15, 4.210e-14 },
  115.         {  920000.0, 2.210e-15, 3.130e-14 },
  116.         {  960000.0, 1.560e-15, 2.360e-14 },
  117.         { 1000000.0, 1.150e-15, 1.810e-14 }
  118.     };
  119.     // CHECKSTYLE: resume NoWhitespaceAfter check

  121.     /** Cosine exponent from 2 to 6 according to inclination. */
  122.     private double n;

  124.     /** Earth body shape. */
  125.     private OneAxisEllipsoid earth;

  127.     /** Density table. */
  128.     private double[][] tabAltRho;

  130.     /** Simple constructor for Modified Harris-Priester atmosphere model.
  131.      *  <p>The cosine exponent value is set to 4 by default.</p>
  132.      *  <p>The default embedded density table is the one given in the referenced
  133.      *  book from Montenbruck &amp; Gill. It is given for mean solar activity and
  134.      *  spreads over 100 to 1000 km.</p>
  135.      * @param sun the sun position
  136.      * @param earth the earth body shape
  137.      */
  138.     public HarrisPriester(final ExtendedPositionProvider sun,
  139.                           final OneAxisEllipsoid earth) {
  140.         this(sun, earth, ALT_RHO, N_DEFAULT);
  141.     }

  143.     /** Constructor for Modified Harris-Priester atmosphere model.
  144.      *  <p>Recommanded values for the cosine exponent spread over the range
  145.      *  2, for low inclination orbits, to 6, for polar orbits.</p>
  146.      *  <p> The default embedded density table is the one given in the referenced
  147.      *  book from Montenbruck &amp; Gill. It is given for mean solar activity and
  148.      *  spreads over 100 to 1000 km. </p>
  149.      *  @param sun the sun position
  150.      * @param earth the earth body shape
  151.      * @param n the cosine exponent
  152.      */
  153.     public HarrisPriester(final ExtendedPositionProvider sun,
  154.                           final OneAxisEllipsoid earth,
  155.                           final double n) {
  156.         this(sun, earth, ALT_RHO, n);
  157.     }

  159.     /** Constructor for Modified Harris-Priester atmosphere model.
  160.      *  <p>The provided density table must be an array such as:
  161.      *  <ul>
  162.      *   <li>tabAltRho[][0] = altitude (m)</li>
  163.      *   <li>tabAltRho[][1] = min density (kg/m³)</li>
  164.      *   <li>tabAltRho[][2] = max density (kg/m³)</li>
  165.      *  </ul>
  166.      *  <p> The altitude must be increasing without limitation in range. The
  167.      *  internal density table is a copy of the provided one.
  168.      *
  169.      *  <p>The cosine exponent value is set to 4 by default.</p>
  170.      * @param sun the sun position
  171.      * @param earth the earth body shape
  172.      * @param tabAltRho the density table
  173.      */
  174.     public HarrisPriester(final ExtendedPositionProvider sun,
  175.                           final OneAxisEllipsoid earth,
  176.                           final double[][] tabAltRho) {
  177.         this(sun, earth, tabAltRho, N_DEFAULT);
  178.     }

  180.     /** Constructor for Modified Harris-Priester atmosphere model.
  181.      *  <p>Recommanded values for the cosine exponent spread over the range
  182.      *  2, for low inclination orbits, to 6, for polar orbits.</p>
  183.      *  <p>The provided density table must be an array such as:
  184.      *  <ul>
  185.      *   <li>tabAltRho[][0] = altitude (m)</li>
  186.      *   <li>tabAltRho[][1] = min density (kg/m³)</li>
  187.      *   <li>tabAltRho[][2] = max density (kg/m³)</li>
  188.      *  </ul>
  189.      *  <p> The altitude must be increasing without limitation in range. The
  190.      *  internal density table is a copy of the provided one.
  191.      *
  192.      *  @param sun the sun position
  193.      * @param earth the earth body shape
  194.      * @param tabAltRho the density table
  195.      * @param n the cosine exponent
  196.      */
  197.     public HarrisPriester(final ExtendedPositionProvider sun,
  198.                           final OneAxisEllipsoid earth,
  199.                           final double[][] tabAltRho,
  200.                           final double n) {
  201.         super(sun);
  202.         this.earth = earth;
  203.         setTabDensity(tabAltRho);
  204.         setN(n);
  205.     }

  207.     /** {@inheritDoc} */
  208.     public Frame getFrame() {
  209.         return earth.getBodyFrame();
  210.     }

  212.     /** Set parameter N, the cosine exponent.
  213.      *  @param n the cosine exponent
  214.      */
  215.     private void setN(final double n) {
  216.         this.n = n;
  217.     }

  219.     /** Set a user define density table to deal with different solar activities.
  220.      *  @param tab density vs. altitude table
  221.      */
  222.     private void setTabDensity(final double[][] tab) {
  223.         this.tabAltRho = new double[tab.length][];
  224.         for (int i = 0; i < tab.length; i++) {
  225.             this.tabAltRho[i] = tab[i].clone();
  226.         }
  227.     }

  229.     /** Get the current density table.
  230.      *  <p>The density table is an array such as:
  231.      *  <ul>
  232.      *   <li>tabAltRho[][0] = altitude (m)</li>
  233.      *   <li>tabAltRho[][1] = min density (kg/m³)</li>
  234.      *   <li>tabAltRho[][2] = max density (kg/m³)</li>
  235.      *  </ul>
  236.      *  <p> The altitude must be increasing without limitation in range.
  237.      *
  238.      *  <p>
  239.      *  The returned density table is a copy of the current one.
  240.      *  </p>
  241.      *  @return density vs. altitude table
  242.      */
  243.     public double[][] getTabDensity() {
  244.         final double[][] copy = new double[tabAltRho.length][];
  245.         for (int i = 0; i < tabAltRho.length; i++) {
  246.             copy[i] = tabAltRho[i].clone();
  247.         }
  248.         return copy;
  249.     }

  251.     /** Get the minimal altitude for the model.
  252.      * <p>No computation is possible below this altitude.</p>
  253.      *  @return the minimal altitude (m)
  254.      */
  255.     public double getMinAlt() {
  256.         return tabAltRho[0][0];
  257.     }

  259.     /** Get the maximal altitude for the model.
  260.      * <p>Above this altitude, density is assumed to be zero.</p>
  261.      *  @return the maximal altitude (m)
  262.      */
  263.     public double getMaxAlt() {
  264.         return tabAltRho[tabAltRho.length - 1][0];
  265.     }

  267.     /** Get the local density.
  268.      * @param sunInEarth position of the Sun in Earth frame (m)
  269.      * @param posInEarth target position in Earth frame (m)
  270.      * @return the local density (kg/m³)
  271.      */
  272.     public double getDensity(final Vector3D sunInEarth, final Vector3D posInEarth) {

  274.         final double posAlt = getHeight(posInEarth);
  275.         // Check for height boundaries
  276.         if (posAlt < getMinAlt()) {
  277.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
  278.         }
  279.         if (posAlt > getMaxAlt()) {
  280.             return 0.;
  281.         }

  283.         // Diurnal bulge apex direction
  284.         final Vector3D sunDir = sunInEarth.normalize();
  285.         final Vector3D bulDir = new Vector3D(sunDir.getX() * SCLAG.cos() - sunDir.getY() * SCLAG.sin(),
  286.                                              sunDir.getX() * SCLAG.sin() + sunDir.getY() * SCLAG.cos(),
  287.                                              sunDir.getZ());

  289.         // Cosine of angle Psi between the diurnal bulge apex and the satellite
  290.         final double cosPsi = bulDir.normalize().dotProduct(posInEarth.normalize());
  291.         // (1 + cos(Psi))/2 = cos²(Psi/2)
  292.         final double c2Psi2 = (1. + cosPsi) / 2.;
  293.         final double cPsi2  = FastMath.sqrt(c2Psi2);
  294.         final double cosPow = (cPsi2 > MIN_COS) ? c2Psi2 * FastMath.pow(cPsi2, n - 2) : 0.;

  296.         // Search altitude index in density table
  297.         int ia = 0;
  298.         while (ia < tabAltRho.length - 2 && posAlt > tabAltRho[ia + 1][0]) {
  299.             ia++;
  300.         }

  302.         // Fractional satellite height
  303.         final double dH = (tabAltRho[ia][0] - posAlt) / (tabAltRho[ia][0] - tabAltRho[ia + 1][0]);

  305.         // Min exponential density interpolation
  306.         final double rhoMin = tabAltRho[ia][1] * FastMath.pow(tabAltRho[ia + 1][1] / tabAltRho[ia][1], dH);

  308.         if (Precision.equals(cosPow, 0.)) {
  309.             return rhoMin;
  310.         } else {
  311.             // Max exponential density interpolation
  312.             final double rhoMax = tabAltRho[ia][2] * FastMath.pow(tabAltRho[ia + 1][2] / tabAltRho[ia][2], dH);
  313.             return rhoMin + (rhoMax - rhoMin) * cosPow;
  314.         }

  316.     }

  318.     /** Get the local density.
  319.      * @param sunInEarth position of the Sun in Earth frame (m)
  320.      * @param posInEarth target position in Earth frame (m)
  321.      * @return the local density (kg/m³)
  322.      * @param <T> instance of CalculusFieldElement&lt;T&gt;
  323.      */
  324.     public <T extends CalculusFieldElement<T>> T getDensity(final FieldVector3D<T> sunInEarth, final FieldVector3D<T> posInEarth) {
  325.         final T zero = posInEarth.getX().getField().getZero();
  326.         final T posAlt = getHeight(posInEarth);
  327.         // Check for height boundaries
  328.         if (posAlt.getReal() < getMinAlt()) {
  329.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
  330.         }
  331.         if (posAlt.getReal() > getMaxAlt()) {
  332.             return zero;
  333.         }

  335.         // Diurnal bulge apex direction
  336.         final FieldVector3D<T> sunDir = sunInEarth.normalize();
  337.         final FieldVector3D<T> bulDir = new FieldVector3D<>(sunDir.getX().multiply(SCLAG.cos()).subtract(sunDir.getY().multiply(SCLAG.sin())),
  338.                                              sunDir.getX().multiply(SCLAG.sin()).add(sunDir.getY().multiply(SCLAG.cos())),
  339.                                              sunDir.getZ());

  341.         // Cosine of angle Psi between the diurnal bulge apex and the satellite
  342.         final T cosPsi = posInEarth.normalize().dotProduct(bulDir.normalize());
  343.         // (1 + cos(Psi))/2 = cos²(Psi/2)
  344.         final T c2Psi2 = cosPsi.add(1.).divide(2);
  345.         final T cPsi2  = c2Psi2.sqrt();
  346.         final T cosPow = (cPsi2.getReal() > MIN_COS) ? c2Psi2.multiply(cPsi2.pow(n - 2)) : zero;

  348.         // Search altitude index in density table
  349.         int ia = 0;
  350.         while (ia < tabAltRho.length - 2 && posAlt.getReal() > tabAltRho[ia + 1][0]) {
  351.             ia++;
  352.         }

  354.         // Fractional satellite height
  355.         final T dH = posAlt.negate().add(tabAltRho[ia][0]).divide(tabAltRho[ia][0] - tabAltRho[ia + 1][0]);

  357.         // Min exponential density interpolation
  358.         final T rhoMin = zero.newInstance(tabAltRho[ia + 1][1] / tabAltRho[ia][1]).pow(dH).multiply(tabAltRho[ia][1]);

  360.         if (Precision.equals(cosPow.getReal(), 0.)) {
  361.             return rhoMin;
  362.         } else {
  363.             // Max exponential density interpolation
  364.             final T rhoMax = zero.newInstance(tabAltRho[ia + 1][2] / tabAltRho[ia][2]).pow(dH).multiply(tabAltRho[ia][2]);
  365.             return rhoMin.add(rhoMax.subtract(rhoMin).multiply(cosPow));
  366.         }

  368.     }

  370.     /** Get the local density at some position.
  371.      * @param date current date
  372.      * @param position current position
  373.      * @param frame the frame in which is defined the position
  374.      * @return local density (kg/m³)
  375.           *            or if altitude is below the model minimal altitude
  376.      */
  377.     public double getDensity(final AbsoluteDate date, final Vector3D position, final Frame frame) {

  379.         // Sun position in earth frame
  380.         final Vector3D sunInEarth = getSunPosition(date, earth.getBodyFrame());

  382.         // Target position in earth frame
  383.         final Vector3D posInEarth = frame
  384.                 .getStaticTransformTo(earth.getBodyFrame(), date)
  385.                 .transformPosition(position);

  387.         return getDensity(sunInEarth, posInEarth);
  388.     }

  390.     /** Get the local density at some position.
  391.      * @param date current date
  392.      * @param position current position
  393.      * @param <T> implements a CalculusFieldElement
  394.      * @param frame the frame in which is defined the position
  395.      * @return local density (kg/m³)
  396.           *            or if altitude is below the model minimal altitude
  397.      */
  398.     public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
  399.                                                             final FieldVector3D<T> position, final Frame frame) {
  400.         // Sun position in earth frame
  401.         final FieldVector3D<T> sunInEarth = getSunPosition(date, earth.getBodyFrame());

  403.         // Target position in earth frame
  404.         final FieldVector3D<T> posInEarth = frame
  405.                 .getStaticTransformTo(earth.getBodyFrame(), date.toAbsoluteDate())
  406.                 .transformPosition(position);

  408.         return getDensity(sunInEarth, posInEarth);
  409.     }

  411.     /** Get the height above the Earth for the given position.
  412.      *  <p>
  413.      *  The height computation is an approximation valid for the considered atmosphere.
  414.      *  </p>
  415.      *  @param position current position in Earth frame
  416.      *  @return height (m)
  417.      */
  418.     private double getHeight(final Vector3D position) {
  419.         final double a    = earth.getEquatorialRadius();
  420.         final double f    = earth.getFlattening();
  421.         final double e2   = f * (2. - f);
  422.         final double r    = position.getNorm();
  423.         final double sl   = position.getZ() / r;
  424.         final double cl2  = 1. - sl * sl;
  425.         final double coef = FastMath.sqrt((1. - e2) / (1. - e2 * cl2));

  427.         return r - a * coef;
  428.     }

  430.     /** Get the height above the Earth for the given position.
  431.      *  <p>
  432.      *  The height computation is an approximation valid for the considered atmosphere.
  433.      *  </p>
  434.      *  @param position current position in Earth frame
  435.      *  @param <T> instance of CalculusFieldElement<T>
  436.      *  @return height (m)
  437.      */
  438.     private <T extends CalculusFieldElement<T>> T getHeight(final FieldVector3D<T> position) {
  439.         final double a    = earth.getEquatorialRadius();
  440.         final double f    = earth.getFlattening();
  441.         final double e2   = f * (2. - f);
  442.         final T r    = position.getNorm();
  443.         final T sl   = position.getZ().divide(r);
  444.         final T cl2  = sl.square().negate().add(1.);
  445.         final T coef = cl2.multiply(-e2).add(1.).reciprocal().multiply(1. - e2).sqrt();

  447.         return r.subtract(coef.multiply(a));
  448.     }

  450. }
Created with JaCoCo 0.8.12.202403310830