001/*
002 * Import from fr.geo.convert package, a geographic coordinates converter.
003 * (https://www.i3s.unice.fr/~johan/gps/)
004 * License: GPL. For details, see LICENSE file.
005 * Copyright (C) 2002 Johan Montagnat (johan@creatis.insa-lyon.fr)
006 */
007package org.openstreetmap.josm.data.projection;
008
009import org.openstreetmap.josm.data.coor.LatLon;
010
011/**
012 * Reference ellipsoids.
013 */
014public final class Ellipsoid {
015
016    /**
017     * Airy 1830
018     */
019    public static final Ellipsoid Airy = Ellipsoid.createAb(6377563.396, 6356256.910);
020
021    /**
022     * Modified Airy 1849
023     */
024    public static final Ellipsoid AiryMod = Ellipsoid.createAb(6377340.189, 6356034.446);
025
026    /**
027     * Australian National Spheroid (Australian Natl & S. Amer. 1969)
028     * same as GRS67 Modified
029     */
030    public static final Ellipsoid AustSA = Ellipsoid.createArf(6378160.0, 298.25);
031
032    /**
033     * Bessel 1841 ellipsoid
034     */
035    public static final Ellipsoid Bessel1841 = Ellipsoid.createArf(6377397.155, 299.1528128);
036
037    /**
038     * Bessel 1841 (Namibia)
039     */
040    public static final Ellipsoid BesselNamibia = Ellipsoid.createArf(6377483.865, 299.1528128);
041
042    /**
043     * Clarke 1866 ellipsoid
044     */
045    public static final Ellipsoid Clarke1866 = Ellipsoid.createAb(6378206.4, 6356583.8);
046
047    /**
048     * Clarke 1880 (modified)
049     */
050    public static final Ellipsoid Clarke1880 = Ellipsoid.createArf(6378249.145, 293.4663);
051
052    /**
053     * Clarke 1880 IGN (French national geographic institute)
054     */
055    public static final Ellipsoid ClarkeIGN = Ellipsoid.createAb(6378249.2, 6356515.0);
056
057    /**
058     * Everest (Sabah & Sarawak)
059     */
060    public static final Ellipsoid EverestSabahSarawak = Ellipsoid.createArf(6377298.556, 300.8017);
061
062    /**
063     * GRS67 ellipsoid
064     */
065    public static final Ellipsoid GRS67 = Ellipsoid.createArf(6378160.0, 298.247167427);
066
067    /**
068     * GRS80 ellipsoid
069     */
070    public static final Ellipsoid GRS80 = Ellipsoid.createArf(6378137.0, 298.257222101);
071
072    /**
073     * Hayford's ellipsoid 1909 (ED50 system)
074     * Also known as International 1924
075     * Proj.4 code: intl
076     */
077    public static final Ellipsoid Hayford = Ellipsoid.createArf(6378388.0, 297.0);
078
079    /**
080     * Helmert 1906
081     */
082    public static final Ellipsoid Helmert = Ellipsoid.createArf(6378200.0, 298.3);
083
084    /**
085     * Krassowsky 1940 ellipsoid
086     */
087    public static final Ellipsoid Krassowsky = Ellipsoid.createArf(6378245.0, 298.3);
088
089    /**
090     * WGS66 ellipsoid
091     */
092    public static final Ellipsoid WGS66 = Ellipsoid.createArf(6378145.0, 298.25);
093
094    /**
095     * WGS72 ellipsoid
096     */
097    public static final Ellipsoid WGS72 = Ellipsoid.createArf(6378135.0, 298.26);
098
099    /**
100     * WGS84 ellipsoid
101     */
102    public static final Ellipsoid WGS84 = Ellipsoid.createArf(6378137.0, 298.257223563);
103
104    /**
105     * half long axis
106     */
107    public final double a;
108
109    /**
110     * half short axis
111     */
112    public final double b;
113
114    /**
115     * first eccentricity:
116     * sqrt(a*a - b*b) / a
117     */
118    public final double e;
119
120    /**
121     * first eccentricity squared:
122     * (a*a - b*b) / (a*a)
123     */
124    public final double e2;
125
126    /**
127     * square of the second eccentricity:
128     * (a*a - b*b) / (b*b)
129     */
130    public final double eb2;
131
132    /**
133     * if ellipsoid is spherical, i.e. the major and minor semiaxis are
134     * the same
135     */
136    public final boolean spherical;
137
138    /**
139     * private constructur - use one of the create_* methods
140     *
141     * @param a semimajor radius of the ellipsoid axis
142     * @param b semiminor radius of the ellipsoid axis
143     * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a)))
144     * @param e2 first eccentricity squared
145     * @param eb2 square of the second eccentricity
146     * @param sperical if the ellipsoid is sphere
147     */
148    private Ellipsoid(double a, double b, double e, double e2, double eb2, boolean sperical) {
149        this.a = a;
150        this.b = b;
151        this.e = e;
152        this.e2 = e2;
153        this.eb2 = eb2;
154        this.spherical = sperical;
155    }
156
157    /**
158     * create a new ellipsoid
159     *
160     * @param a semimajor radius of the ellipsoid axis (in meters)
161     * @param b semiminor radius of the ellipsoid axis (in meters)
162     * @return the new ellipsoid
163     */
164    public static Ellipsoid createAb(double a, double b) {
165        double e2 = (a*a - b*b) / (a*a);
166        double e = Math.sqrt(e2);
167        double eb2 = e2 / (1.0 - e2);
168        return new Ellipsoid(a, b, e, e2, eb2, a == b);
169    }
170
171    /**
172     * create a new ellipsoid
173     *
174     * @param a semimajor radius of the ellipsoid axis (in meters)
175     * @param es first eccentricity squared
176     * @return the new ellipsoid
177     */
178    public static Ellipsoid createAes(double a, double es) {
179        double b = a * Math.sqrt(1.0 - es);
180        double e = Math.sqrt(es);
181        double eb2 = es / (1.0 - es);
182        return new Ellipsoid(a, b, e, es, eb2, es == 0);
183    }
184
185    /**
186     * create a new ellipsoid
187     *
188     * @param a semimajor radius of the ellipsoid axis (in meters)
189     * @param f flattening ( = (a - b) / a)
190     * @return the new ellipsoid
191     */
192    public static Ellipsoid createAf(double a, double f) {
193        double b = a * (1.0 - f);
194        double e2 = f * (2 - f);
195        double e = Math.sqrt(e2);
196        double eb2 = e2 / (1.0 - e2);
197        return new Ellipsoid(a, b, e, e2, eb2, f == 0);
198    }
199
200    /**
201     * create a new ellipsoid
202     *
203     * @param a semimajor radius of the ellipsoid axis (in meters)
204     * @param rf inverse flattening
205     * @return the new ellipsoid
206     */
207    public static Ellipsoid createArf(double a, double rf) {
208        return createAf(a, 1.0 / rf);
209    }
210
211    @Override
212    public String toString() {
213        return "Ellipsoid{a="+a+", b="+b+'}';
214    }
215
216    /**
217     * Returns the <i>radius of curvature in the prime vertical</i>
218     * for this reference ellipsoid at the specified latitude.
219     *
220     * @param phi The local latitude (radians).
221     * @return The radius of curvature in the prime vertical (meters).
222     */
223    public double verticalRadiusOfCurvature(final double phi) {
224        return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
225    }
226
227    private static double sqr(final double x) {
228        return x * x;
229    }
230
231    /**
232     *  Returns the meridional arc, the true meridional distance on the
233     * ellipsoid from the equator to the specified latitude, in meters.
234     *
235     * @param phi   The local latitude (in radians).
236     * @return  The meridional arc (in meters).
237     */
238    public double meridionalArc(final double phi) {
239        final double sin2Phi = Math.sin(2.0 * phi);
240        final double sin4Phi = Math.sin(4.0 * phi);
241        final double sin6Phi = Math.sin(6.0 * phi);
242        final double sin8Phi = Math.sin(8.0 * phi);
243        // TODO . calculate 'f'
244        //double f = 1.0 / 298.257222101; // GRS80
245        double f = 1.0 / 298.257223563; // WGS84
246        final double n = f / (2.0 - f);
247        final double n2 = n * n;
248        final double n3 = n2 * n;
249        final double n4 = n3 * n;
250        final double n5 = n4 * n;
251        final double n1n2 = n - n2;
252        final double n2n3 = n2 - n3;
253        final double n3n4 = n3 - n4;
254        final double n4n5 = n4 - n5;
255        final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
256        final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
257        final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
258        final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
259        final double ep = (315.0 / 512.0) * a * (n4n5);
260        return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
261    }
262
263    /**
264     *  Returns the <i>radius of curvature in the meridian</i>
265     *  for this reference ellipsoid at the specified latitude.
266     *
267     * @param phi The local latitude (in radians).
268     * @return  The radius of curvature in the meridian (in meters).
269     */
270    public double meridionalRadiusOfCurvature(final double phi) {
271        return verticalRadiusOfCurvature(phi)
272        / (1.0 + eb2 * sqr(Math.cos(phi)));
273    }
274
275    /**
276     * Returns isometric latitude of phi on given first eccentricity (e)
277     * @param phi The local latitude (radians).
278     * @param e first eccentricity
279     * @return isometric latitude of phi on first eccentricity (e)
280     */
281    public double latitudeIsometric(double phi, double e) {
282        double v1 = 1-e*Math.sin(phi);
283        double v2 = 1+e*Math.sin(phi);
284        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
285    }
286
287    /**
288     * Returns isometric latitude of phi on first eccentricity (e)
289     * @param phi The local latitude (radians).
290     * @return isometric latitude of phi on first eccentricity (e)
291     */
292    public double latitudeIsometric(double phi) {
293        double v1 = 1-e*Math.sin(phi);
294        double v2 = 1+e*Math.sin(phi);
295        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
296    }
297
298    /**
299     * Returns geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision
300     * @param latIso isometric latitude
301     * @param e first eccentricity
302     * @param epsilon epsilon precision
303     * @return geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision
304     */
305    public double latitude(double latIso, double e, double epsilon) {
306        double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
307        double lati = lat0;
308        double lati1 = 1.0; // random value to start the iterative processus
309        while (Math.abs(lati1-lati) >= epsilon) {
310            lati = lati1;
311            double v1 = 1+e*Math.sin(lati);
312            double v2 = 1-e*Math.sin(lati);
313            lati1 = 2*Math.atan(Math.pow(v1/v2, e/2)*Math.exp(latIso))-Math.PI/2;
314        }
315        return lati1;
316    }
317
318    /**
319     * convert cartesian coordinates to ellipsoidal coordinates
320     *
321     * @param xyz the coordinates in meters (X, Y, Z)
322     * @return The corresponding latitude and longitude in degrees
323     */
324    public LatLon cart2LatLon(double ... xyz) {
325        return cart2LatLon(xyz, 1e-11);
326    }
327
328    public LatLon cart2LatLon(double[] xyz, double epsilon) {
329        double norm = Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
330        double lg = 2.0 * Math.atan(xyz[1] / (xyz[0] + norm));
331        double lt = Math.atan(xyz[2] / (norm * (1.0 - (a * e2 / Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2])))));
332        double delta = 1.0;
333        while (delta > epsilon) {
334            double s2 = Math.sin(lt);
335            s2 *= s2;
336            double l = Math.atan((xyz[2] / norm)
337                    / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2)))));
338            delta = Math.abs(l - lt);
339            lt = l;
340        }
341        return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg));
342    }
343
344    /**
345     * convert ellipsoidal coordinates to cartesian coordinates
346     *
347     * @param coord The Latitude and longitude in degrees
348     * @return the corresponding (X, Y Z) cartesian coordinates in meters.
349     */
350    public double[] latLon2Cart(LatLon coord) {
351        double phi = Math.toRadians(coord.lat());
352        double lambda = Math.toRadians(coord.lon());
353
354        double rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2));
355        double[] xyz = new double[3];
356        xyz[0] = rn * Math.cos(phi) * Math.cos(lambda);
357        xyz[1] = rn * Math.cos(phi) * Math.sin(lambda);
358        xyz[2] = rn * (1 - e2) * Math.sin(phi);
359
360        return xyz;
361    }
362}