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

Last change on this file since 6463 was 6463, checked in by bastiK, 10 years ago

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