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

Last change on this file since 11266 was 10748, checked in by Don-vip, 8 years ago

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