ReferenceEllipsoid.java

  1. /* Contributed in the public domain.
  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;

  19. import org.hipparchus.util.FastMath;
  20. import org.orekit.bodies.OneAxisEllipsoid;
  21. import org.orekit.frames.Frame;
  22. import org.orekit.utils.Constants;

  24. /**
  25.  * A Reference Ellipsoid for use in geodesy. The ellipsoid defines an
  26.  * ellipsoidal potential called the normal potential, and its gradient, normal
  27.  * gravity.
  28.  *
  29.  * <p> These parameters are needed to define the normal potential:
  30.  *
  31.  *
  32.  * <ul> <li>a, semi-major axis</li>
  33.  *
  34.  * <li>f, flattening</li>
  35.  *
  36.  * <li>GM, the gravitational parameter</li>
  37.  *
  38.  * <li>&omega;, the spin rate</li> </ul>
  39.  *
  40.  * <p> References:
  41.  *
  42.  * <ol> <li>Martin Losch, Verena Seufer. How to Compute Geoid Undulations (Geoid
  43.  * Height Relative to a Given Reference Ellipsoid) from Spherical Harmonic
  44.  * Coefficients for Satellite Altimetry Applications. , 2003. <a
  45.  * href="http://mitgcm.org/~mlosch/geoidcookbook.pdf" >mitgcm.org/~mlosch/geoidcookbook.pdf</a></li>
  46.  *
  47.  * <li>Weikko A. Heiskanen, Helmut Moritz. Physical Geodesy. W. H. Freeman and
  48.  * Company, 1967. (especially sections 2.13 and equation 2-144)</li>
  49.  *
  50.  * <li>Department of Defense World Geodetic System 1984. 2000. NIMA TR 8350.2
  51.  * Third Edition, Amendment 1.</li> </ol>
  52.  *
  53.  * @author Evan Ward
  54.  * @author Guylaine Prat
  55.  */
  56. public class ReferenceEllipsoid extends OneAxisEllipsoid implements EarthShape {

  58.     /** uid is date of last modification. */
  59.     private static final long serialVersionUID = 20150311L;

  61.     /** the gravitational parameter of the ellipsoid, in m<sup>3</sup>/s<sup>2</sup>. */
  62.     private final double GM;
  63.     /** the rotation rate of the ellipsoid, in rad/s. */
  64.     private final double spin;

  66.     /**
  67.      * Creates a new geodetic Reference Ellipsoid from four defining
  68.      * parameters.
  69.      *
  70.      * @param ae        Equatorial radius, in m
  71.      * @param f         flattening of the ellipsoid.
  72.      * @param bodyFrame the frame to attach to the ellipsoid. The origin is at
  73.      *                  the center of mass, the z axis is the minor axis.
  74.      * @param GM        gravitational parameter, in m<sup>3</sup>/s<sup>2</sup>
  75.      * @param spin      &omega; in rad/s
  76.      */
  77.     public ReferenceEllipsoid(final double ae,
  78.                               final double f,
  79.                               final Frame bodyFrame,
  80.                               final double GM,
  81.                               final double spin) {
  82.         super(ae, f, bodyFrame);
  83.         this.GM = GM;
  84.         this.spin = spin;
  85.     }

  87.     /**
  88.      * Gets the gravitational parameter that is part of the definition of the
  89.      * reference ellipsoid.
  90.      *
  91.      * @return GM in m<sup>3</sup>/s<sup>2</sup>
  92.      */
  93.     public double getGM() {
  94.         return this.GM;
  95.     }

  97.     /**
  98.      * Gets the rotation of the ellipsoid about its axis.
  99.      *
  100.      * @return &omega; in rad/s
  101.      */
  102.     public double getSpin() {
  103.         return this.spin;
  104.     }

  106.     /**
  107.      * Get the radius of this ellipsoid at the poles.
  108.      *
  109.      * @return the polar radius, in meters
  110.      * @see #getEquatorialRadius()
  111.      */
  112.     public double getPolarRadius() {
  113.         // use the definition of flattening: f = (a-b)/a
  114.         final double a = this.getEquatorialRadius();
  115.         final double f = this.getFlattening();
  116.         return a - f * a;
  117.     }

  119.     /**
  120.      * Gets the normal gravity, that is gravity just due to the reference
  121.      * ellipsoid's potential. The normal gravity only depends on latitude
  122.      * because the ellipsoid is axis symmetric.
  123.      *
  124.      * <p> The normal gravity is a vector, having both magnitude and direction.
  125.      * This method only give the magnitude.
  126.      *
  127.      * @param latitude geodetic latitude, in radians. That is the angle between
  128.      *                 the local normal on the ellipsoid and the equatorial
  129.      *                 plane.
  130.      * @return the normal gravity, &gamma;, at the given latitude in
  131.      * m/s<sup>2</sup>. This is the acceleration felt by a mass at rest on the
  132.      * surface of the reference ellipsoid.
  133.      */
  134.     public double getNormalGravity(final double latitude) {
  135.         /*
  136.          * Uses the equations from [2] as compiled in [1]. See Class comment.
  137.          */

  139.         final double a  = this.getEquatorialRadius();
  140.         final double f  = this.getFlattening();

  142.         // define derived constants, move to constructor for more speed
  143.         // semi-minor axis
  144.         final double b = a * (1 - f);
  145.         final double a2 = a * a;
  146.         final double b2 = b * b;
  147.         // linear eccentricity
  148.         final double E = FastMath.sqrt(a2 - b2);
  149.         // first numerical eccentricity
  150.         final double e = E / a;
  151.         // second numerical eccentricity
  152.         final double eprime = E / b;
  153.         // an abbreviation for a common term
  154.         final double m = this.spin * this.spin * a2 * b / this.GM;
  155.         // gravity at equator
  156.         final double ya = this.GM / (a * b) *
  157.                 (1 - 3. / 2. * m - 3. / 14. * eprime * m);
  158.         // gravity at the poles
  159.         final double yb = this.GM / a2 * (1 + m + 3. / 7. * eprime * m);
  160.         // another abbreviation for a common term
  161.         final double kappa = (b * yb - a * ya) / (a * ya);

  163.         // calculate normal gravity at the given latitude.
  164.         final double sin  = FastMath.sin(latitude);
  165.         final double sin2 = sin * sin;
  166.         return ya * (1 + kappa * sin2) / FastMath.sqrt(1 - e * e * sin2);
  167.     }

  169.     /**
  170.      * Get the fully normalized coefficient C<sub>2n,0</sub> for the normal
  171.      * gravity potential.
  172.      *
  173.      * @param n index in C<sub>2n,0</sub>, n &gt;= 1.
  174.      * @return normalized C<sub>2n,0</sub> of the ellipsoid
  175.      * @see "Department of Defense World Geodetic System 1984. 2000. NIMA TR
  176.      * 8350.2 Third Edition, Amendment 1."
  177.      * @see "DMA TR 8350.2. 1984."
  178.      */
  179.     public double getC2n0(final int n) {
  180.         // parameter check
  181.         if (n < 1) {
  182.             throw new IllegalArgumentException("Expected n < 1, got n=" + n);
  183.         }

  185.         final double a = this.getEquatorialRadius();
  186.         final double f = this.getFlattening();
  187.         // define derived constants, move to constructor for more speed
  188.         // semi-minor axis
  189.         final double b = a * (1 - f);
  190.         final double a2 = a * a;
  191.         final double b2 = b * b;
  192.         // linear eccentricity
  193.         final double E = FastMath.sqrt(a2 - b2);
  194.         // first numerical eccentricity
  195.         final double e = E / a;
  196.         // an abbreviation for a common term
  197.         final double m = this.spin * this.spin * a2 * b / this.GM;

  199.         /*
  200.          * derive C2 using a linear approximation, good to ~1e-9, eq 2.118 in
  201.          * Heiskanen & Moritz[2]. See comment for ReferenceEllipsoid
  202.          */
  203.         final double J2 = 2. / 3. * f - 1. / 3. * m - 1. / 3. * f * f + 2. / 21. * f * m;
  204.         final double C2 = -J2 / FastMath.sqrt(5);

  206.         // eq 3-62 in chapter 3 of DMA TR 8350.2, calculated by scaling C2,0
  207.         return (((n & 0x1) == 0) ? 3 : -3) * FastMath.pow(e, 2 * n) *
  208.                 (1 - n - FastMath.pow(5, 3. / 2.) * n * C2 / (e * e)) /
  209.                 ((2 * n + 1) * (2 * n + 3) * FastMath.sqrt(4 * n + 1));
  210.     }

  212.     @Override
  213.     public ReferenceEllipsoid getEllipsoid() {
  214.         return this;
  215.     }

  217.     /**
  218.      * Get the WGS84 ellipsoid, attached to the given body frame.
  219.      *
  220.      * @param bodyFrame the earth centered fixed frame
  221.      * @return a WGS84 reference ellipsoid
  222.      */
  223.     public static ReferenceEllipsoid getWgs84(final Frame bodyFrame) {
  224.         return new ReferenceEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
  225.                 Constants.WGS84_EARTH_FLATTENING, bodyFrame,
  226.                 Constants.WGS84_EARTH_MU,
  227.                 Constants.WGS84_EARTH_ANGULAR_VELOCITY);
  228.     }

  230.     /**
  231.      * Get the GRS80 ellipsoid, attached to the given body frame.
  232.      *
  233.      * @param bodyFrame the earth centered fixed frame
  234.      * @return a GRS80 reference ellipsoid
  235.      */
  236.     public static ReferenceEllipsoid getGrs80(final Frame bodyFrame) {
  237.         return new ReferenceEllipsoid(
  238.                 Constants.GRS80_EARTH_EQUATORIAL_RADIUS,
  239.                 Constants.GRS80_EARTH_FLATTENING,
  240.                 bodyFrame,
  241.                 Constants.GRS80_EARTH_MU,
  242.                 Constants.GRS80_EARTH_ANGULAR_VELOCITY
  243.         );
  244.     }

  246.     /**
  247.      * Get the IERS96 ellipsoid, attached to the given body frame.
  248.      *
  249.      * @param bodyFrame the earth centered fixed frame
  250.      * @return an IERS96 reference ellipsoid
  251.      */
  252.     public static ReferenceEllipsoid getIers96(final Frame bodyFrame) {
  253.         return new ReferenceEllipsoid(Constants.IERS96_EARTH_EQUATORIAL_RADIUS,
  254.                 Constants.IERS96_EARTH_FLATTENING, bodyFrame,
  255.                 Constants.IERS96_EARTH_MU,
  256.                 Constants.IERS96_EARTH_ANGULAR_VELOCITY);
  257.     }

  259.     /**
  260.      * Get the IERS2003 ellipsoid, attached to the given body frame.
  261.      *
  262.      * @param bodyFrame the earth centered fixed frame
  263.      * @return an IERS2003 reference ellipsoid
  264.      */
  265.     public static ReferenceEllipsoid getIers2003(final Frame bodyFrame) {
  266.         return new ReferenceEllipsoid(Constants.IERS2003_EARTH_EQUATORIAL_RADIUS,
  267.                 Constants.IERS2003_EARTH_FLATTENING, bodyFrame,
  268.                 Constants.IERS2003_EARTH_MU,
  269.                 Constants.IERS2003_EARTH_ANGULAR_VELOCITY);
  270.     }

  272.     /**
  273.      * Get the IERS2010 ellipsoid, attached to the given body frame.
  274.      *
  275.      * @param bodyFrame the earth centered fixed frame
  276.      * @return an IERS2010 reference ellipsoid
  277.      */
  278.     public static ReferenceEllipsoid getIers2010(final Frame bodyFrame) {
  279.         return new ReferenceEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
  280.                 Constants.IERS2010_EARTH_FLATTENING, bodyFrame,
  281.                 Constants.IERS2010_EARTH_MU,
  282.                 Constants.IERS2010_EARTH_ANGULAR_VELOCITY);
  283.     }
  284. }
Created with JaCoCo 0.8.12.202403310830