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

Last change on this file since 4275 was 3530, checked in by stoecker, 14 years ago

fix array preferences
  • Property svn:eol-style set to native
File size: 7.2 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 class Ellipsoid {
16 /**
17 * Clarke's ellipsoid (NTF system)
18 */
19 public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
20 /**
21 * Hayford's ellipsoid (ED50 system)
22 */
23 public static final Ellipsoid hayford =
24 new Ellipsoid(6378388.0, 6356911.9461);
25 /**
26 * GRS80 ellipsoid
27 */
28 public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.3141);
29
30 /**
31 * WGS84 ellipsoid
32 */
33 public static final Ellipsoid WGS84 = new Ellipsoid(6378137.0, 6356752.3142);
34
35 /**
36 * Bessel 1841 ellipsoid
37 */
38 public static final Ellipsoid Bessel1841 = new Ellipsoid(6377397.155, 6356078.962822);
39
40 /**
41 * half long axis
42 */
43 public final double a;
44 /**
45 * half short axis
46 */
47 public final double b;
48 /**
49 * first eccentricity
50 */
51 public final double e;
52 /**
53 * first eccentricity squared
54 */
55 public final double e2;
56
57 /**
58 * square of the second eccentricity
59 */
60 public final double eb2;
61
62 /**
63 * create a new ellipsoid and precompute its parameters
64 *
65 * @param a ellipsoid long axis (in meters)
66 * @param b ellipsoid short axis (in meters)
67 */
68 public Ellipsoid(double a, double b) {
69 this.a = a;
70 this.b = b;
71 e2 = (a*a - b*b) / (a*a);
72 e = Math.sqrt(e2);
73 eb2 = e2 / (1.0 - e2);
74 }
75
76 /**
77 * Returns the <i>radius of curvature in the prime vertical</i>
78 * for this reference ellipsoid at the specified latitude.
79 *
80 * @param phi The local latitude (radians).
81 * @return The radius of curvature in the prime vertical (meters).
82 */
83 public double verticalRadiusOfCurvature(final double phi) {
84 return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
85 }
86
87 private static double sqr(final double x) {
88 return x * x;
89 }
90
91 /**
92 * Returns the meridional arc, the true meridional distance on the
93 * ellipsoid from the equator to the specified latitude, in meters.
94 *
95 * @param phi The local latitude (in radians).
96 * @return The meridional arc (in meters).
97 */
98 public double meridionalArc(final double phi) {
99 final double sin2Phi = Math.sin(2.0 * phi);
100 final double sin4Phi = Math.sin(4.0 * phi);
101 final double sin6Phi = Math.sin(6.0 * phi);
102 final double sin8Phi = Math.sin(8.0 * phi);
103 // TODO . calculate 'f'
104 //double f = 1.0 / 298.257222101; // GRS80
105 double f = 1.0 / 298.257223563; // WGS84
106 final double n = f / (2.0 - f);
107 final double n2 = n * n;
108 final double n3 = n2 * n;
109 final double n4 = n3 * n;
110 final double n5 = n4 * n;
111 final double n1n2 = n - n2;
112 final double n2n3 = n2 - n3;
113 final double n3n4 = n3 - n4;
114 final double n4n5 = n4 - n5;
115 final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
116 final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
117 final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
118 final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
119 final double ep = (315.0 / 512.0) * a * (n4n5);
120 return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
121 }
122
123 /**
124 * Returns the <i>radius of curvature in the meridian<i>
125 * for this reference ellipsoid at the specified latitude.
126 *
127 * @param phi The local latitude (in radians).
128 * @return The radius of curvature in the meridian (in meters).
129 */
130 public double meridionalRadiusOfCurvature(final double phi) {
131 return verticalRadiusOfCurvature(phi)
132 / (1.0 + eb2 * sqr(Math.cos(phi)));
133 }
134
135 /**
136 * Returns isometric latitude of phi on given first eccentricity (e)
137 * @param phi The local latitude (radians).
138 * @return isometric latitude of phi on first eccentricity (e)
139 */
140 public double latitudeIsometric(double phi, double e) {
141 double v1 = 1-e*Math.sin(phi);
142 double v2 = 1+e*Math.sin(phi);
143 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
144 }
145
146 /**
147 * Returns isometric latitude of phi on first eccentricity (e)
148 * @param phi The local latitude (radians).
149 * @return isometric latitude of phi on first eccentricity (e)
150 */
151 public double latitudeIsometric(double phi) {
152 double v1 = 1-e*Math.sin(phi);
153 double v2 = 1+e*Math.sin(phi);
154 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
155 }
156
157 /*
158 * Returns geographic latitude of isometric latitude of first eccentricity (e)
159 * and epsilon precision
160 */
161 public double latitude(double latIso, double e, double epsilon) {
162 double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
163 double lati = lat0;
164 double lati1 = 1.0; // random value to start the iterative processus
165 while(Math.abs(lati1-lati)>=epsilon) {
166 lati = lati1;
167 double v1 = 1+e*Math.sin(lati);
168 double v2 = 1-e*Math.sin(lati);
169 lati1 = 2*Math.atan(Math.pow(v1/v2,e/2)*Math.exp(latIso))-Math.PI/2;
170 }
171 return lati1;
172 }
173
174 /**
175 * convert cartesian coordinates to ellipsoidal coordinates
176 *
177 * @param XYZ the coordinates in meters (X, Y, Z)
178 * @return The corresponding latitude and longitude in degrees
179 */
180 public LatLon cart2LatLon(double[] XYZ) {
181 return cart2LatLon(XYZ, 1e-11);
182 }
183 public LatLon cart2LatLon(double[] XYZ, double epsilon) {
184 double norm = Math.sqrt(XYZ[0] * XYZ[0] + XYZ[1] * XYZ[1]);
185 double lg = 2.0 * Math.atan(XYZ[1] / (XYZ[0] + norm));
186 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])))));
187 double delta = 1.0;
188 while (delta > epsilon) {
189 double s2 = Math.sin(lt);
190 s2 *= s2;
191 double l = Math.atan((XYZ[2] / norm)
192 / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2)))));
193 delta = Math.abs(l - lt);
194 lt = l;
195 }
196 return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg));
197 }
198
199 /**
200 * convert ellipsoidal coordinates to cartesian coordinates
201 *
202 * @param coord The Latitude and longitude in degrees
203 * @return the corresponding (X, Y Z) cartesian coordinates in meters.
204 */
205 public double[] latLon2Cart(LatLon coord) {
206 double phi = Math.toRadians(coord.lat());
207 double lambda = Math.toRadians(coord.lon());
208
209 double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2));
210 double[] XYZ = new double[3];
211 XYZ[0] = Rn * Math.cos(phi) * Math.cos(lambda);
212 XYZ[1] = Rn * Math.cos(phi) * Math.sin(lambda);
213 XYZ[2] = Rn * (1 - e2) * Math.sin(phi);
214
215 return XYZ;
216 }
217}
Note: See TracBrowser for help on using the repository browser.