001//License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import static java.lang.Math.PI; 005import static java.lang.Math.abs; 006import static java.lang.Math.asin; 007import static java.lang.Math.atan; 008import static java.lang.Math.atan2; 009import static java.lang.Math.cos; 010import static java.lang.Math.exp; 011import static java.lang.Math.log; 012import static java.lang.Math.pow; 013import static java.lang.Math.sin; 014import static java.lang.Math.sqrt; 015import static java.lang.Math.tan; 016import static java.lang.Math.toRadians; 017import static org.openstreetmap.josm.tools.I18n.tr; 018 019import org.openstreetmap.josm.data.projection.Ellipsoid; 020import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 021 022/** 023 * Projection for the SwissGrid CH1903 / L03, see http://en.wikipedia.org/wiki/Swiss_coordinate_system. 024 * 025 * Calculations were originally based on simple formula from 026 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf 027 * 028 * August 2010 update to this formula (rigorous formulas) 029 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf 030 */ 031public class SwissObliqueMercator implements Proj { 032 033 private Ellipsoid ellps; 034 private double kR; 035 private double alpha; 036 private double b0; 037 private double K; 038 039 private static final double EPSILON = 1e-11; 040 041 @Override 042 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 043 if (params.lat_0 == null) 044 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0")); 045 ellps = params.ellps; 046 initialize(params.lat_0); 047 } 048 049 private void initialize(double lat_0) { 050 double phi0 = toRadians(lat_0); 051 kR = sqrt(1 - ellps.e2) / (1 - (ellps.e2 * pow(sin(phi0), 2))); 052 alpha = sqrt(1 + (ellps.eb2 * pow(cos(phi0), 4))); 053 b0 = asin(sin(phi0) / alpha); 054 K = log(tan(PI / 4 + b0 / 2)) - alpha 055 * log(tan(PI / 4 + phi0 / 2)) + alpha * ellps.e / 2 056 * log((1 + ellps.e * sin(phi0)) / (1 - ellps.e * sin(phi0))); 057 } 058 059 @Override 060 public String getName() { 061 return tr("Swiss Oblique Mercator"); 062 } 063 064 @Override 065 public String getProj4Id() { 066 return "somerc"; 067 } 068 069 @Override 070 public double[] project(double phi, double lambda) { 071 072 double S = alpha * log(tan(PI / 4 + phi / 2)) - alpha * ellps.e / 2 073 * log((1 + ellps.e * sin(phi)) / (1 - ellps.e * sin(phi))) + K; 074 double b = 2 * (atan(exp(S)) - PI / 4); 075 double l = alpha * lambda; 076 077 double lb = atan2(sin(l), sin(b0) * tan(b) + cos(b0) * cos(l)); 078 double bb = asin(cos(b0) * sin(b) - sin(b0) * cos(b) * cos(l)); 079 080 double y = kR * lb; 081 double x = kR / 2 * log((1 + sin(bb)) / (1 - sin(bb))); 082 083 return new double[] { y, x }; 084 } 085 086 @Override 087 public double[] invproject(double y, double x) { 088 double lb = y / kR; 089 double bb = 2 * (atan(exp(x / kR)) - PI / 4); 090 091 double b = asin(cos(b0) * sin(bb) + sin(b0) * cos(bb) * cos(lb)); 092 double l = atan2(sin(lb), cos(b0) * cos(lb) - sin(b0) * tan(bb)); 093 094 double lambda = l / alpha; 095 double phi = b; 096 double S = 0; 097 098 double prevPhi = -1000; 099 int iteration = 0; 100 // iteration to finds S and phi 101 while (abs(phi - prevPhi) > EPSILON) { 102 if (++iteration > 30) 103 throw new RuntimeException("Two many iterations"); 104 prevPhi = phi; 105 S = 1 / alpha * (log(tan(PI / 4 + b / 2)) - K) + ellps.e 106 * log(tan(PI / 4 + asin(ellps.e * sin(phi)) / 2)); 107 phi = 2 * atan(exp(S)) - PI / 2; 108 } 109 return new double[] { phi, lambda }; 110 } 111 112}