ReferenceEllipsoid.java

  1. /* Contributed in the public domain.
  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.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="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.  */
  55. public class ReferenceEllipsoid extends OneAxisEllipsoid implements EarthShape {

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  245. }
Created with JaCoCo 0.7.9.201702052155