source: josm/trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java@ 8470

Last change on this file since 8470 was 8451, checked in by Don-vip, 9 years ago

fix #11512 - Add most commonly used ellipsoids (patch by BathoryPeter)
  • Property svn:eol-style set to native
File size: 10.4 KB
Line 
1/*
2 * Import from fr.geo.convert package, a geographic coordinates converter.
3 * (https://www.i3s.unice.fr/~johan/gps/)
4 * License: GPL. For details, see LICENSE file.
5 * Copyright (C) 2002 Johan Montagnat (johan@creatis.insa-lyon.fr)
6 */
7package org.openstreetmap.josm.data.projection;
8
9import org.openstreetmap.josm.data.coor.LatLon;
10
11/**
12 * Reference ellipsoids.
13 */
14public final class Ellipsoid {
15
16 /**
17 * Airy 1830
18 */
19 public static final Ellipsoid Airy = Ellipsoid.create_a_b(6377563.396, 6356256.910);
20
21 /**
22 * Modified Airy 1849
23 */
24 public static final Ellipsoid AiryMod = Ellipsoid.create_a_b(6377340.189, 6356034.446);
25
26 /**
27 * Australian National Spheroid (Australian Natl & S. Amer. 1969)
28 * same as GRS67 Modified
29 */
30 public static final Ellipsoid AustSA = Ellipsoid.create_a_rf(6378160.0, 298.25);
31
32 /**
33 * Bessel 1841 ellipsoid
34 */
35 public static final Ellipsoid Bessel1841 = Ellipsoid.create_a_rf(6377397.155, 299.1528128);
36
37 /**
38 * Clarke 1866 ellipsoid
39 */
40 public static final Ellipsoid Clarke1866 = Ellipsoid.create_a_b(6378206.4, 6356583.8);
41
42 /**
43 * Clarke 1880 IGN (French national geographic institute)
44 */
45 public static final Ellipsoid ClarkeIGN = Ellipsoid.create_a_b(6378249.2, 6356515.0);
46
47 /**
48 * GRS67 ellipsoid
49 */
50 public static final Ellipsoid GRS67 = Ellipsoid.create_a_rf(6378160.0, 298.247167427);
51
52 /**
53 * GRS80 ellipsoid
54 */
55 public static final Ellipsoid GRS80 = Ellipsoid.create_a_rf(6378137.0, 298.257222101);
56
57 /**
58 * Hayford's ellipsoid 1909 (ED50 system)
59 * Also known as International 1924
60 * Proj.4 code: intl
61 */
62 public static final Ellipsoid Hayford = Ellipsoid.create_a_rf(6378388.0, 297.0);
63
64 /**
65 * Helmert 1906
66 */
67 public static final Ellipsoid Helmert = Ellipsoid.create_a_rf(6378200.0, 298.3);
68
69 /**
70 * Krassowsky 1940 ellipsoid
71 */
72 public static final Ellipsoid Krassowsky = Ellipsoid.create_a_rf(6378245.0, 298.3);
73
74 /**
75 * WGS72 ellipsoid
76 */
77 public static final Ellipsoid WGS72 = Ellipsoid.create_a_rf(6378135.0, 298.26);
78
79 /**
80 * WGS84 ellipsoid
81 */
82 public static final Ellipsoid WGS84 = Ellipsoid.create_a_rf(6378137.0, 298.257223563);
83
84
85 /**
86 * half long axis
87 */
88 public final double a;
89
90 /**
91 * half short axis
92 */
93 public final double b;
94
95 /**
96 * first eccentricity
97 */
98 public final double e;
99
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}
Note: See TracBrowser for help on using the repository browser.