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