HarrisPriester.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.apache.commons.math3.util.Precision;
  22. import org.orekit.bodies.OneAxisEllipsoid;
  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.PVCoordinates;
  29. import org.orekit.utils.PVCoordinatesProvider;


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

  48.     /** Serializable UID.*/
  49.     private static final long serialVersionUID = 2772347498196369601L;

  51.     // Constants :

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

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

  59.     /** Lag angle for diurnal bulge. */
  60.     private static final double LAG = FastMath.toRadians(30.0);
  61.     /** Lag angle cosine. */
  62.     private static final double COSLAG = FastMath.cos(LAG);
  63.     /** Lag angle sine. */
  64.     private static final double SINLAG = FastMath.sin(LAG);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  272.     /** Get the local density.
  273.      * @param sunInEarth position of the Sun in Earth frame (m)
  274.      * @param posInEarth target position in Earth frame (m)
  275.      * @return the local density (kg/m<sup>3</sup>)
  276.      * @exception OrekitException if altitude is below the model minimal altitude
  277.      */
  278.     public double getDensity(final Vector3D sunInEarth, final Vector3D posInEarth)
  279.         throws OrekitException {

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

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

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

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

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

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

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

  323.     }

  325.     /** Get the local density at some position.
  326.      * @param date current date
  327.      * @param position current position
  328.      * @param frame the frame in which is defined the position
  329.      * @return local density (kg/m<sup>3</sup>)
  330.      * @exception OrekitException if some frame conversion cannot be performed
  331.      *            or if altitude is below the model minimal altitude
  332.      */
  333.     public double getDensity(final AbsoluteDate date, final Vector3D position, final Frame frame)
  334.         throws OrekitException {

  336.         // Sun position in earth frame
  337.         final Vector3D sunInEarth = sun.getPVCoordinates(date, earth.getBodyFrame()).getPosition();

  339.         // Target position in earth frame
  340.         final Vector3D posInEarth = frame.getTransformTo(earth.getBodyFrame(), date).transformPosition(position);

  342.         return getDensity(sunInEarth, posInEarth);
  343.     }

  345.     /** Get the inertial velocity of atmosphere molecules.
  346.      * <p>
  347.      * Here the case is simplified : atmosphere is supposed to have a null velocity
  348.      * in earth frame.
  349.      * </p>
  350.      * @param date current date
  351.      * @param position current position
  352.      * @param frame the frame in which is defined the position
  353.      * @return velocity (m/s) (defined in the same frame as the position)
  354.      * @exception OrekitException if some frame conversion cannot be performed
  355.      */
  356.     public Vector3D getVelocity(final AbsoluteDate date, final Vector3D position, final Frame frame)
  357.         throws OrekitException {
  358.         final Transform bodyToFrame = earth.getBodyFrame().getTransformTo(frame, date);
  359.         final Vector3D posInBody    = bodyToFrame.getInverse().transformPosition(position);
  360.         final PVCoordinates pvBody  = new PVCoordinates(posInBody, Vector3D.ZERO);
  361.         final PVCoordinates pvFrame = bodyToFrame.transformPVCoordinates(pvBody);
  362.         return pvFrame.getVelocity();
  363.     }

  365.     /** Get the height above the Earth for the given position.
  366.      *  <p>
  367.      *  The height computation is an approximation valid for the considered atmosphere.
  368.      *  </p>
  369.      *  @param position current position in Earth frame
  370.      *  @return height (m)
  371.      */
  372.     private double getHeight(final Vector3D position) {
  373.         final double a    = earth.getEquatorialRadius();
  374.         final double f    = earth.getFlattening();
  375.         final double e2   = f * (2. - f);
  376.         final double r    = position.getNorm();
  377.         final double sl   = position.getZ() / r;
  378.         final double cl2  = 1. - sl * sl;
  379.         final double coef = FastMath.sqrt((1. - e2) / (1. - e2 * cl2));

  381.         return r - a * coef;
  382.     }

  384. }
Created with JaCoCo 0.6.3.201306030806