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