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