HarrisPriester.java

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


  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 implements Atmosphere {

  50.     /** Serializable UID.*/
  51.     private static final long serialVersionUID = 2772347498196369601L;

  53.     // Constants :

  55.     /** Default cosine exponent value. */
  56.     private static final int N_DEFAULT = 4;

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

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

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

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

  124.     /** Cosine exponent from 2 to 6 according to inclination. */
  125.     private double n;

  127.     /** Sun position. */
  128.     private PVCoordinatesProvider sun;

  130.     /** Earth body shape. */
  131.     private OneAxisEllipsoid earth;

  133.     /** Density table. */
  134.     private double[][] tabAltRho;

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

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

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

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

  213.     /** {@inheritDoc} */
  214.     public Frame getFrame() {
  215.         return earth.getBodyFrame();
  216.     }

  218.     /** Set parameter N, the cosine exponent.
  219.      *  @param n the cosine exponent
  220.      */
  221.     private void setN(final double n) {
  222.         this.n = n;
  223.     }

  225.     /** Set a user define density table to deal with different solar activities.
  226.      *  @param tab density vs. altitude table
  227.      */
  228.     private void setTabDensity(final double[][] tab) {
  229.         this.tabAltRho = new double[tab.length][];
  230.         for (int i = 0; i < tab.length; i++) {
  231.             this.tabAltRho[i] = tab[i].clone();
  232.         }
  233.     }

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

  257.     /** Get the minimal altitude for the model.
  258.      * <p>No computation is possible below this altitude.</p>
  259.      *  @return the minimal altitude (m)
  260.      */
  261.     public double getMinAlt() {
  262.         return tabAltRho[0][0];
  263.     }

  265.     /** Get the maximal altitude for the model.
  266.      * <p>Above this altitude, density is assumed to be zero.</p>
  267.      *  @return the maximal altitude (m)
  268.      */
  269.     public double getMaxAlt() {
  270.         return tabAltRho[tabAltRho.length - 1][0];
  271.     }

  273.     /** Get the local density.
  274.      * @param sunInEarth position of the Sun in Earth frame (m)
  275.      * @param posInEarth target position in Earth frame (m)
  276.      * @return the local density (kg/m³)
  277.      */
  278.     public double getDensity(final Vector3D sunInEarth, final Vector3D posInEarth) {

  280.         final double posAlt = getHeight(posInEarth);
  281.         // Check for height boundaries
  282.         if (posAlt < getMinAlt()) {
  283.             throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
  284.         }
  285.         if (posAlt > getMaxAlt()) {
  286.             return 0.;
  287.         }

  289.         // Diurnal bulge apex direction
  290.         final Vector3D sunDir = sunInEarth.normalize();
  291.         final Vector3D bulDir = new Vector3D(sunDir.getX() * SCLAG.cos() - sunDir.getY() * SCLAG.sin(),
  292.                                              sunDir.getX() * SCLAG.sin() + sunDir.getY() * SCLAG.cos(),
  293.                                              sunDir.getZ());

  295.         // Cosine of angle Psi between the diurnal bulge apex and the satellite
  296.         final double cosPsi = bulDir.normalize().dotProduct(posInEarth.normalize());
  297.         // (1 + cos(Psi))/2 = cos²(Psi/2)
  298.         final double c2Psi2 = (1. + cosPsi) / 2.;
  299.         final double cPsi2  = FastMath.sqrt(c2Psi2);
  300.         final double cosPow = (cPsi2 > MIN_COS) ? c2Psi2 * FastMath.pow(cPsi2, n - 2) : 0.;

  302.         // Search altitude index in density table
  303.         int ia = 0;
  304.         while (ia < tabAltRho.length - 2 && posAlt > tabAltRho[ia + 1][0]) {
  305.             ia++;
  306.         }

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

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

  314.         if (Precision.equals(cosPow, 0.)) {
  315.             return rhoMin;
  316.         } else {
  317.             // Max exponential density interpolation
  318.             final double rhoMax = tabAltRho[ia][2] * FastMath.pow(tabAltRho[ia + 1][2] / tabAltRho[ia][2], dH);
  319.             return rhoMin + (rhoMax - rhoMin) * cosPow;
  320.         }

  322.     }

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

  341.         // Diurnal bulge apex direction
  342.         final Vector3D sunDir = sunInEarth.normalize();
  343.         final Vector3D bulDir = new Vector3D(sunDir.getX() * SCLAG.cos() - sunDir.getY() * SCLAG.sin(),
  344.                                              sunDir.getX() * SCLAG.sin() + sunDir.getY() * SCLAG.cos(),
  345.                                              sunDir.getZ());

  347.         // Cosine of angle Psi between the diurnal bulge apex and the satellite
  348.         final T cosPsi = posInEarth.normalize().dotProduct(bulDir.normalize());
  349.         // (1 + cos(Psi))/2 = cos²(Psi/2)
  350.         final T c2Psi2 = cosPsi.add(1.).divide(2);
  351.         final T cPsi2  = c2Psi2.sqrt();
  352.         final T cosPow = (cPsi2.getReal() > MIN_COS) ? c2Psi2.multiply(cPsi2.pow(n - 2)) : zero;

  354.         // Search altitude index in density table
  355.         int ia = 0;
  356.         while (ia < tabAltRho.length - 2 && posAlt.getReal() > tabAltRho[ia + 1][0]) {
  357.             ia++;
  358.         }

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

  363.         // Min exponential density interpolation
  364.         final T rhoMin = zero.add(tabAltRho[ia + 1][1] / tabAltRho[ia][1]).pow(dH).multiply(tabAltRho[ia][1]);

  366.         if (Precision.equals(cosPow.getReal(), 0.)) {
  367.             return zero.add(rhoMin);
  368.         } else {
  369.             // Max exponential density interpolation
  370.             final T rhoMax = zero.add(tabAltRho[ia + 1][2] / tabAltRho[ia][2]).pow(dH).multiply(tabAltRho[ia][2]);
  371.             return rhoMin.add(rhoMax.subtract(rhoMin).multiply(cosPow));
  372.         }

  374.     }

  376.     /** Get the local density at some position.
  377.      * @param date current date
  378.      * @param position current position
  379.      * @param frame the frame in which is defined the position
  380.      * @return local density (kg/m³)
  381.           *            or if altitude is below the model minimal altitude
  382.      */
  383.     public double getDensity(final AbsoluteDate date, final Vector3D position, final Frame frame) {

  385.         // Sun position in earth frame
  386.         final Vector3D sunInEarth = sun.getPosition(date, earth.getBodyFrame());

  388.         // Target position in earth frame
  389.         final Vector3D posInEarth = frame
  390.                 .getStaticTransformTo(earth.getBodyFrame(), date)
  391.                 .transformPosition(position);

  393.         return getDensity(sunInEarth, posInEarth);
  394.     }

  396.     /** Get the local density at some position.
  397.      * @param date current date
  398.      * @param position current position
  399.      * @param <T> implements a CalculusFieldElement
  400.      * @param frame the frame in which is defined the position
  401.      * @return local density (kg/m³)
  402.           *            or if altitude is below the model minimal altitude
  403.      */
  404.     public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
  405.                                                         final FieldVector3D<T> position,
  406.                                                         final Frame frame) {
  407.         // Sun position in earth frame
  408.         final Vector3D sunInEarth = sun.getPosition(date.toAbsoluteDate(), earth.getBodyFrame());

  410.         // Target position in earth frame
  411.         final FieldVector3D<T> posInEarth = frame
  412.                 .getStaticTransformTo(earth.getBodyFrame(), date.toAbsoluteDate())
  413.                 .transformPosition(position);

  415.         return getDensity(sunInEarth, posInEarth);
  416.     }

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

  434.         return r - a * coef;
  435.     }

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

  454.         return r.subtract(coef.multiply(a));
  455.     }

  457. }
Created with JaCoCo 0.8.11.202310140853