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