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

Last change on this file since 7936 was 7402, checked in by Don-vip, 10 years ago

fix some Sonar issues
  • Property svn:eol-style set to native
File size: 9.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 */
7
8package org.openstreetmap.josm.data.projection;
9
10import org.openstreetmap.josm.data.coor.LatLon;
11
12/**
13 * Reference ellipsoids.
14 */
15public final class Ellipsoid {
16
17 /**
18 * Clarke 1880 IGN (French national geographic institute)
19 */
20 public static final Ellipsoid clarkeIGN = Ellipsoid.create_a_b(6378249.2, 6356515.0);
21
22 /**
23 * Hayford's ellipsoid 1909 (ED50 system)<br>
24 * Proj.4 code: intl
25 */
26 public static final Ellipsoid hayford = Ellipsoid.create_a_rf(6378388.0, 297.0);
27
28 /**
29 * GRS67 ellipsoid
30 */
31 public static final Ellipsoid GRS67 = Ellipsoid.create_a_rf(6378160.0, 298.247167472);
32
33 /**
34 * GRS80 ellipsoid
35 */
36 public static final Ellipsoid GRS80 = Ellipsoid.create_a_rf(6378137.0, 298.257222101);
37
38 /**
39 * WGS84 ellipsoid
40 */
41 public static final Ellipsoid WGS84 = Ellipsoid.create_a_rf(6378137.0, 298.257223563);
42
43 /**
44 * Bessel 1841 ellipsoid
45 */
46 public static final Ellipsoid Bessel1841 = Ellipsoid.create_a_rf(6377397.155, 299.1528128);
47
48 /**
49 * half long axis
50 */
51 public final double a;
52
53 /**
54 * half short axis
55 */
56 public final double b;
57
58 /**
59 * first eccentricity
60 */
61 public final double e;
62
63 /**
64 * first eccentricity squared
65 */
66 public final double e2;
67
68 /**
69 * square of the second eccentricity
70 */
71 public final double eb2;
72
73 /**
74 * private constructur - use one of the create_* methods
75 *
76 * @param a semimajor radius of the ellipsoid axis
77 * @param b semiminor radius of the ellipsoid axis
78 * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a)))
79 * @param e2 first eccentricity squared
80 * @param eb2 square of the second eccentricity
81 */
82 private Ellipsoid(double a, double b, double e, double e2, double eb2) {
83 this.a = a;
84 this.b = b;
85 this.e = e;
86 this.e2 = e2;
87 this.eb2 = eb2;
88 }
89
90 /**
91 * create a new ellipsoid
92 *
93 * @param a semimajor radius of the ellipsoid axis (in meters)
94 * @param b semiminor radius of the ellipsoid axis (in meters)
95 * @return the new ellipsoid
96 */
97 public static Ellipsoid create_a_b(double a, double b) {
98 double e2 = (a*a - b*b) / (a*a);
99 double e = Math.sqrt(e2);
100 double eb2 = e2 / (1.0 - e2);
101 return new Ellipsoid(a, b, e, e2, eb2);
102 }
103
104 /**
105 * create a new ellipsoid
106 *
107 * @param a semimajor radius of the ellipsoid axis (in meters)
108 * @param es first eccentricity squared
109 * @return the new ellipsoid
110 */
111 public static Ellipsoid create_a_es(double a, double es) {
112 double b = a * Math.sqrt(1.0 - es);
113 double e = Math.sqrt(es);
114 double eb2 = es / (1.0 - es);
115 return new Ellipsoid(a, b, e, es, eb2);
116 }
117
118 /**
119 * create a new ellipsoid
120 *
121 * @param a semimajor radius of the ellipsoid axis (in meters)
122 * @param f flattening ( = (a - b) / a)
123 * @return the new ellipsoid
124 */
125 public static Ellipsoid create_a_f(double a, double f) {
126 double b = a * (1.0 - f);
127 double e2 = f * (2 - f);
128 double e = Math.sqrt(e2);
129 double eb2 = e2 / (1.0 - e2);
130 return new Ellipsoid(a, b, e, e2, eb2);
131 }
132
133 /**
134 * create a new ellipsoid
135 *
136 * @param a semimajor radius of the ellipsoid axis (in meters)
137 * @param rf inverse flattening
138 * @return the new ellipsoid
139 */
140 public static Ellipsoid create_a_rf(double a, double rf) {
141 return create_a_f(a, 1.0 / rf);
142 }
143
144 @Override
145 public String toString() {
146 return "Ellipsoid{a="+a+", b="+b+"}";
147 }
148
149 /**
150 * Returns the <i>radius of curvature in the prime vertical</i>
151 * for this reference ellipsoid at the specified latitude.
152 *
153 * @param phi The local latitude (radians).
154 * @return The radius of curvature in the prime vertical (meters).
155 */
156 public double verticalRadiusOfCurvature(final double phi) {
157 return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
158 }
159
160 private static double sqr(final double x) {
161 return x * x;
162 }
163
164 /**
165 * Returns the meridional arc, the true meridional distance on the
166 * ellipsoid from the equator to the specified latitude, in meters.
167 *
168 * @param phi The local latitude (in radians).
169 * @return The meridional arc (in meters).
170 */
171 public double meridionalArc(final double phi) {
172 final double sin2Phi = Math.sin(2.0 * phi);
173 final double sin4Phi = Math.sin(4.0 * phi);
174 final double sin6Phi = Math.sin(6.0 * phi);
175 final double sin8Phi = Math.sin(8.0 * phi);
176 // TODO . calculate 'f'
177 //double f = 1.0 / 298.257222101; // GRS80
178 double f = 1.0 / 298.257223563; // WGS84
179 final double n = f / (2.0 - f);
180 final double n2 = n * n;
181 final double n3 = n2 * n;
182 final double n4 = n3 * n;
183 final double n5 = n4 * n;
184 final double n1n2 = n - n2;
185 final double n2n3 = n2 - n3;
186 final double n3n4 = n3 - n4;
187 final double n4n5 = n4 - n5;
188 final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
189 final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
190 final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
191 final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
192 final double ep = (315.0 / 512.0) * a * (n4n5);
193 return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
194 }
195
196 /**
197 * Returns the <i>radius of curvature in the meridian</i>
198 * for this reference ellipsoid at the specified latitude.
199 *
200 * @param phi The local latitude (in radians).
201 * @return The radius of curvature in the meridian (in meters).
202 */
203 public double meridionalRadiusOfCurvature(final double phi) {
204 return verticalRadiusOfCurvature(phi)
205 / (1.0 + eb2 * sqr(Math.cos(phi)));
206 }
207
208 /**
209 * Returns isometric latitude of phi on given first eccentricity (e)
210 * @param phi The local latitude (radians).
211 * @return isometric latitude of phi on first eccentricity (e)
212 */
213 public double latitudeIsometric(double phi, double e) {
214 double v1 = 1-e*Math.sin(phi);
215 double v2 = 1+e*Math.sin(phi);
216 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
217 }
218
219 /**
220 * Returns isometric latitude of phi on first eccentricity (e)
221 * @param phi The local latitude (radians).
222 * @return isometric latitude of phi on first eccentricity (e)
223 */
224 public double latitudeIsometric(double phi) {
225 double v1 = 1-e*Math.sin(phi);
226 double v2 = 1+e*Math.sin(phi);
227 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
228 }
229
230 /**
231 * Returns geographic latitude of isometric latitude of first eccentricity (e)
232 * and epsilon precision
233 * @return geographic latitude of isometric latitude of first eccentricity (e)
234 * and epsilon precision
235 */
236 public double latitude(double latIso, double e, double epsilon) {
237 double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
238 double lati = lat0;
239 double lati1 = 1.0; // random value to start the iterative processus
240 while(Math.abs(lati1-lati)>=epsilon) {
241 lati = lati1;
242 double v1 = 1+e*Math.sin(lati);
243 double v2 = 1-e*Math.sin(lati);
244 lati1 = 2*Math.atan(Math.pow(v1/v2,e/2)*Math.exp(latIso))-Math.PI/2;
245 }
246 return lati1;
247 }
248
249 /**
250 * convert cartesian coordinates to ellipsoidal coordinates
251 *
252 * @param xyz the coordinates in meters (X, Y, Z)
253 * @return The corresponding latitude and longitude in degrees
254 */
255 public LatLon cart2LatLon(double[] xyz) {
256 return cart2LatLon(xyz, 1e-11);
257 }
258
259 public LatLon cart2LatLon(double[] xyz, double epsilon) {
260 double norm = Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
261 double lg = 2.0 * Math.atan(xyz[1] / (xyz[0] + norm));
262 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])))));
263 double delta = 1.0;
264 while (delta > epsilon) {
265 double s2 = Math.sin(lt);
266 s2 *= s2;
267 double l = Math.atan((xyz[2] / norm)
268 / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2)))));
269 delta = Math.abs(l - lt);
270 lt = l;
271 }
272 return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg));
273 }
274
275 /**
276 * convert ellipsoidal coordinates to cartesian coordinates
277 *
278 * @param coord The Latitude and longitude in degrees
279 * @return the corresponding (X, Y Z) cartesian coordinates in meters.
280 */
281 public double[] latLon2Cart(LatLon coord) {
282 double phi = Math.toRadians(coord.lat());
283 double lambda = Math.toRadians(coord.lon());
284
285 double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2));
286 double[] xyz = new double[3];
287 xyz[0] = Rn * Math.cos(phi) * Math.cos(lambda);
288 xyz[1] = Rn * Math.cos(phi) * Math.sin(lambda);
289 xyz[2] = Rn * (1 - e2) * Math.sin(phi);
290
291 return xyz;
292 }
293}
Note: See TracBrowser for help on using the repository browser.