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 * the reference ellipsoids 014 */ 015public final class Ellipsoid { 016 /** 017 * Clarke 1880 IGN (French national geographic institute) 018 */ 019 public static final Ellipsoid clarkeIGN = Ellipsoid.create_a_b(6378249.2, 6356515.0); 020 /** 021 * Hayford's ellipsoid 1909 (ED50 system) 022 * Proj.4 code: intl 023 */ 024 public static final Ellipsoid hayford = Ellipsoid.create_a_rf(6378388.0, 297.0); 025 /** 026 * GRS67 ellipsoid 027 */ 028 public static final Ellipsoid GRS67 = Ellipsoid.create_a_rf(6378160.0, 298.247167472); 029 /** 030 * GRS80 ellipsoid 031 */ 032 public static final Ellipsoid GRS80 = Ellipsoid.create_a_rf(6378137.0, 298.257222101); 033 034 /** 035 * WGS84 ellipsoid 036 */ 037 public static final Ellipsoid WGS84 = Ellipsoid.create_a_rf(6378137.0, 298.257223563); 038 039 /** 040 * Bessel 1841 ellipsoid 041 */ 042 public static final Ellipsoid Bessel1841 = Ellipsoid.create_a_rf(6377397.155, 299.1528128); 043 044 /** 045 * half long axis 046 */ 047 public final double a; 048 /** 049 * half short axis 050 */ 051 public final double b; 052 /** 053 * first eccentricity 054 */ 055 public final double e; 056 /** 057 * first eccentricity squared 058 */ 059 public final double e2; 060 061 /** 062 * square of the second eccentricity 063 */ 064 public final double eb2; 065 066 /** 067 * private constructur - use one of the create_* methods 068 * 069 * @param a semimajor radius of the ellipsoid axis 070 * @param b semiminor radius of the ellipsoid axis 071 * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a))) 072 * @param e2 first eccentricity squared 073 * @param eb2 square of the second eccentricity 074 */ 075 private Ellipsoid(double a, double b, double e, double e2, double eb2) { 076 this.a = a; 077 this.b = b; 078 this.e = e; 079 this.e2 = e2; 080 this.eb2 = eb2; 081 } 082 083 /** 084 * create a new ellipsoid 085 * 086 * @param a semimajor radius of the ellipsoid axis (in meters) 087 * @param b semiminor radius of the ellipsoid axis (in meters) 088 */ 089 public static Ellipsoid create_a_b(double a, double b) { 090 double e2 = (a*a - b*b) / (a*a); 091 double e = Math.sqrt(e2); 092 double eb2 = e2 / (1.0 - e2); 093 return new Ellipsoid(a, b, e, e2, eb2); 094 } 095 096 /** 097 * create a new ellipsoid 098 * 099 * @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 246 public LatLon cart2LatLon(double[] xyz, double epsilon) { 247 double norm = Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]); 248 double lg = 2.0 * Math.atan(xyz[1] / (xyz[0] + norm)); 249 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]))))); 250 double delta = 1.0; 251 while (delta > epsilon) { 252 double s2 = Math.sin(lt); 253 s2 *= s2; 254 double l = Math.atan((xyz[2] / norm) 255 / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2))))); 256 delta = Math.abs(l - lt); 257 lt = l; 258 } 259 return new LatLon(Math.toDegrees(lt), Math.toDegrees(lg)); 260 } 261 262 /** 263 * convert ellipsoidal coordinates to cartesian coordinates 264 * 265 * @param coord The Latitude and longitude in degrees 266 * @return the corresponding (X, Y Z) cartesian coordinates in meters. 267 */ 268 public double[] latLon2Cart(LatLon coord) { 269 double phi = Math.toRadians(coord.lat()); 270 double lambda = Math.toRadians(coord.lon()); 271 272 double Rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2)); 273 double[] xyz = new double[3]; 274 xyz[0] = Rn * Math.cos(phi) * Math.cos(lambda); 275 xyz[1] = Rn * Math.cos(phi) * Math.sin(lambda); 276 xyz[2] = Rn * (1 - e2) * Math.sin(phi); 277 278 return xyz; 279 } 280}