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

Last change on this file since 13423 was 13423, checked in by bastiK, 6 years ago

see #15880 - add Fischer ellipsoid
  • Property svn:eol-style set to native
File size: 11.7 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;
10import org.openstreetmap.josm.tools.Utils;
11
12/**
13 * Reference ellipsoids.
14 */
15public final class Ellipsoid {
16
17 /**
18 * Airy 1830
19 */
20 public static final Ellipsoid Airy = Ellipsoid.createAb(6377563.396, 6356256.910);
21
22 /**
23 * Modified Airy 1849
24 */
25 public static final Ellipsoid AiryMod = Ellipsoid.createAb(6377340.189, 6356034.446);
26
27 /**
28 * Australian National Spheroid (Australian Natl & S. Amer. 1969)
29 * same as GRS67 Modified
30 */
31 public static final Ellipsoid AustSA = Ellipsoid.createArf(6378160.0, 298.25);
32
33 /**
34 * Bessel 1841 ellipsoid
35 */
36 public static final Ellipsoid Bessel1841 = Ellipsoid.createArf(6377397.155, 299.1528128);
37
38 /**
39 * Bessel 1841 (Namibia)
40 */
41 public static final Ellipsoid BesselNamibia = Ellipsoid.createArf(6377483.865, 299.1528128);
42
43 /**
44 * Clarke 1866 ellipsoid
45 */
46 public static final Ellipsoid Clarke1866 = Ellipsoid.createAb(6378206.4, 6356583.8);
47
48 /**
49 * Clarke 1880 (modified)
50 */
51 public static final Ellipsoid Clarke1880 = Ellipsoid.createArf(6378249.145, 293.4663);
52
53 /**
54 * Clarke 1880 IGN (French national geographic institute)
55 */
56 public static final Ellipsoid ClarkeIGN = Ellipsoid.createAb(6378249.2, 6356515.0);
57
58 /**
59 * Everest (Sabah & Sarawak)
60 */
61 public static final Ellipsoid EverestSabahSarawak = Ellipsoid.createArf(6377298.556, 300.8017);
62
63 /**
64 * Fischer (Mercury Datum) 1960
65 */
66 public static final Ellipsoid Fischer = Ellipsoid.createArf(6378166., 298.3);
67
68 /**
69 * Modified Fischer 1960
70 */
71 public static final Ellipsoid FischerMod = Ellipsoid.createArf(6378155., 298.3);
72
73 /**
74 * GRS67 ellipsoid
75 */
76 public static final Ellipsoid GRS67 = Ellipsoid.createArf(6378160.0, 298.247167427);
77
78 /**
79 * GRS80 ellipsoid
80 */
81 public static final Ellipsoid GRS80 = Ellipsoid.createArf(6378137.0, 298.257222101);
82
83 /**
84 * Hayford's ellipsoid 1909 (ED50 system)
85 * Also known as International 1924
86 * Proj.4 code: intl
87 */
88 public static final Ellipsoid Hayford = Ellipsoid.createArf(6378388.0, 297.0);
89
90 /**
91 * Helmert 1906
92 */
93 public static final Ellipsoid Helmert = Ellipsoid.createArf(6378200.0, 298.3);
94
95 /**
96 * Krassowsky 1940 ellipsoid
97 */
98 public static final Ellipsoid Krassowsky = Ellipsoid.createArf(6378245.0, 298.3);
99
100 /**
101 * WGS66 ellipsoid
102 */
103 public static final Ellipsoid WGS66 = Ellipsoid.createArf(6378145.0, 298.25);
104
105 /**
106 * WGS72 ellipsoid
107 */
108 public static final Ellipsoid WGS72 = Ellipsoid.createArf(6378135.0, 298.26);
109
110 /**
111 * WGS84 ellipsoid
112 */
113 public static final Ellipsoid WGS84 = Ellipsoid.createArf(6378137.0, 298.257223563);
114
115 /**
116 * half long axis
117 */
118 public final double a;
119
120 /**
121 * half short axis
122 */
123 public final double b;
124
125 /**
126 * first eccentricity:
127 * sqrt(a*a - b*b) / a
128 */
129 public final double e;
130
131 /**
132 * first eccentricity squared:
133 * (a*a - b*b) / (a*a)
134 */
135 public final double e2;
136
137 /**
138 * square of the second eccentricity:
139 * (a*a - b*b) / (b*b)
140 */
141 public final double eb2;
142
143 /**
144 * if ellipsoid is spherical, i.e. the major and minor semiaxis are
145 * the same
146 */
147 public final boolean spherical;
148
149 /**
150 * private constructur - use one of the create_* methods
151 *
152 * @param a semimajor radius of the ellipsoid axis
153 * @param b semiminor radius of the ellipsoid axis
154 * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a)))
155 * @param e2 first eccentricity squared
156 * @param eb2 square of the second eccentricity
157 * @param sperical if the ellipsoid is sphere
158 */
159 private Ellipsoid(double a, double b, double e, double e2, double eb2, boolean sperical) {
160 this.a = a;
161 this.b = b;
162 this.e = e;
163 this.e2 = e2;
164 this.eb2 = eb2;
165 this.spherical = sperical;
166 }
167
168 /**
169 * create a new ellipsoid
170 *
171 * @param a semimajor radius of the ellipsoid axis (in meters)
172 * @param b semiminor radius of the ellipsoid axis (in meters)
173 * @return the new ellipsoid
174 */
175 public static Ellipsoid createAb(double a, double b) {
176 double e2 = (a*a - b*b) / (a*a);
177 double e = Math.sqrt(e2);
178 double eb2 = e2 / (1.0 - e2);
179 return new Ellipsoid(a, b, e, e2, eb2, a == b);
180 }
181
182 /**
183 * create a new ellipsoid
184 *
185 * @param a semimajor radius of the ellipsoid axis (in meters)
186 * @param es first eccentricity squared
187 * @return the new ellipsoid
188 */
189 public static Ellipsoid createAes(double a, double es) {
190 double b = a * Math.sqrt(1.0 - es);
191 double e = Math.sqrt(es);
192 double eb2 = es / (1.0 - es);
193 return new Ellipsoid(a, b, e, es, eb2, es == 0);
194 }
195
196 /**
197 * create a new ellipsoid
198 *
199 * @param a semimajor radius of the ellipsoid axis (in meters)
200 * @param f flattening ( = (a - b) / a)
201 * @return the new ellipsoid
202 */
203 public static Ellipsoid createAf(double a, double f) {
204 double b = a * (1.0 - f);
205 double e2 = f * (2 - f);
206 double e = Math.sqrt(e2);
207 double eb2 = e2 / (1.0 - e2);
208 return new Ellipsoid(a, b, e, e2, eb2, f == 0);
209 }
210
211 /**
212 * create a new ellipsoid
213 *
214 * @param a semimajor radius of the ellipsoid axis (in meters)
215 * @param rf inverse flattening
216 * @return the new ellipsoid
217 */
218 public static Ellipsoid createArf(double a, double rf) {
219 return createAf(a, 1.0 / rf);
220 }
221
222 @Override
223 public String toString() {
224 return "Ellipsoid{a="+a+", b="+b+'}';
225 }
226
227 /**
228 * Returns the <i>radius of curvature in the prime vertical</i>
229 * for this reference ellipsoid at the specified latitude.
230 *
231 * @param phi The local latitude (radians).
232 * @return The radius of curvature in the prime vertical (meters).
233 */
234 public double verticalRadiusOfCurvature(final double phi) {
235 return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
236 }
237
238 private static double sqr(final double x) {
239 return x * x;
240 }
241
242 /**
243 * Returns the meridional arc, the true meridional distance on the
244 * ellipsoid from the equator to the specified latitude, in meters.
245 *
246 * @param phi The local latitude (in radians).
247 * @return The meridional arc (in meters).
248 */
249 public double meridionalArc(final double phi) {
250 final double sin2Phi = Math.sin(2.0 * phi);
251 final double sin4Phi = Math.sin(4.0 * phi);
252 final double sin6Phi = Math.sin(6.0 * phi);
253 final double sin8Phi = Math.sin(8.0 * phi);
254 // TODO . calculate 'f'
255 //double f = 1.0 / 298.257222101; // GRS80
256 double f = 1.0 / 298.257223563; // WGS84
257 final double n = f / (2.0 - f);
258 final double n2 = n * n;
259 final double n3 = n2 * n;
260 final double n4 = n3 * n;
261 final double n5 = n4 * n;
262 final double n1n2 = n - n2;
263 final double n2n3 = n2 - n3;
264 final double n3n4 = n3 - n4;
265 final double n4n5 = n4 - n5;
266 final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
267 final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
268 final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
269 final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
270 final double ep = (315.0 / 512.0) * a * (n4n5);
271 return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
272 }
273
274 /**
275 * Returns the <i>radius of curvature in the meridian</i>
276 * for this reference ellipsoid at the specified latitude.
277 *
278 * @param phi The local latitude (in radians).
279 * @return The radius of curvature in the meridian (in meters).
280 */
281 public double meridionalRadiusOfCurvature(final double phi) {
282 return verticalRadiusOfCurvature(phi)
283 / (1.0 + eb2 * sqr(Math.cos(phi)));
284 }
285
286 /**
287 * Returns isometric latitude of phi on given first eccentricity (e)
288 * @param phi The local latitude (radians).
289 * @param e first eccentricity
290 * @return isometric latitude of phi on first eccentricity (e)
291 */
292 public double latitudeIsometric(double phi, double e) {
293 double v1 = 1-e*Math.sin(phi);
294 double v2 = 1+e*Math.sin(phi);
295 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
296 }
297
298 /**
299 * Returns isometric latitude of phi on first eccentricity (e)
300 * @param phi The local latitude (radians).
301 * @return isometric latitude of phi on first eccentricity (e)
302 */
303 public double latitudeIsometric(double phi) {
304 double v1 = 1-e*Math.sin(phi);
305 double v2 = 1+e*Math.sin(phi);
306 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
307 }
308
309 /**
310 * Returns geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision
311 * @param latIso isometric latitude
312 * @param e first eccentricity
313 * @param epsilon epsilon precision
314 * @return geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision
315 */
316 public double latitude(double latIso, double e, double epsilon) {
317 double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
318 double lati = lat0;
319 double lati1 = 1.0; // random value to start the iterative processus
320 while (Math.abs(lati1-lati) >= epsilon) {
321 lati = lati1;
322 double v1 = 1+e*Math.sin(lati);
323 double v2 = 1-e*Math.sin(lati);
324 lati1 = 2*Math.atan(Math.pow(v1/v2, e/2)*Math.exp(latIso))-Math.PI/2;
325 }
326 return lati1;
327 }
328
329 /**
330 * convert cartesian coordinates to ellipsoidal coordinates
331 *
332 * @param xyz the coordinates in meters (X, Y, Z)
333 * @return The corresponding latitude and longitude in degrees
334 */
335 public LatLon cart2LatLon(double... xyz) {
336 return cart2LatLon(xyz, 1e-11);
337 }
338
339 public LatLon cart2LatLon(double[] xyz, double epsilon) {
340 double norm = Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
341 double lg = 2.0 * Math.atan(xyz[1] / (xyz[0] + norm));
342 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])))));
343 double delta = 1.0;
344 while (delta > epsilon) {
345 double s2 = Math.sin(lt);
346 s2 *= s2;
347 double l = Math.atan((xyz[2] / norm)
348 / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2)))));
349 delta = Math.abs(l - lt);
350 lt = l;
351 }
352 return new LatLon(Utils.toDegrees(lt), Utils.toDegrees(lg));
353 }
354
355 /**
356 * convert ellipsoidal coordinates to cartesian coordinates
357 *
358 * @param coord The Latitude and longitude in degrees
359 * @return the corresponding (X, Y Z) cartesian coordinates in meters.
360 */
361 public double[] latLon2Cart(LatLon coord) {
362 double phi = Utils.toRadians(coord.lat());
363 double lambda = Utils.toRadians(coord.lon());
364
365 double rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2));
366 double[] xyz = new double[3];
367 xyz[0] = rn * Math.cos(phi) * Math.cos(lambda);
368 xyz[1] = rn * Math.cos(phi) * Math.sin(lambda);
369 xyz[2] = rn * (1 - e2) * Math.sin(phi);
370
371 return xyz;
372 }
373}
Note: See TracBrowser for help on using the repository browser.