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 */ 007 008package org.openstreetmap.josm.data.projection; 009 010import org.openstreetmap.josm.data.coor.LatLon; 011 012/** 013 * Reference ellipsoids. 014 */ 015public final class Ellipsoid { 016 017 /** 018 * Clarke 1880 IGN (French national geographic institute) 019 */ 020 public static final Ellipsoid clarkeIGN = Ellipsoid.create_a_b(6378249.2, 6356515.0); 021 022 /** 023 * Hayford's ellipsoid 1909 (ED50 system)<br> 024 * Proj.4 code: intl 025 */ 026 public static final Ellipsoid hayford = Ellipsoid.create_a_rf(6378388.0, 297.0); 027 028 /** 029 * GRS67 ellipsoid 030 */ 031 public static final Ellipsoid GRS67 = Ellipsoid.create_a_rf(6378160.0, 298.247167472); 032 033 /** 034 * GRS80 ellipsoid 035 */ 036 public static final Ellipsoid GRS80 = Ellipsoid.create_a_rf(6378137.0, 298.257222101); 037 038 /** 039 * WGS84 ellipsoid 040 */ 041 public static final Ellipsoid WGS84 = Ellipsoid.create_a_rf(6378137.0, 298.257223563); 042 043 /** 044 * Bessel 1841 ellipsoid 045 */ 046 public static final Ellipsoid Bessel1841 = Ellipsoid.create_a_rf(6377397.155, 299.1528128); 047 048 /** 049 * half long axis 050 */ 051 public final double a; 052 053 /** 054 * half short axis 055 */ 056 public final double b; 057 058 /** 059 * first eccentricity 060 */ 061 public final double e; 062 063 /** 064 * first eccentricity squared 065 */ 066 public final double e2; 067 068 /** 069 * square of the second eccentricity 070 */ 071 public final double eb2; 072 073 /** 074 * private constructur - use one of the create_* methods 075 * 076 * @param a semimajor radius of the ellipsoid axis 077 * @param b semiminor radius of the ellipsoid axis 078 * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a))) 079 * @param e2 first eccentricity squared 080 * @param eb2 square of the second eccentricity 081 */ 082 private Ellipsoid(double a, double b, double e, double e2, double eb2) { 083 this.a = a; 084 this.b = b; 085 this.e = e; 086 this.e2 = e2; 087 this.eb2 = eb2; 088 } 089 090 /** 091 * create a new ellipsoid 092 * 093 * @param a semimajor radius of the ellipsoid axis (in meters) 094 * @param b semiminor radius of the ellipsoid axis (in meters) 095 * @return the new ellipsoid 096 */ 097 public static Ellipsoid create_a_b(double a, double b) { 098 double e2 = (a*a - b*b) / (a*a); 099 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}