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}