HarrisPriester.java
/* Copyright 2002-2023 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.models.earth.atmosphere;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
import org.hipparchus.util.SinCos;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.PVCoordinatesProvider;
/** This atmosphere model is the realization of the Modified Harris-Priester model.
* <p>
* This model is a static one that takes into account the diurnal density bulge.
* It doesn't need any space weather data but a density vs. altitude table, which
* depends on solar activity.
* </p>
* <p>
* The implementation relies on the book:<br>
* <b>Satellite Orbits</b><br>
* <i>Oliver Montenbruck, Eberhard Gill</i><br>
* Springer 2005
* </p>
* @author Pascal Parraud
*/
public class HarrisPriester implements Atmosphere {
/** Serializable UID.*/
private static final long serialVersionUID = 2772347498196369601L;
// Constants :
/** Default cosine exponent value. */
private static final int N_DEFAULT = 4;
/** Minimal value for calculating poxer of cosine. */
private static final double MIN_COS = 1.e-12;
/** Lag angle for diurnal bulge. */
private static final double LAG = FastMath.toRadians(30.0);
/** Lag angle sine and cosine. */
private static final SinCos SCLAG = FastMath.sinCos(LAG);
// CHECKSTYLE: stop NoWhitespaceAfter check
/** Harris-Priester min-max density (kg/m3) vs. altitude (m) table.
* These data are valid for a mean solar activity. */
private static final double[][] ALT_RHO = {
{ 100000.0, 4.974e-07, 4.974e-07 },
{ 120000.0, 2.490e-08, 2.490e-08 },
{ 130000.0, 8.377e-09, 8.710e-09 },
{ 140000.0, 3.899e-09, 4.059e-09 },
{ 150000.0, 2.122e-09, 2.215e-09 },
{ 160000.0, 1.263e-09, 1.344e-09 },
{ 170000.0, 8.008e-10, 8.758e-10 },
{ 180000.0, 5.283e-10, 6.010e-10 },
{ 190000.0, 3.617e-10, 4.297e-10 },
{ 200000.0, 2.557e-10, 3.162e-10 },
{ 210000.0, 1.839e-10, 2.396e-10 },
{ 220000.0, 1.341e-10, 1.853e-10 },
{ 230000.0, 9.949e-11, 1.455e-10 },
{ 240000.0, 7.488e-11, 1.157e-10 },
{ 250000.0, 5.709e-11, 9.308e-11 },
{ 260000.0, 4.403e-11, 7.555e-11 },
{ 270000.0, 3.430e-11, 6.182e-11 },
{ 280000.0, 2.697e-11, 5.095e-11 },
{ 290000.0, 2.139e-11, 4.226e-11 },
{ 300000.0, 1.708e-11, 3.526e-11 },
{ 320000.0, 1.099e-11, 2.511e-11 },
{ 340000.0, 7.214e-12, 1.819e-11 },
{ 360000.0, 4.824e-12, 1.337e-11 },
{ 380000.0, 3.274e-12, 9.955e-12 },
{ 400000.0, 2.249e-12, 7.492e-12 },
{ 420000.0, 1.558e-12, 5.684e-12 },
{ 440000.0, 1.091e-12, 4.355e-12 },
{ 460000.0, 7.701e-13, 3.362e-12 },
{ 480000.0, 5.474e-13, 2.612e-12 },
{ 500000.0, 3.916e-13, 2.042e-12 },
{ 520000.0, 2.819e-13, 1.605e-12 },
{ 540000.0, 2.042e-13, 1.267e-12 },
{ 560000.0, 1.488e-13, 1.005e-12 },
{ 580000.0, 1.092e-13, 7.997e-13 },
{ 600000.0, 8.070e-14, 6.390e-13 },
{ 620000.0, 6.012e-14, 5.123e-13 },
{ 640000.0, 4.519e-14, 4.121e-13 },
{ 660000.0, 3.430e-14, 3.325e-13 },
{ 680000.0, 2.632e-14, 2.691e-13 },
{ 700000.0, 2.043e-14, 2.185e-13 },
{ 720000.0, 1.607e-14, 1.779e-13 },
{ 740000.0, 1.281e-14, 1.452e-13 },
{ 760000.0, 1.036e-14, 1.190e-13 },
{ 780000.0, 8.496e-15, 9.776e-14 },
{ 800000.0, 7.069e-15, 8.059e-14 },
{ 840000.0, 4.680e-15, 5.741e-14 },
{ 880000.0, 3.200e-15, 4.210e-14 },
{ 920000.0, 2.210e-15, 3.130e-14 },
{ 960000.0, 1.560e-15, 2.360e-14 },
{ 1000000.0, 1.150e-15, 1.810e-14 }
};
// CHECKSTYLE: resume NoWhitespaceAfter check
/** Cosine exponent from 2 to 6 according to inclination. */
private double n;
/** Sun position. */
private PVCoordinatesProvider sun;
/** Earth body shape. */
private OneAxisEllipsoid earth;
/** Density table. */
private double[][] tabAltRho;
/** Simple constructor for Modified Harris-Priester atmosphere model.
* <p>The cosine exponent value is set to 4 by default.</p>
* <p>The default embedded density table is the one given in the referenced
* book from Montenbruck & Gill. It is given for mean solar activity and
* spreads over 100 to 1000 km.</p>
* @param sun the sun position
* @param earth the earth body shape
*/
public HarrisPriester(final PVCoordinatesProvider sun,
final OneAxisEllipsoid earth) {
this(sun, earth, ALT_RHO, N_DEFAULT);
}
/** Constructor for Modified Harris-Priester atmosphere model.
* <p>Recommanded values for the cosine exponent spread over the range
* 2, for low inclination orbits, to 6, for polar orbits.</p>
* <p> The default embedded density table is the one given in the referenced
* book from Montenbruck & Gill. It is given for mean solar activity and
* spreads over 100 to 1000 km. </p>
* @param sun the sun position
* @param earth the earth body shape
* @param n the cosine exponent
*/
public HarrisPriester(final PVCoordinatesProvider sun,
final OneAxisEllipsoid earth,
final double n) {
this(sun, earth, ALT_RHO, n);
}
/** Constructor for Modified Harris-Priester atmosphere model.
* <p>The provided density table must be an array such as:
* <ul>
* <li>tabAltRho[][0] = altitude (m)</li>
* <li>tabAltRho[][1] = min density (kg/m³)</li>
* <li>tabAltRho[][2] = max density (kg/m³)</li>
* </ul>
* <p> The altitude must be increasing without limitation in range. The
* internal density table is a copy of the provided one.
*
* <p>The cosine exponent value is set to 4 by default.</p>
* @param sun the sun position
* @param earth the earth body shape
* @param tabAltRho the density table
*/
public HarrisPriester(final PVCoordinatesProvider sun,
final OneAxisEllipsoid earth,
final double[][] tabAltRho) {
this(sun, earth, tabAltRho, N_DEFAULT);
}
/** Constructor for Modified Harris-Priester atmosphere model.
* <p>Recommanded values for the cosine exponent spread over the range
* 2, for low inclination orbits, to 6, for polar orbits.</p>
* <p>The provided density table must be an array such as:
* <ul>
* <li>tabAltRho[][0] = altitude (m)</li>
* <li>tabAltRho[][1] = min density (kg/m³)</li>
* <li>tabAltRho[][2] = max density (kg/m³)</li>
* </ul>
* <p> The altitude must be increasing without limitation in range. The
* internal density table is a copy of the provided one.
*
* @param sun the sun position
* @param earth the earth body shape
* @param tabAltRho the density table
* @param n the cosine exponent
*/
public HarrisPriester(final PVCoordinatesProvider sun,
final OneAxisEllipsoid earth,
final double[][] tabAltRho,
final double n) {
this.sun = sun;
this.earth = earth;
setTabDensity(tabAltRho);
setN(n);
}
/** {@inheritDoc} */
public Frame getFrame() {
return earth.getBodyFrame();
}
/** Set parameter N, the cosine exponent.
* @param n the cosine exponent
*/
private void setN(final double n) {
this.n = n;
}
/** Set a user define density table to deal with different solar activities.
* @param tab density vs. altitude table
*/
private void setTabDensity(final double[][] tab) {
this.tabAltRho = new double[tab.length][];
for (int i = 0; i < tab.length; i++) {
this.tabAltRho[i] = tab[i].clone();
}
}
/** Get the current density table.
* <p>The density table is an array such as:
* <ul>
* <li>tabAltRho[][0] = altitude (m)</li>
* <li>tabAltRho[][1] = min density (kg/m³)</li>
* <li>tabAltRho[][2] = max density (kg/m³)</li>
* </ul>
* <p> The altitude must be increasing without limitation in range.
*
* <p>
* The returned density table is a copy of the current one.
* </p>
* @return density vs. altitude table
*/
public double[][] getTabDensity() {
final double[][] copy = new double[tabAltRho.length][];
for (int i = 0; i < tabAltRho.length; i++) {
copy[i] = tabAltRho[i].clone();
}
return copy;
}
/** Get the minimal altitude for the model.
* <p>No computation is possible below this altitude.</p>
* @return the minimal altitude (m)
*/
public double getMinAlt() {
return tabAltRho[0][0];
}
/** Get the maximal altitude for the model.
* <p>Above this altitude, density is assumed to be zero.</p>
* @return the maximal altitude (m)
*/
public double getMaxAlt() {
return tabAltRho[tabAltRho.length - 1][0];
}
/** Get the local density.
* @param sunInEarth position of the Sun in Earth frame (m)
* @param posInEarth target position in Earth frame (m)
* @return the local density (kg/m³)
*/
public double getDensity(final Vector3D sunInEarth, final Vector3D posInEarth) {
final double posAlt = getHeight(posInEarth);
// Check for height boundaries
if (posAlt < getMinAlt()) {
throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
}
if (posAlt > getMaxAlt()) {
return 0.;
}
// Diurnal bulge apex direction
final Vector3D sunDir = sunInEarth.normalize();
final Vector3D bulDir = new Vector3D(sunDir.getX() * SCLAG.cos() - sunDir.getY() * SCLAG.sin(),
sunDir.getX() * SCLAG.sin() + sunDir.getY() * SCLAG.cos(),
sunDir.getZ());
// Cosine of angle Psi between the diurnal bulge apex and the satellite
final double cosPsi = bulDir.normalize().dotProduct(posInEarth.normalize());
// (1 + cos(Psi))/2 = cos²(Psi/2)
final double c2Psi2 = (1. + cosPsi) / 2.;
final double cPsi2 = FastMath.sqrt(c2Psi2);
final double cosPow = (cPsi2 > MIN_COS) ? c2Psi2 * FastMath.pow(cPsi2, n - 2) : 0.;
// Search altitude index in density table
int ia = 0;
while (ia < tabAltRho.length - 2 && posAlt > tabAltRho[ia + 1][0]) {
ia++;
}
// Fractional satellite height
final double dH = (tabAltRho[ia][0] - posAlt) / (tabAltRho[ia][0] - tabAltRho[ia + 1][0]);
// Min exponential density interpolation
final double rhoMin = tabAltRho[ia][1] * FastMath.pow(tabAltRho[ia + 1][1] / tabAltRho[ia][1], dH);
if (Precision.equals(cosPow, 0.)) {
return rhoMin;
} else {
// Max exponential density interpolation
final double rhoMax = tabAltRho[ia][2] * FastMath.pow(tabAltRho[ia + 1][2] / tabAltRho[ia][2], dH);
return rhoMin + (rhoMax - rhoMin) * cosPow;
}
}
/** Get the local density.
* @param sunInEarth position of the Sun in Earth frame (m)
* @param posInEarth target position in Earth frame (m)
* @return the local density (kg/m³)
* @param <T> instance of CalculusFieldElement<T>
*/
public <T extends CalculusFieldElement<T>> T getDensity(final Vector3D sunInEarth, final FieldVector3D<T> posInEarth) {
final T zero = posInEarth.getX().getField().getZero();
final T posAlt = getHeight(posInEarth);
// Check for height boundaries
if (posAlt.getReal() < getMinAlt()) {
throw new OrekitException(OrekitMessages.ALTITUDE_BELOW_ALLOWED_THRESHOLD, posAlt, getMinAlt());
}
if (posAlt.getReal() > getMaxAlt()) {
return zero;
}
// Diurnal bulge apex direction
final Vector3D sunDir = sunInEarth.normalize();
final Vector3D bulDir = new Vector3D(sunDir.getX() * SCLAG.cos() - sunDir.getY() * SCLAG.sin(),
sunDir.getX() * SCLAG.sin() + sunDir.getY() * SCLAG.cos(),
sunDir.getZ());
// Cosine of angle Psi between the diurnal bulge apex and the satellite
final T cosPsi = posInEarth.normalize().dotProduct(bulDir.normalize());
// (1 + cos(Psi))/2 = cos²(Psi/2)
final T c2Psi2 = cosPsi.add(1.).divide(2);
final T cPsi2 = c2Psi2.sqrt();
final T cosPow = (cPsi2.getReal() > MIN_COS) ? c2Psi2.multiply(cPsi2.pow(n - 2)) : zero;
// Search altitude index in density table
int ia = 0;
while (ia < tabAltRho.length - 2 && posAlt.getReal() > tabAltRho[ia + 1][0]) {
ia++;
}
// Fractional satellite height
final T dH = posAlt.negate().add(tabAltRho[ia][0]).divide(tabAltRho[ia][0] - tabAltRho[ia + 1][0]);
// Min exponential density interpolation
final T rhoMin = zero.add(tabAltRho[ia + 1][1] / tabAltRho[ia][1]).pow(dH).multiply(tabAltRho[ia][1]);
if (Precision.equals(cosPow.getReal(), 0.)) {
return zero.add(rhoMin);
} else {
// Max exponential density interpolation
final T rhoMax = zero.add(tabAltRho[ia + 1][2] / tabAltRho[ia][2]).pow(dH).multiply(tabAltRho[ia][2]);
return rhoMin.add(rhoMax.subtract(rhoMin).multiply(cosPow));
}
}
/** Get the local density at some position.
* @param date current date
* @param position current position
* @param frame the frame in which is defined the position
* @return local density (kg/m³)
* or if altitude is below the model minimal altitude
*/
public double getDensity(final AbsoluteDate date, final Vector3D position, final Frame frame) {
// Sun position in earth frame
final Vector3D sunInEarth = sun.getPosition(date, earth.getBodyFrame());
// Target position in earth frame
final Vector3D posInEarth = frame
.getStaticTransformTo(earth.getBodyFrame(), date)
.transformPosition(position);
return getDensity(sunInEarth, posInEarth);
}
/** Get the local density at some position.
* @param date current date
* @param position current position
* @param <T> implements a CalculusFieldElement
* @param frame the frame in which is defined the position
* @return local density (kg/m³)
* or if altitude is below the model minimal altitude
*/
public <T extends CalculusFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date,
final FieldVector3D<T> position,
final Frame frame) {
// Sun position in earth frame
final Vector3D sunInEarth = sun.getPosition(date.toAbsoluteDate(), earth.getBodyFrame());
// Target position in earth frame
final FieldVector3D<T> posInEarth = frame
.getStaticTransformTo(earth.getBodyFrame(), date.toAbsoluteDate())
.transformPosition(position);
return getDensity(sunInEarth, posInEarth);
}
/** Get the height above the Earth for the given position.
* <p>
* The height computation is an approximation valid for the considered atmosphere.
* </p>
* @param position current position in Earth frame
* @return height (m)
*/
private double getHeight(final Vector3D position) {
final double a = earth.getEquatorialRadius();
final double f = earth.getFlattening();
final double e2 = f * (2. - f);
final double r = position.getNorm();
final double sl = position.getZ() / r;
final double cl2 = 1. - sl * sl;
final double coef = FastMath.sqrt((1. - e2) / (1. - e2 * cl2));
return r - a * coef;
}
/** Get the height above the Earth for the given position.
* <p>
* The height computation is an approximation valid for the considered atmosphere.
* </p>
* @param position current position in Earth frame
* @param <T> instance of CalculusFieldElement<T>
* @return height (m)
*/
private <T extends CalculusFieldElement<T>> T getHeight(final FieldVector3D<T> position) {
final double a = earth.getEquatorialRadius();
final double f = earth.getFlattening();
final double e2 = f * (2. - f);
final T r = position.getNorm();
final T sl = position.getZ().divide(r);
final T cl2 = sl.multiply(sl).negate().add(1.);
final T coef = cl2.multiply(-e2).add(1.).reciprocal().multiply(1. - e2).sqrt();
return r.subtract(coef.multiply(a));
}
}