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 /** the gravitational parameter of the ellipsoid, in m<sup>3</sup>/s<sup>2</sup>. */
59 private final double GM;
60 /** the rotation rate of the ellipsoid, in rad/s. */
61 private final double spin;
62
63 /**
64 * Creates a new geodetic Reference Ellipsoid from four defining
65 * parameters.
66 *
67 * @param ae Equatorial radius, in m
68 * @param f flattening of the ellipsoid.
69 * @param bodyFrame the frame to attach to the ellipsoid. The origin is at
70 * the center of mass, the z axis is the minor axis.
71 * @param GM gravitational parameter, in m<sup>3</sup>/s<sup>2</sup>
72 * @param spin ω in rad/s
73 */
74 public ReferenceEllipsoid(final double ae,
75 final double f,
76 final Frame bodyFrame,
77 final double GM,
78 final double spin) {
79 super(ae, f, bodyFrame);
80 this.GM = GM;
81 this.spin = spin;
82 }
83
84 /**
85 * Gets the gravitational parameter that is part of the definition of the
86 * reference ellipsoid.
87 *
88 * @return GM in m<sup>3</sup>/s<sup>2</sup>
89 */
90 public double getGM() {
91 return this.GM;
92 }
93
94 /**
95 * Gets the rotation of the ellipsoid about its axis.
96 *
97 * @return ω in rad/s
98 */
99 public double getSpin() {
100 return this.spin;
101 }
102
103 /**
104 * Get the radius of this ellipsoid at the poles.
105 *
106 * @return the polar radius, in meters
107 * @see #getEquatorialRadius()
108 */
109 public double getPolarRadius() {
110 // use the definition of flattening: f = (a-b)/a
111 final double a = this.getEquatorialRadius();
112 final double f = this.getFlattening();
113 return a - f * a;
114 }
115
116 /**
117 * Gets the normal gravity, that is gravity just due to the reference
118 * ellipsoid's potential. The normal gravity only depends on latitude
119 * because the ellipsoid is axis symmetric.
120 *
121 * <p> The normal gravity is a vector, having both magnitude and direction.
122 * This method only give the magnitude.
123 *
124 * @param latitude geodetic latitude, in radians. That is the angle between
125 * the local normal on the ellipsoid and the equatorial
126 * plane.
127 * @return the normal gravity, γ, at the given latitude in
128 * m/s<sup>2</sup>. This is the acceleration felt by a mass at rest on the
129 * surface of the reference ellipsoid.
130 */
131 public double getNormalGravity(final double latitude) {
132 /*
133 * Uses the equations from [2] as compiled in [1]. See Class comment.
134 */
135
136 final double a = this.getEquatorialRadius();
137 final double f = this.getFlattening();
138
139 // define derived constants, move to constructor for more speed
140 // semi-minor axis
141 final double b = a * (1 - f);
142 final double a2 = a * a;
143 final double b2 = b * b;
144 // linear eccentricity
145 final double E = FastMath.sqrt(a2 - b2);
146 // first numerical eccentricity
147 final double e = E / a;
148 // second numerical eccentricity
149 final double eprime = E / b;
150 // an abbreviation for a common term
151 final double m = this.spin * this.spin * a2 * b / this.GM;
152 // gravity at equator
153 final double ya = this.GM / (a * b) *
154 (1 - 3. / 2. * m - 3. / 14. * eprime * m);
155 // gravity at the poles
156 final double yb = this.GM / a2 * (1 + m + 3. / 7. * eprime * m);
157 // another abbreviation for a common term
158 final double kappa = (b * yb - a * ya) / (a * ya);
159
160 // calculate normal gravity at the given latitude.
161 final double sin = FastMath.sin(latitude);
162 final double sin2 = sin * sin;
163 return ya * (1 + kappa * sin2) / FastMath.sqrt(1 - e * e * sin2);
164 }
165
166 /**
167 * Get the fully normalized coefficient C<sub>2n,0</sub> for the normal
168 * gravity potential.
169 *
170 * @param n index in C<sub>2n,0</sub>, n >= 1.
171 * @return normalized C<sub>2n,0</sub> of the ellipsoid
172 * @see "Department of Defense World Geodetic System 1984. 2000. NIMA TR
173 * 8350.2 Third Edition, Amendment 1."
174 * @see "DMA TR 8350.2. 1984."
175 */
176 public double getC2n0(final int n) {
177 // parameter check
178 if (n < 1) {
179 throw new IllegalArgumentException("Expected n < 1, got n=" + n);
180 }
181
182 final double a = this.getEquatorialRadius();
183 final double f = this.getFlattening();
184 // define derived constants, move to constructor for more speed
185 // semi-minor axis
186 final double b = a * (1 - f);
187 final double a2 = a * a;
188 final double b2 = b * b;
189 // linear eccentricity
190 final double E = FastMath.sqrt(a2 - b2);
191 // first numerical eccentricity
192 final double e = E / a;
193 // an abbreviation for a common term
194 final double m = this.spin * this.spin * a2 * b / this.GM;
195
196 /*
197 * derive C2 using a linear approximation, good to ~1e-9, eq 2.118 in
198 * Heiskanen & Moritz[2]. See comment for ReferenceEllipsoid
199 */
200 final double J2 = 2. / 3. * f - 1. / 3. * m - 1. / 3. * f * f + 2. / 21. * f * m;
201 final double C2 = -J2 / FastMath.sqrt(5);
202
203 // eq 3-62 in chapter 3 of DMA TR 8350.2, calculated by scaling C2,0
204 return (((n & 0x1) == 0) ? 3 : -3) * FastMath.pow(e, 2 * n) *
205 (1 - n - FastMath.pow(5, 3. / 2.) * n * C2 / (e * e)) /
206 ((2 * n + 1) * (2 * n + 3) * FastMath.sqrt(4 * n + 1));
207 }
208
209 @Override
210 public ReferenceEllipsoid getEllipsoid() {
211 return this;
212 }
213
214 /**
215 * Get the WGS84 ellipsoid, attached to the given body frame.
216 *
217 * @param bodyFrame the earth centered fixed frame
218 * @return a WGS84 reference ellipsoid
219 */
220 public static ReferenceEllipsoid getWgs84(final Frame bodyFrame) {
221 return new ReferenceEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
222 Constants.WGS84_EARTH_FLATTENING, bodyFrame,
223 Constants.WGS84_EARTH_MU,
224 Constants.WGS84_EARTH_ANGULAR_VELOCITY);
225 }
226
227 /**
228 * Get the GRS80 ellipsoid, attached to the given body frame.
229 *
230 * @param bodyFrame the earth centered fixed frame
231 * @return a GRS80 reference ellipsoid
232 */
233 public static ReferenceEllipsoid getGrs80(final Frame bodyFrame) {
234 return new ReferenceEllipsoid(
235 Constants.GRS80_EARTH_EQUATORIAL_RADIUS,
236 Constants.GRS80_EARTH_FLATTENING,
237 bodyFrame,
238 Constants.GRS80_EARTH_MU,
239 Constants.GRS80_EARTH_ANGULAR_VELOCITY
240 );
241 }
242
243 /**
244 * Get the IERS96 ellipsoid, attached to the given body frame.
245 *
246 * @param bodyFrame the earth centered fixed frame
247 * @return an IERS96 reference ellipsoid
248 */
249 public static ReferenceEllipsoid getIers96(final Frame bodyFrame) {
250 return new ReferenceEllipsoid(Constants.IERS96_EARTH_EQUATORIAL_RADIUS,
251 Constants.IERS96_EARTH_FLATTENING, bodyFrame,
252 Constants.IERS96_EARTH_MU,
253 Constants.IERS96_EARTH_ANGULAR_VELOCITY);
254 }
255
256 /**
257 * Get the IERS2003 ellipsoid, attached to the given body frame.
258 *
259 * @param bodyFrame the earth centered fixed frame
260 * @return an IERS2003 reference ellipsoid
261 */
262 public static ReferenceEllipsoid getIers2003(final Frame bodyFrame) {
263 return new ReferenceEllipsoid(Constants.IERS2003_EARTH_EQUATORIAL_RADIUS,
264 Constants.IERS2003_EARTH_FLATTENING, bodyFrame,
265 Constants.IERS2003_EARTH_MU,
266 Constants.IERS2003_EARTH_ANGULAR_VELOCITY);
267 }
268
269 /**
270 * Get the IERS2010 ellipsoid, attached to the given body frame.
271 *
272 * @param bodyFrame the earth centered fixed frame
273 * @return an IERS2010 reference ellipsoid
274 */
275 public static ReferenceEllipsoid getIers2010(final Frame bodyFrame) {
276 return new ReferenceEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
277 Constants.IERS2010_EARTH_FLATTENING, bodyFrame,
278 Constants.IERS2010_EARTH_MU,
279 Constants.IERS2010_EARTH_ANGULAR_VELOCITY);
280 }
281 }