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 org.openstreetmap.josm.tools.I18n.tr;
017
018import org.openstreetmap.josm.data.Bounds;
019import org.openstreetmap.josm.data.projection.Ellipsoid;
020import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
021import org.openstreetmap.josm.tools.JosmRuntimeException;
022import org.openstreetmap.josm.tools.Utils;
023
024// CHECKSTYLE.OFF: LineLength
025
026/**
027 * Projection for the SwissGrid CH1903 / L03, see <a href="https://en.wikipedia.org/wiki/Swiss_coordinate_system">Wikipedia article</a>.<br>
028 *
029 * Calculations were originally based on
030 * <a href="http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf">
031 * simple formula</a>.<br>
032 *
033 * August 2010 update to
034 * <a href="http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf">
035 * this formula (rigorous formulas)</a>.
036 */
037public class SwissObliqueMercator extends AbstractProj {
038
039    // CHECKSTYLE.ON: LineLength
040
041    private Ellipsoid ellps;
042    private double kR;
043    private double alpha;
044    private double b0;
045    private double k;
046    private double phi0;
047
048    private static final double EPSILON = 1e-11;
049
050    @Override
051    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
052        super.initialize(params);
053        if (params.lat0 == null)
054            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
055        ellps = params.ellps;
056        initialize(params.lat0);
057    }
058
059    private void initialize(double lat0) {
060        phi0 = Utils.toRadians(lat0);
061        kR = sqrt(1 - ellps.e2) / (1 - (ellps.e2 * pow(sin(phi0), 2)));
062        alpha = sqrt(1 + (ellps.eb2 * pow(cos(phi0), 4)));
063        b0 = asin(sin(phi0) / alpha);
064        k = log(tan(PI / 4 + b0 / 2)) - alpha
065            * log(tan(PI / 4 + phi0 / 2)) + alpha * ellps.e / 2
066            * log((1 + ellps.e * sin(phi0)) / (1 - ellps.e * sin(phi0)));
067    }
068
069    @Override
070    public String getName() {
071        return tr("Swiss Oblique Mercator");
072    }
073
074    @Override
075    public String getProj4Id() {
076        return "somerc";
077    }
078
079    @Override
080    public double[] project(double phi, double lambda) {
081        double s = alpha * log(tan(PI / 4 + phi / 2)) - alpha * ellps.e / 2
082            * log((1 + ellps.e * sin(phi)) / (1 - ellps.e * sin(phi))) + k;
083        double b = 2 * (atan(exp(s)) - PI / 4);
084        double l = alpha * lambda;
085
086        double lb = atan2(sin(l), sin(b0) * tan(b) + cos(b0) * cos(l));
087        double bb = asin(cos(b0) * sin(b) - sin(b0) * cos(b) * cos(l));
088
089        double y = kR * lb;
090        double x = kR / 2 * log((1 + sin(bb)) / (1 - sin(bb)));
091
092        return new double[] {y, x};
093    }
094
095    @Override
096    public double[] invproject(double y, double x) {
097        double lb = y / kR;
098        double bb = 2 * (atan(exp(x / kR)) - PI / 4);
099
100        double b = asin(cos(b0) * sin(bb) + sin(b0) * cos(bb) * cos(lb));
101        double l = atan2(sin(lb), cos(b0) * cos(lb) - sin(b0) * tan(bb));
102
103        double lambda = l / alpha;
104        double phi = b;
105
106        double prevPhi = -1000;
107        int iteration = 0;
108        // iteration to finds S and phi
109        while (abs(phi - prevPhi) > EPSILON) {
110            if (++iteration > 30)
111                throw new JosmRuntimeException("Too many iterations");
112            prevPhi = phi;
113            double s = 1 / alpha * (log(tan(PI / 4 + b / 2)) - k) + ellps.e
114            * log(tan(PI / 4 + asin(ellps.e * sin(phi)) / 2));
115            phi = 2 * atan(exp(s)) - PI / 2;
116        }
117        return new double[] {phi, lambda};
118    }
119
120    @Override
121    public Bounds getAlgorithmBounds() {
122        if (phi0 > 0) {
123            return new Bounds(-10, -40, 85, 40, false);
124        } else {
125            return new Bounds(-85, -40, 10, 40, false);
126        }
127    }
128}