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 */
009package org.jscience.geography.coordinates.crs;
010
011import javax.measure.Measure;
012import javax.measure.Measurable;
013import javax.measure.quantity.*;
014import 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 */
045public class ReferenceEllipsoid {
046
047    /**
048     * The World Geodetic System 1984 reference ellipsoid.
049     */
050    public static final ReferenceEllipsoid WGS84
051        = new ReferenceEllipsoid(6378137.0, 298.257223563);
052    /**
053     * Geodetic Reference System 1980 ellipsoid.
054     */
055    public static final ReferenceEllipsoid GRS80
056        = new ReferenceEllipsoid(6378137.0, 298.257222101);
057    /**
058     * The World Geodetic System 1972 reference ellipsoid.
059     */
060    public static final ReferenceEllipsoid WGS72
061        = new ReferenceEllipsoid(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
067        = new ReferenceEllipsoid(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, double inverseFlattening) {
094        this.a = semimajorAxis;
095        f = 1.0 / inverseFlattening;
096        b = semimajorAxis * (1.0 - f);
097        ea2 = f * (2.0 - f);
098        e = Math.sqrt(ea2);
099        eb2 = ea2 / (1.0 - ea2);
100    }
101
102    private static double sqr(final double x) {
103        return x * x;
104    }
105
106    /**
107     * Returns the semimajor or equatorial radius of this reference ellipsoid.
108     *
109     * @return The semimajor radius.
110     */
111    public Measurable<Length> getSemimajorAxis() {
112        if (_semimajorAxis == null) {
113            _semimajorAxis = Measure.valueOf(a, SI.METRE);
114        }
115        return _semimajorAxis;
116    }
117
118    /**
119     * Returns the semiminor or polar radius of this reference ellipsoid.
120     *
121     * @return  The semiminor radius.
122     */
123    public Measurable<Length> getsSemiminorAxis() {
124        if (_semiminorAxis == null) {
125            _semiminorAxis = Measure.valueOf(b, SI.METRE);
126        }
127        return _semiminorAxis;
128    }
129
130    /**
131     * Returns the flattening or ellipticity of this reference ellipsoid.
132     *
133     * @return The flattening.
134     */
135    public double getFlattening() {
136        return f;
137    }
138
139    /**
140     * Returns the (first) eccentricity of this reference ellipsoid.
141     *
142     * @return The eccentricity.
143     */
144    public double getEccentricity() {
145        return e;
146    }
147
148    /**
149     * Returns the square of the (first) eccentricity. This number is frequently
150     * used in ellipsoidal calculations.
151     *
152     * @return The square of the eccentricity.
153     */
154    public double getEccentricitySquared() {
155        return ea2;
156    }
157
158    /**
159     * Returns the square of the second eccentricity of this reference ellipsoid.
160     * This number is frequently used in ellipsoidal calculations.
161     *
162     * @return The square of the second eccentricity.
163     */
164    public double getSecondEccentricitySquared() {
165        return eb2;
166    }
167
168    /**
169      * Returns the <i>radius of curvature in the prime vertical</i>
170      * for this reference ellipsoid at the specified latitude.
171      *
172      * @param phi The local latitude (radians).
173      * @return The radius of curvature in the prime vertical (meters).
174      */
175     public double verticalRadiusOfCurvature(final double phi) {
176         return a / Math.sqrt(1.0 - (ea2 * sqr(Math.sin(phi))));
177     }
178
179    /**
180      * Returns the <i>radius of curvature in the prime vertical</i>
181      * for this reference ellipsoid at the specified latitude.
182      *
183      * @param latitude The local latitude.
184      * @return The radius of curvature in the prime vertical.
185      */
186     public Measurable<Length> verticalRadiusOfCurvature(final Measurable<Angle> latitude) {
187         return Measure.valueOf(verticalRadiusOfCurvature(latitude.doubleValue(SI.RADIAN)), SI.METRE);
188     }
189
190    /**
191     *  Returns the <i>radius of curvature in the meridian<i>
192     *  for this reference ellipsoid at the specified latitude.
193     *
194     * @param phi The local latitude (in radians).
195     * @return  The radius of curvature in the meridian (in meters).
196     */
197    public double meridionalRadiusOfCurvature(final double phi) {
198        return verticalRadiusOfCurvature(phi)
199               / (1.0 + eb2 * sqr(Math.cos(phi)));
200    }
201
202    /**
203     *  Returns the <i>radius of curvature in the meridian<i>
204     *  for this reference ellipsoid at the specified latitude.
205     *
206     * @param latitude The local latitude (in radians).
207     * @return  The radius of curvature in the meridian (in meters).
208     */
209    public Measurable<Length> meridionalRadiusOfCurvature(final Measurable<Angle> latitude) {
210        return Measure.valueOf(meridionalRadiusOfCurvature(latitude.doubleValue(SI.RADIAN)), SI.METRE);
211    }
212
213    /**
214     *  Returns the meridional arc, the true meridional distance on the
215     * ellipsoid from the equator to the specified latitude, in meters.
216     *
217     * @param phi   The local latitude (in radians).
218     * @return  The meridional arc (in meters).
219     */
220    public double meridionalArc(final double phi) {
221        final double sin2Phi = Math.sin(2.0 * phi);
222        final double sin4Phi = Math.sin(4.0 * phi);
223        final double sin6Phi = Math.sin(6.0 * phi);
224        final double sin8Phi = Math.sin(8.0 * phi);
225        final double n = f / (2.0 - f);
226        final double n2 = n * n;
227        final double n3 = n2 * n;
228        final double n4 = n3 * n;
229        final double n5 = n4 * n;
230        final double n1n2 = n - n2;
231        final double n2n3 = n2 - n3;
232        final double n3n4 = n3 - n4;
233        final double n4n5 = n4 - n5;
234        final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
235        final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
236        final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
237        final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
238        final double ep = (315.0 / 512.0) * a * (n4n5);
239        return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
240    }
241
242    /**
243     *  Returns the meridional arc, the true meridional distance on the
244     * ellipsoid from the equator to the specified latitude.
245     *
246     * @param latitude   The local latitude.
247     * @return  The meridional arc.
248     */
249    public Measurable<Length> meridionalArc(final Measurable<Angle> latitude) {
250        return Measure.valueOf(meridionalArc(latitude.doubleValue(SI.RADIAN)), SI.METRE);
251    }
252
253}