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