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.create_a_b(6377563.396, 6356256.910);
020
021    /**
022     * Modified Airy 1849
023     */
024    public static final Ellipsoid AiryMod = Ellipsoid.create_a_b(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.create_a_rf(6378160.0, 298.25);
031
032    /**
033     * Bessel 1841 ellipsoid
034     */
035    public static final Ellipsoid Bessel1841 = Ellipsoid.create_a_rf(6377397.155, 299.1528128);
036
037    /**
038     * Clarke 1866 ellipsoid
039     */
040    public static final Ellipsoid Clarke1866 = Ellipsoid.create_a_b(6378206.4, 6356583.8);
041
042    /**
043     * Clarke 1880 IGN (French national geographic institute)
044     */
045    public static final Ellipsoid ClarkeIGN = Ellipsoid.create_a_b(6378249.2, 6356515.0);
046
047    /**
048     * GRS67 ellipsoid
049     */
050    public static final Ellipsoid GRS67 = Ellipsoid.create_a_rf(6378160.0, 298.247167427);
051
052    /**
053     * GRS80 ellipsoid
054     */
055    public static final Ellipsoid GRS80 = Ellipsoid.create_a_rf(6378137.0, 298.257222101);
056
057    /**
058     * Hayford's ellipsoid 1909 (ED50 system)
059     * Also known as International 1924
060     * Proj.4 code: intl
061     */
062    public static final Ellipsoid Hayford = Ellipsoid.create_a_rf(6378388.0, 297.0);
063
064    /**
065     * Helmert 1906
066     */
067    public static final Ellipsoid Helmert = Ellipsoid.create_a_rf(6378200.0, 298.3);
068
069    /**
070     * Krassowsky 1940 ellipsoid
071     */
072    public static final Ellipsoid Krassowsky = Ellipsoid.create_a_rf(6378245.0, 298.3);
073
074    /**
075     * WGS72 ellipsoid
076     */
077    public static final Ellipsoid WGS72 = Ellipsoid.create_a_rf(6378135.0, 298.26);
078
079    /**
080     * WGS84 ellipsoid
081     */
082    public static final Ellipsoid WGS84 = Ellipsoid.create_a_rf(6378137.0, 298.257223563);
083
084
085    /**
086     * half long axis
087     */
088    public final double a;
089
090    /**
091     * half short axis
092     */
093    public final double b;
094
095    /**
096     * first eccentricity
097     */
098    public final double e;
099
100    /**
101     * first eccentricity squared
102     */
103    public final double e2;
104
105    /**
106     * square of the second eccentricity
107     */
108    public final double eb2;
109
110    /**
111     * private constructur - use one of the create_* methods
112     *
113     * @param a semimajor radius of the ellipsoid axis
114     * @param b semiminor radius of the ellipsoid axis
115     * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a)))
116     * @param e2 first eccentricity squared
117     * @param eb2 square of the second eccentricity
118     */
119    private Ellipsoid(double a, double b, double e, double e2, double eb2) {
120        this.a = a;
121        this.b = b;
122        this.e = e;
123        this.e2 = e2;
124        this.eb2 = eb2;
125    }
126
127    /**
128     * create a new ellipsoid
129     *
130     * @param a semimajor radius of the ellipsoid axis (in meters)
131     * @param b semiminor radius of the ellipsoid axis (in meters)
132     * @return the new ellipsoid
133     */
134    public static Ellipsoid create_a_b(double a, double b) {
135        double e2 = (a*a - b*b) / (a*a);
136        double e = Math.sqrt(e2);
137        double eb2 = e2 / (1.0 - e2);
138        return new Ellipsoid(a, b, e, e2, eb2);
139    }
140
141    /**
142     * create a new ellipsoid
143     *
144     * @param a semimajor radius of the ellipsoid axis (in meters)
145     * @param es first eccentricity squared
146     * @return the new ellipsoid
147     */
148    public static Ellipsoid create_a_es(double a, double es) {
149        double b = a * Math.sqrt(1.0 - es);
150        double e = Math.sqrt(es);
151        double eb2 = es / (1.0 - es);
152        return new Ellipsoid(a, b, e, es, eb2);
153    }
154
155    /**
156     * create a new ellipsoid
157     *
158     * @param a semimajor radius of the ellipsoid axis (in meters)
159     * @param f flattening ( = (a - b) / a)
160     * @return the new ellipsoid
161     */
162    public static Ellipsoid create_a_f(double a, double f) {
163        double b = a * (1.0 - f);
164        double e2 = f * (2 - f);
165        double e = Math.sqrt(e2);
166        double eb2 = e2 / (1.0 - e2);
167        return new Ellipsoid(a, b, e, e2, eb2);
168    }
169
170    /**
171     * create a new ellipsoid
172     *
173     * @param a semimajor radius of the ellipsoid axis (in meters)
174     * @param rf inverse flattening
175     * @return the new ellipsoid
176     */
177    public static Ellipsoid create_a_rf(double a, double rf) {
178        return create_a_f(a, 1.0 / rf);
179    }
180
181    @Override
182    public String toString() {
183        return "Ellipsoid{a="+a+", b="+b+'}';
184    }
185
186    /**
187     * Returns the <i>radius of curvature in the prime vertical</i>
188     * for this reference ellipsoid at the specified latitude.
189     *
190     * @param phi The local latitude (radians).
191     * @return The radius of curvature in the prime vertical (meters).
192     */
193    public double verticalRadiusOfCurvature(final double phi) {
194        return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
195    }
196
197    private static double sqr(final double x) {
198        return x * x;
199    }
200
201    /**
202     *  Returns the meridional arc, the true meridional distance on the
203     * ellipsoid from the equator to the specified latitude, in meters.
204     *
205     * @param phi   The local latitude (in radians).
206     * @return  The meridional arc (in meters).
207     */
208    public double meridionalArc(final double phi) {
209        final double sin2Phi = Math.sin(2.0 * phi);
210        final double sin4Phi = Math.sin(4.0 * phi);
211        final double sin6Phi = Math.sin(6.0 * phi);
212        final double sin8Phi = Math.sin(8.0 * phi);
213        // TODO . calculate 'f'
214        //double f = 1.0 / 298.257222101; // GRS80
215        double f = 1.0 / 298.257223563; // WGS84
216        final double n = f / (2.0 - f);
217        final double n2 = n * n;
218        final double n3 = n2 * n;
219        final double n4 = n3 * n;
220        final double n5 = n4 * n;
221        final double n1n2 = n - n2;
222        final double n2n3 = n2 - n3;
223        final double n3n4 = n3 - n4;
224        final double n4n5 = n4 - n5;
225        final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
226        final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
227        final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
228        final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
229        final double ep = (315.0 / 512.0) * a * (n4n5);
230        return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
231    }
232
233    /**
234     *  Returns the <i>radius of curvature in the meridian</i>
235     *  for this reference ellipsoid at the specified latitude.
236     *
237     * @param phi The local latitude (in radians).
238     * @return  The radius of curvature in the meridian (in meters).
239     */
240    public double meridionalRadiusOfCurvature(final double phi) {
241        return verticalRadiusOfCurvature(phi)
242        / (1.0 + eb2 * sqr(Math.cos(phi)));
243    }
244
245    /**
246     * Returns isometric latitude of phi on given first eccentricity (e)
247     * @param phi The local latitude (radians).
248     * @return isometric latitude of phi on first eccentricity (e)
249     */
250    public double latitudeIsometric(double phi, double e) {
251        double v1 = 1-e*Math.sin(phi);
252        double v2 = 1+e*Math.sin(phi);
253        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
254    }
255
256    /**
257     * Returns isometric latitude of phi on first eccentricity (e)
258     * @param phi The local latitude (radians).
259     * @return isometric latitude of phi on first eccentricity (e)
260     */
261    public double latitudeIsometric(double phi) {
262        double v1 = 1-e*Math.sin(phi);
263        double v2 = 1+e*Math.sin(phi);
264        return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
265    }
266
267    /**
268     * Returns geographic latitude of isometric latitude of first eccentricity (e)
269     * and epsilon precision
270     * @return geographic latitude of isometric latitude of first eccentricity (e)
271     * and epsilon precision
272     */
273    public double latitude(double latIso, double e, double epsilon) {
274        double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
275        double lati = lat0;
276        double lati1 = 1.0; // random value to start the iterative processus
277        while (Math.abs(lati1-lati) >= epsilon) {
278            lati = lati1;
279            double v1 = 1+e*Math.sin(lati);
280            double v2 = 1-e*Math.sin(lati);
281            lati1 = 2*Math.atan(Math.pow(v1/v2, e/2)*Math.exp(latIso))-Math.PI/2;
282        }
283        return lati1;
284    }
285
286    /**
287     * convert cartesian coordinates to ellipsoidal coordinates
288     *
289     * @param xyz the coordinates in meters (X, Y, Z)
290     * @return The corresponding latitude and longitude in degrees
291     */
292    public LatLon cart2LatLon(double[] xyz) {
293        return cart2LatLon(xyz, 1e-11);
294    }
295
296    public LatLon cart2LatLon(double[] xyz, double epsilon) {
297        double norm = Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
298        double lg = 2.0 * Math.atan(xyz[1] / (xyz[0] + norm));
299        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])))));
300        double delta = 1.0;
301        while (delta > epsilon) {
302            double s2 = Math.sin(lt);
303            s2 *= s2;
304            double l = Math.atan((xyz[2] / norm)
305                    / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2)))));
306            delta = Math.abs(l - lt);
307            lt = l;
308        }
309        return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg));
310    }
311
312    /**
313     * convert ellipsoidal coordinates to cartesian coordinates
314     *
315     * @param coord The Latitude and longitude in degrees
316     * @return the corresponding (X, Y Z) cartesian coordinates in meters.
317     */
318    public double[] latLon2Cart(LatLon coord) {
319        double phi = Math.toRadians(coord.lat());
320        double lambda = Math.toRadians(coord.lon());
321
322        double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2));
323        double[] xyz = new double[3];
324        xyz[0] = Rn * Math.cos(phi) * Math.cos(lambda);
325        xyz[1] = Rn * Math.cos(phi) * Math.sin(lambda);
326        xyz[2] = Rn * (1 - e2) * Math.sin(phi);
327
328        return xyz;
329    }
330}