001: /*
002: * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
003: * Copyright (C) 2006 - JScience (http://jscience.org/)
004: * All rights reserved.
005: *
006: * Permission to use, copy, modify, and distribute this software is
007: * freely granted, provided that this notice is preserved.
008: */
009: package org.jscience.geography.coordinates.crs;
010:
011: import javax.measure.Measure;
012: import javax.measure.Measurable;
013: import javax.measure.quantity.*;
014: import javax.measure.unit.SI;
015:
016: /**
017: * <p> The ReferenceEllipsoid class defines a geodetic reference ellipsoid
018: * used as a standard for geodetic measurements. The World Geodetic System
019: * 1984 (WGS84) ellipsoid is the current standard for most geographic and
020: * geodetic coordinate systems, including GPS. The WGS84 ellipsoid is
021: * provided as a static instance of this class.</p>
022: *
023: * <p> The ellipsoid (actually an oblate spheroid) is uniquely specified by
024: * two parameters, the semimajor (or equatorial) radius and the ellipticity
025: * or flattening. In practice, the reciprocal of the flattening is
026: * specified.</p>
027: *
028: * <p> The ellipsoid is an approximation of the shape of the earth. Although
029: * not exact, the ellipsoid is much more accurate than a spherical
030: * approximation and is still mathematically simple. The <i>geoid</i> is
031: * a still closer approximation of the shape of the earth (intended to
032: * represent the mean sea level), and is generally specified by it's
033: * deviation from the ellipsoid.</p>
034: *
035: * <p> Different reference ellipsoids give more or less accurate results at
036: * different locations, so it was previously common for different nations
037: * to use ellipsoids that were more accurate for their areas. More recent
038: * efforts have provided ellipsoids with better overall global accuracy,
039: * such as the WGS84 ellipsiod, and these have now largely supplanted
040: * the others.</p>
041: *
042: * @author Paul D. Anderson
043: * @version 3.0, February 18, 2006
044: */
045: public class ReferenceEllipsoid {
046:
047: /**
048: * The World Geodetic System 1984 reference ellipsoid.
049: */
050: public static final ReferenceEllipsoid WGS84 = new ReferenceEllipsoid(
051: 6378137.0, 298.257223563);
052: /**
053: * Geodetic Reference System 1980 ellipsoid.
054: */
055: public static final ReferenceEllipsoid GRS80 = new ReferenceEllipsoid(
056: 6378137.0, 298.257222101);
057: /**
058: * The World Geodetic System 1972 reference ellipsoid.
059: */
060: public static final ReferenceEllipsoid WGS72 = new ReferenceEllipsoid(
061: 6378135.0, 298.26);
062: /**
063: * The International 1924 reference ellipsoid, one of the earliest
064: * "global" ellipsoids.
065: */
066: public static final ReferenceEllipsoid INTERNATIONAL1924 = new ReferenceEllipsoid(
067: 6378388.0, 297.0);
068:
069: private final double a;
070:
071: private final double b;
072:
073: private final double f;
074:
075: private final double ea2;
076:
077: private final double e;
078:
079: private final double eb2;
080:
081: private Measurable<Length> _semimajorAxis;
082:
083: private Measurable<Length> _semiminorAxis;
084:
085: /**
086: * Constructs an instance of a reference ellipsoid.
087: *
088: * @param semimajorAxis The semimajor or equatorial radius of this
089: * reference ellipsoid, in meters.
090: * @param inverseFlattening The reciprocal of the ellipticity or flattening
091: * of this reference ellipsoid (dimensionless).
092: */
093: public ReferenceEllipsoid(double semimajorAxis,
094: double inverseFlattening) {
095: this .a = semimajorAxis;
096: f = 1.0 / inverseFlattening;
097: b = semimajorAxis * (1.0 - f);
098: ea2 = f * (2.0 - f);
099: e = Math.sqrt(ea2);
100: eb2 = ea2 / (1.0 - ea2);
101: }
102:
103: private static double sqr(final double x) {
104: return x * x;
105: }
106:
107: /**
108: * Returns the semimajor or equatorial radius of this reference ellipsoid.
109: *
110: * @return The semimajor radius.
111: */
112: public Measurable<Length> getSemimajorAxis() {
113: if (_semimajorAxis == null) {
114: _semimajorAxis = Measure.valueOf(a, SI.METRE);
115: }
116: return _semimajorAxis;
117: }
118:
119: /**
120: * Returns the semiminor or polar radius of this reference ellipsoid.
121: *
122: * @return The semiminor radius.
123: */
124: public Measurable<Length> getsSemiminorAxis() {
125: if (_semiminorAxis == null) {
126: _semiminorAxis = Measure.valueOf(b, SI.METRE);
127: }
128: return _semiminorAxis;
129: }
130:
131: /**
132: * Returns the flattening or ellipticity of this reference ellipsoid.
133: *
134: * @return The flattening.
135: */
136: public double getFlattening() {
137: return f;
138: }
139:
140: /**
141: * Returns the (first) eccentricity of this reference ellipsoid.
142: *
143: * @return The eccentricity.
144: */
145: public double getEccentricity() {
146: return e;
147: }
148:
149: /**
150: * Returns the square of the (first) eccentricity. This number is frequently
151: * used in ellipsoidal calculations.
152: *
153: * @return The square of the eccentricity.
154: */
155: public double getEccentricitySquared() {
156: return ea2;
157: }
158:
159: /**
160: * Returns the square of the second eccentricity of this reference ellipsoid.
161: * This number is frequently used in ellipsoidal calculations.
162: *
163: * @return The square of the second eccentricity.
164: */
165: public double getSecondEccentricitySquared() {
166: return eb2;
167: }
168:
169: /**
170: * Returns the <i>radius of curvature in the prime vertical</i>
171: * for this reference ellipsoid at the specified latitude.
172: *
173: * @param phi The local latitude (radians).
174: * @return The radius of curvature in the prime vertical (meters).
175: */
176: public double verticalRadiusOfCurvature(final double phi) {
177: return a / Math.sqrt(1.0 - (ea2 * sqr(Math.sin(phi))));
178: }
179:
180: /**
181: * Returns the <i>radius of curvature in the prime vertical</i>
182: * for this reference ellipsoid at the specified latitude.
183: *
184: * @param latitude The local latitude.
185: * @return The radius of curvature in the prime vertical.
186: */
187: public Measurable<Length> verticalRadiusOfCurvature(
188: final Measurable<Angle> latitude) {
189: return Measure.valueOf(verticalRadiusOfCurvature(latitude
190: .doubleValue(SI.RADIAN)), SI.METRE);
191: }
192:
193: /**
194: * Returns the <i>radius of curvature in the meridian<i>
195: * for this reference ellipsoid at the specified latitude.
196: *
197: * @param phi The local latitude (in radians).
198: * @return The radius of curvature in the meridian (in meters).
199: */
200: public double meridionalRadiusOfCurvature(final double phi) {
201: return verticalRadiusOfCurvature(phi)
202: / (1.0 + eb2 * sqr(Math.cos(phi)));
203: }
204:
205: /**
206: * Returns the <i>radius of curvature in the meridian<i>
207: * for this reference ellipsoid at the specified latitude.
208: *
209: * @param latitude The local latitude (in radians).
210: * @return The radius of curvature in the meridian (in meters).
211: */
212: public Measurable<Length> meridionalRadiusOfCurvature(
213: final Measurable<Angle> latitude) {
214: return Measure.valueOf(meridionalRadiusOfCurvature(latitude
215: .doubleValue(SI.RADIAN)), SI.METRE);
216: }
217:
218: /**
219: * Returns the meridional arc, the true meridional distance on the
220: * ellipsoid from the equator to the specified latitude, in meters.
221: *
222: * @param phi The local latitude (in radians).
223: * @return The meridional arc (in meters).
224: */
225: public double meridionalArc(final double phi) {
226: final double sin2Phi = Math.sin(2.0 * phi);
227: final double sin4Phi = Math.sin(4.0 * phi);
228: final double sin6Phi = Math.sin(6.0 * phi);
229: final double sin8Phi = Math.sin(8.0 * phi);
230: final double n = f / (2.0 - f);
231: final double n2 = n * n;
232: final double n3 = n2 * n;
233: final double n4 = n3 * n;
234: final double n5 = n4 * n;
235: final double n1n2 = n - n2;
236: final double n2n3 = n2 - n3;
237: final double n3n4 = n3 - n4;
238: final double n4n5 = n4 - n5;
239: final double ap = a
240: * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0)
241: * (n4n5));
242: final double bp = (3.0 / 2.0) * a
243: * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
244: final double cp = (15.0 / 16.0) * a
245: * (n2n3 + (3.0 / 4.0) * (n4n5));
246: final double dp = (35.0 / 48.0) * a
247: * (n3n4 + (11.0 / 16.0) * n5);
248: final double ep = (315.0 / 512.0) * a * (n4n5);
249: return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi
250: + ep * sin8Phi;
251: }
252:
253: /**
254: * Returns the meridional arc, the true meridional distance on the
255: * ellipsoid from the equator to the specified latitude.
256: *
257: * @param latitude The local latitude.
258: * @return The meridional arc.
259: */
260: public Measurable<Length> meridionalArc(
261: final Measurable<Angle> latitude) {
262: return Measure.valueOf(meridionalArc(latitude
263: .doubleValue(SI.RADIAN)), SI.METRE);
264: }
265:
266: }
|