Changeset 9974 in josm for trunk/src/org/openstreetmap/josm/data/projection/proj/ObliqueMercator.java
- Timestamp:
- 2016-03-13T00:16:36+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/org/openstreetmap/josm/data/projection/proj/ObliqueMercator.java
r9970 r9974 159 159 * Constants used in the transformation. 160 160 */ 161 private double B, A, E;162 163 /** 164 * Convenience values equal to {@link # A} / {@link #B},165 * {@link # A}×{@link #B}, and {@link #B} / {@link #A}.166 */ 167 private double ArB, AB, BrA;161 private double b, a, e; 162 163 /** 164 * Convenience values equal to {@link #a} / {@link #b}, 165 * {@link #a}×{@link #b}, and {@link #b} / {@link #a}. 166 */ 167 private double arb, ab, bra; 168 168 169 169 /** 170 170 * <var>v</var> values when the input latitude is a pole. 171 171 */ 172 private double v _pole_n, v_pole_s;172 private double vPoleN, vPoleS; 173 173 174 174 /** … … 189 189 * the oblique mercator case. The <var>v</var> value of the central point is 0.0. 190 190 */ 191 private double u _c;191 private double uc; 192 192 193 193 /** … … 227 227 final double con = 1. - e2 * sinph0 * sinph0; 228 228 double temp = cosph0 * cosph0; 229 B= Math.sqrt(1.0 + e2 * (temp * temp) / (1.0 - e2));230 A = B* com / con;231 final double D = B* com / (cosph0 * Math.sqrt(con));232 double F = D * D- 1.0;233 if ( F< 0.0) {234 F= 0.0;229 b = Math.sqrt(1.0 + e2 * (temp * temp) / (1.0 - e2)); 230 a = b * com / con; 231 final double d = b * com / (cosph0 * Math.sqrt(con)); 232 double f = d * d - 1.0; 233 if (f < 0.0) { 234 f = 0.0; 235 235 } else { 236 F = Math.sqrt(F);236 f = Math.sqrt(f); 237 237 if (latCenter < 0.0) { 238 F = -F;239 } 240 } 241 E = F += D;242 E = F * Math.pow(tsfn(latCenter, sinph0), B);238 f = -f; 239 } 240 } 241 e = f += d; 242 e = f * Math.pow(tsfn(latCenter, sinph0), b); 243 243 244 244 /* … … 276 276 * The coefficients for the "two points" case. 277 277 */ 278 final double H = Math.pow(tsfn(lat1, Math.sin(lat1)), B);279 final double L = Math.pow(tsfn(lat2, Math.sin(lat2)), B);280 final double Fp = E / H;281 final double P = (L - H) / (L + H);282 double J = E * E;283 J = (J - L * H) / (J + L * H);278 final double h = Math.pow(tsfn(lat1, Math.sin(lat1)), b); 279 final double l = Math.pow(tsfn(lat2, Math.sin(lat2)), b); 280 final double fp = e / h; 281 final double p = (l - h) / (l + h); 282 double j = e * e; 283 j = (j - l * h) / (j + l * h); 284 284 double diff = lon1 - lon2; 285 285 if (diff < -Math.PI) { … … 289 289 } 290 290 centralMeridian = normalizeLonRad(0.5 * (lon1 + lon2) - 291 Math.atan( J * Math.tan(0.5 * B * (lon1 - lon2)) / P) / B);292 gamma0 = Math.atan(2.0 * Math.sin( B* normalizeLonRad(lon1 - centralMeridian)) /293 ( Fp - 1.0 / Fp));294 azimuth = Math.asin( D* Math.sin(gamma0));291 Math.atan(j * Math.tan(0.5 * b * (lon1 - lon2)) / p) / b); 292 gamma0 = Math.atan(2.0 * Math.sin(b * normalizeLonRad(lon1 - centralMeridian)) / 293 (fp - 1.0 / fp)); 294 azimuth = Math.asin(d * Math.sin(gamma0)); 295 295 rectifiedGridAngle = azimuth; 296 296 } else { … … 315 315 rectifiedGridAngle = azimuth; 316 316 } 317 gamma0 = Math.asin(Math.sin(azimuth) / D);317 gamma0 = Math.asin(Math.sin(azimuth) / d); 318 318 // Check for asin(+-1.00000001) 319 temp = 0.5 * ( F - 1.0 / F) * Math.tan(gamma0);319 temp = 0.5 * (f - 1.0 / f) * Math.tan(gamma0); 320 320 if (Math.abs(temp) > 1.0) { 321 321 if (Math.abs(Math.abs(temp) - 1.0) > EPSILON) { … … 324 324 temp = (temp > 0) ? 1.0 : -1.0; 325 325 } 326 centralMeridian = lonCenter - Math.asin(temp) / B;326 centralMeridian = lonCenter - Math.asin(temp) / b; 327 327 } 328 328 … … 334 334 sinrot = Math.sin(rectifiedGridAngle); 335 335 cosrot = Math.cos(rectifiedGridAngle); 336 ArB = A / B;337 AB = A * B;338 BrA = B / A;339 v _pole_n = ArB* Math.log(Math.tan(0.5 * (Math.PI/2.0 - gamma0)));340 v _pole_s = ArB* Math.log(Math.tan(0.5 * (Math.PI/2.0 + gamma0)));336 arb = a / b; 337 ab = a * b; 338 bra = b / a; 339 vPoleN = arb * Math.log(Math.tan(0.5 * (Math.PI/2.0 - gamma0))); 340 vPoleS = arb * Math.log(Math.tan(0.5 * (Math.PI/2.0 + gamma0))); 341 341 boolean hotine = params.no_off != null && params.no_off; 342 342 if (hotine) { 343 u _c = 0.0;343 uc = 0.0; 344 344 } else { 345 345 if (Math.abs(Math.abs(azimuth) - Math.PI/2.0) < EPSILON_LATITUDE) { … … 348 348 throw new ProjectionConfigurationException("assertion error"); 349 349 } 350 u _c = A* (lonCenter - centralMeridian);350 uc = a * (lonCenter - centralMeridian); 351 351 } else { 352 double u _c = Math.abs(ArB * Math.atan2(Math.sqrt(D * D- 1.0), Math.cos(azimuth)));352 double uC = Math.abs(arb * Math.atan2(Math.sqrt(d * d - 1.0), Math.cos(azimuth))); 353 353 if (latCenter < 0.0) { 354 u _c = -u_c;354 uC = -uC; 355 355 } 356 this.u _c = u_c;356 this.uc = uC; 357 357 } 358 358 } … … 367 367 double u, v; 368 368 if (Math.abs(Math.abs(y) - Math.PI/2.0) > EPSILON) { 369 double Q = E / Math.pow(tsfn(y, Math.sin(y)), B);370 double temp = 1.0 / Q;371 double S = 0.5 * (Q- temp);372 double V = Math.sin( B* x);373 double U = ( S * singamma0 - V * cosgamma0) / (0.5 * (Q+ temp));369 double q = e / Math.pow(tsfn(y, Math.sin(y)), b); 370 double temp = 1.0 / q; 371 double s = 0.5 * (q - temp); 372 double V = Math.sin(b * x); 373 double U = (s * singamma0 - V * cosgamma0) / (0.5 * (q + temp)); 374 374 if (Math.abs(Math.abs(U) - 1.0) < EPSILON) { 375 375 v = 0; // this is actually an error and should be reported to the caller somehow 376 376 } else { 377 v = 0.5 * ArB* Math.log((1.0 - U) / (1.0 + U));378 } 379 temp = Math.cos( B* x);377 v = 0.5 * arb * Math.log((1.0 - U) / (1.0 + U)); 378 } 379 temp = Math.cos(b * x); 380 380 if (Math.abs(temp) < EPSILON_LATITUDE) { 381 u = AB* x;381 u = ab * x; 382 382 } else { 383 u = ArB * Math.atan2(S* cosgamma0 + V * singamma0, temp);383 u = arb * Math.atan2(s * cosgamma0 + V * singamma0, temp); 384 384 } 385 385 } else { 386 v = y > 0 ? v _pole_n : v_pole_s;387 u = ArB* y;388 } 389 u -= u _c;386 v = y > 0 ? vPoleN : vPoleS; 387 u = arb * y; 388 } 389 u -= uc; 390 390 x = v * cosrot + u * sinrot; 391 391 y = u * cosrot - v * sinrot; … … 396 396 public double[] invproject(double x, double y) { 397 397 double v = x * cosrot - y * sinrot; 398 double u = y * cosrot + x * sinrot + u _c;399 double Qp = Math.exp(-BrA* v);400 double temp = 1.0 / Qp;401 double Sp = 0.5 * (Qp - temp);402 double Vp = Math.sin(BrA* u);403 double Up = (Vp * cosgamma0 + Sp * singamma0) / (0.5 * (Qp + temp));404 if (Math.abs(Math.abs( Up) - 1.0) < EPSILON) {398 double u = y * cosrot + x * sinrot + uc; 399 double qp = Math.exp(-bra * v); 400 double temp = 1.0 / qp; 401 double sp = 0.5 * (qp - temp); 402 double vp = Math.sin(bra * u); 403 double up = (vp * cosgamma0 + sp * singamma0) / (0.5 * (qp + temp)); 404 if (Math.abs(Math.abs(up) - 1.0) < EPSILON) { 405 405 x = 0.0; 406 y = Up < 0.0 ? -Math.PI / 2.0 : Math.PI / 2.0;406 y = up < 0.0 ? -Math.PI / 2.0 : Math.PI / 2.0; 407 407 } else { 408 y = Math.pow( E / Math.sqrt((1. + Up) / (1. - Up)), 1.0 / B); //calculate t408 y = Math.pow(e / Math.sqrt((1. + up) / (1. - up)), 1.0 / b); //calculate t 409 409 y = cphi2(y); 410 x = -Math.atan2( Sp * cosgamma0 - Vp * singamma0, Math.cos(BrA * u)) / B;410 x = -Math.atan2(sp * cosgamma0 - vp * singamma0, Math.cos(bra * u)) / b; 411 411 } 412 412 return new double[] {y, x};
Note:
See TracChangeset
for help on using the changeset viewer.