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 }