001/* 002 * Import from fr.geo.convert package, a geographic coordinates converter. 003 * (https://www.i3s.unice.fr/~johan/gps/) 004 * License: GPL. For details, see LICENSE file. 005 * Copyright (C) 2002 Johan Montagnat (johan@creatis.insa-lyon.fr) 006 */ 007package org.openstreetmap.josm.data.projection; 008 009import org.openstreetmap.josm.data.coor.LatLon; 010 011/** 012 * Reference ellipsoids. 013 */ 014public final class Ellipsoid { 015 016 /** 017 * Airy 1830 018 */ 019 public static final Ellipsoid Airy = Ellipsoid.create_a_b(6377563.396, 6356256.910); 020 021 /** 022 * Modified Airy 1849 023 */ 024 public static final Ellipsoid AiryMod = Ellipsoid.create_a_b(6377340.189, 6356034.446); 025 026 /** 027 * Australian National Spheroid (Australian Natl & S. Amer. 1969) 028 * same as GRS67 Modified 029 */ 030 public static final Ellipsoid AustSA = Ellipsoid.create_a_rf(6378160.0, 298.25); 031 032 /** 033 * Bessel 1841 ellipsoid 034 */ 035 public static final Ellipsoid Bessel1841 = Ellipsoid.create_a_rf(6377397.155, 299.1528128); 036 037 /** 038 * Bessel 1841 (Namibia) 039 */ 040 public static final Ellipsoid BesselNamibia = Ellipsoid.create_a_rf(6377483.865, 299.1528128); 041 042 /** 043 * Clarke 1866 ellipsoid 044 */ 045 public static final Ellipsoid Clarke1866 = Ellipsoid.create_a_b(6378206.4, 6356583.8); 046 047 /** 048 * Clarke 1880 (modified) 049 */ 050 public static final Ellipsoid Clarke1880 = Ellipsoid.create_a_rf(6378249.145, 293.4663); 051 052 /** 053 * Clarke 1880 IGN (French national geographic institute) 054 */ 055 public static final Ellipsoid ClarkeIGN = Ellipsoid.create_a_b(6378249.2, 6356515.0); 056 057 /** 058 * Everest (Sabah & Sarawak) 059 */ 060 public static final Ellipsoid EverestSabahSarawak = Ellipsoid.create_a_rf(6377298.556, 300.8017); 061 062 /** 063 * GRS67 ellipsoid 064 */ 065 public static final Ellipsoid GRS67 = Ellipsoid.create_a_rf(6378160.0, 298.247167427); 066 067 /** 068 * GRS80 ellipsoid 069 */ 070 public static final Ellipsoid GRS80 = Ellipsoid.create_a_rf(6378137.0, 298.257222101); 071 072 /** 073 * Hayford's ellipsoid 1909 (ED50 system) 074 * Also known as International 1924 075 * Proj.4 code: intl 076 */ 077 public static final Ellipsoid Hayford = Ellipsoid.create_a_rf(6378388.0, 297.0); 078 079 /** 080 * Helmert 1906 081 */ 082 public static final Ellipsoid Helmert = Ellipsoid.create_a_rf(6378200.0, 298.3); 083 084 /** 085 * Krassowsky 1940 ellipsoid 086 */ 087 public static final Ellipsoid Krassowsky = Ellipsoid.create_a_rf(6378245.0, 298.3); 088 089 /** 090 * WGS66 ellipsoid 091 */ 092 public static final Ellipsoid WGS66 = Ellipsoid.create_a_rf(6378145.0, 298.25); 093 094 /** 095 * WGS72 ellipsoid 096 */ 097 public static final Ellipsoid WGS72 = Ellipsoid.create_a_rf(6378135.0, 298.26); 098 099 /** 100 * WGS84 ellipsoid 101 */ 102 public static final Ellipsoid WGS84 = Ellipsoid.create_a_rf(6378137.0, 298.257223563); 103 104 105 /** 106 * half long axis 107 */ 108 public final double a; 109 110 /** 111 * half short axis 112 */ 113 public final double b; 114 115 /** 116 * first eccentricity: 117 * sqrt(a*a - b*b) / a 118 */ 119 public final double e; 120 121 /** 122 * first eccentricity squared: 123 * (a*a - b*b) / (a*a) 124 */ 125 public final double e2; 126 127 /** 128 * square of the second eccentricity: 129 * (a*a - b*b) / (b*b) 130 */ 131 public final double eb2; 132 133 /** 134 * if ellipsoid is spherical, i.e. the major and minor semiaxis are 135 * the same 136 */ 137 public final boolean spherical; 138 139 /** 140 * private constructur - use one of the create_* methods 141 * 142 * @param a semimajor radius of the ellipsoid axis 143 * @param b semiminor radius of the ellipsoid axis 144 * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a))) 145 * @param e2 first eccentricity squared 146 * @param eb2 square of the second eccentricity 147 * @param sperical if the ellipsoid is sphere 148 */ 149 private Ellipsoid(double a, double b, double e, double e2, double eb2, boolean sperical) { 150 this.a = a; 151 this.b = b; 152 this.e = e; 153 this.e2 = e2; 154 this.eb2 = eb2; 155 this.spherical = sperical; 156 } 157 158 /** 159 * create a new ellipsoid 160 * 161 * @param a semimajor radius of the ellipsoid axis (in meters) 162 * @param b semiminor radius of the ellipsoid axis (in meters) 163 * @return the new ellipsoid 164 */ 165 public static Ellipsoid create_a_b(double a, double b) { 166 double e2 = (a*a - b*b) / (a*a); 167 double e = Math.sqrt(e2); 168 double eb2 = e2 / (1.0 - e2); 169 return new Ellipsoid(a, b, e, e2, eb2, a == b); 170 } 171 172 /** 173 * create a new ellipsoid 174 * 175 * @param a semimajor radius of the ellipsoid axis (in meters) 176 * @param es first eccentricity squared 177 * @return the new ellipsoid 178 */ 179 public static Ellipsoid create_a_es(double a, double es) { 180 double b = a * Math.sqrt(1.0 - es); 181 double e = Math.sqrt(es); 182 double eb2 = es / (1.0 - es); 183 return new Ellipsoid(a, b, e, es, eb2, es == 0); 184 } 185 186 /** 187 * create a new ellipsoid 188 * 189 * @param a semimajor radius of the ellipsoid axis (in meters) 190 * @param f flattening ( = (a - b) / a) 191 * @return the new ellipsoid 192 */ 193 public static Ellipsoid create_a_f(double a, double f) { 194 double b = a * (1.0 - f); 195 double e2 = f * (2 - f); 196 double e = Math.sqrt(e2); 197 double eb2 = e2 / (1.0 - e2); 198 return new Ellipsoid(a, b, e, e2, eb2, f == 0); 199 } 200 201 /** 202 * create a new ellipsoid 203 * 204 * @param a semimajor radius of the ellipsoid axis (in meters) 205 * @param rf inverse flattening 206 * @return the new ellipsoid 207 */ 208 public static Ellipsoid create_a_rf(double a, double rf) { 209 return create_a_f(a, 1.0 / rf); 210 } 211 212 @Override 213 public String toString() { 214 return "Ellipsoid{a="+a+", b="+b+'}'; 215 } 216 217 /** 218 * Returns the <i>radius of curvature in the prime vertical</i> 219 * for this reference ellipsoid at the specified latitude. 220 * 221 * @param phi The local latitude (radians). 222 * @return The radius of curvature in the prime vertical (meters). 223 */ 224 public double verticalRadiusOfCurvature(final double phi) { 225 return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi)))); 226 } 227 228 private static double sqr(final double x) { 229 return x * x; 230 } 231 232 /** 233 * Returns the meridional arc, the true meridional distance on the 234 * ellipsoid from the equator to the specified latitude, in meters. 235 * 236 * @param phi The local latitude (in radians). 237 * @return The meridional arc (in meters). 238 */ 239 public double meridionalArc(final double phi) { 240 final double sin2Phi = Math.sin(2.0 * phi); 241 final double sin4Phi = Math.sin(4.0 * phi); 242 final double sin6Phi = Math.sin(6.0 * phi); 243 final double sin8Phi = Math.sin(8.0 * phi); 244 // TODO . calculate 'f' 245 //double f = 1.0 / 298.257222101; // GRS80 246 double f = 1.0 / 298.257223563; // WGS84 247 final double n = f / (2.0 - f); 248 final double n2 = n * n; 249 final double n3 = n2 * n; 250 final double n4 = n3 * n; 251 final double n5 = n4 * n; 252 final double n1n2 = n - n2; 253 final double n2n3 = n2 - n3; 254 final double n3n4 = n3 - n4; 255 final double n4n5 = n4 - n5; 256 final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5)); 257 final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5); 258 final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5)); 259 final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5); 260 final double ep = (315.0 / 512.0) * a * (n4n5); 261 return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi; 262 } 263 264 /** 265 * Returns the <i>radius of curvature in the meridian</i> 266 * for this reference ellipsoid at the specified latitude. 267 * 268 * @param phi The local latitude (in radians). 269 * @return The radius of curvature in the meridian (in meters). 270 */ 271 public double meridionalRadiusOfCurvature(final double phi) { 272 return verticalRadiusOfCurvature(phi) 273 / (1.0 + eb2 * sqr(Math.cos(phi))); 274 } 275 276 /** 277 * Returns isometric latitude of phi on given first eccentricity (e) 278 * @param phi The local latitude (radians). 279 * @param e first eccentricity 280 * @return isometric latitude of phi on first eccentricity (e) 281 */ 282 public double latitudeIsometric(double phi, double e) { 283 double v1 = 1-e*Math.sin(phi); 284 double v2 = 1+e*Math.sin(phi); 285 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2)); 286 } 287 288 /** 289 * Returns isometric latitude of phi on first eccentricity (e) 290 * @param phi The local latitude (radians). 291 * @return isometric latitude of phi on first eccentricity (e) 292 */ 293 public double latitudeIsometric(double phi) { 294 double v1 = 1-e*Math.sin(phi); 295 double v2 = 1+e*Math.sin(phi); 296 return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2)); 297 } 298 299 /** 300 * Returns geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision 301 * @param latIso isometric latitude 302 * @param e first eccentricity 303 * @param epsilon epsilon precision 304 * @return geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision 305 */ 306 public double latitude(double latIso, double e, double epsilon) { 307 double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2; 308 double lati = lat0; 309 double lati1 = 1.0; // random value to start the iterative processus 310 while (Math.abs(lati1-lati) >= epsilon) { 311 lati = lati1; 312 double v1 = 1+e*Math.sin(lati); 313 double v2 = 1-e*Math.sin(lati); 314 lati1 = 2*Math.atan(Math.pow(v1/v2, e/2)*Math.exp(latIso))-Math.PI/2; 315 } 316 return lati1; 317 } 318 319 /** 320 * convert cartesian coordinates to ellipsoidal coordinates 321 * 322 * @param xyz the coordinates in meters (X, Y, Z) 323 * @return The corresponding latitude and longitude in degrees 324 */ 325 public LatLon cart2LatLon(double[] xyz) { 326 return cart2LatLon(xyz, 1e-11); 327 } 328 329 public LatLon cart2LatLon(double[] xyz, double epsilon) { 330 double norm = Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]); 331 double lg = 2.0 * Math.atan(xyz[1] / (xyz[0] + norm)); 332 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]))))); 333 double delta = 1.0; 334 while (delta > epsilon) { 335 double s2 = Math.sin(lt); 336 s2 *= s2; 337 double l = Math.atan((xyz[2] / norm) 338 / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2))))); 339 delta = Math.abs(l - lt); 340 lt = l; 341 } 342 return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg)); 343 } 344 345 /** 346 * convert ellipsoidal coordinates to cartesian coordinates 347 * 348 * @param coord The Latitude and longitude in degrees 349 * @return the corresponding (X, Y Z) cartesian coordinates in meters. 350 */ 351 public double[] latLon2Cart(LatLon coord) { 352 double phi = Math.toRadians(coord.lat()); 353 double lambda = Math.toRadians(coord.lon()); 354 355 double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2)); 356 double[] xyz = new double[3]; 357 xyz[0] = Rn * Math.cos(phi) * Math.cos(lambda); 358 xyz[1] = Rn * Math.cos(phi) * Math.sin(lambda); 359 xyz[2] = Rn * (1 - e2) * Math.sin(phi); 360 361 return xyz; 362 } 363}