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