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;
18
19 import org.hipparchus.util.FastMath;
20 import org.orekit.bodies.OneAxisEllipsoid;
21 import org.orekit.frames.Frame;
22 import org.orekit.utils.Constants;
23
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>ω, 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 {
57
58 /** uid is date of last modification. */
59 private static final long serialVersionUID = 20150311L;
60
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;
65
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 ω 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 }
86
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 }
96
97 /**
98 * Gets the rotation of the ellipsoid about its axis.
99 *
100 * @return ω in rad/s
101 */
102 public double getSpin() {
103 return this.spin;
104 }
105
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 }
118
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, γ, 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 */
138
139 final double a = this.getEquatorialRadius();
140 final double f = this.getFlattening();
141
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);
162
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 }
168
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 >= 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 }
184
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;
198
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);
205
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 }
211
212 @Override
213 public ReferenceEllipsoid getEllipsoid() {
214 return this;
215 }
216
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 }
229
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 }
245
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 }
258
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 }
271
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 }