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 }