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.cos;
009import static java.lang.Math.exp;
010import static java.lang.Math.log;
011import static java.lang.Math.pow;
012import static java.lang.Math.sin;
013import static java.lang.Math.sqrt;
014import static java.lang.Math.tan;
015import static java.lang.Math.toRadians;
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;
021
022/**
023 * Implementation of the stereographic double projection,
024 * also known as Oblique Stereographic and the Schreiber double stereographic projection.
025 *
026 * @author vholten
027 *
028 * Source: IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2,
029 * Sec. 1.3.7.1 Oblique and Equatorial Stereographic, http://www.epsg.org/GuidanceNotes
030 */
031public class DoubleStereographic implements Proj {
032
033    private Ellipsoid ellps;
034    private double e;
035    private double n;
036    private double c;
037    private double chi0;
038    private double R;
039
040    private static final double EPSILON = 1e-12;
041
042    @Override
043    public String getName() {
044        return tr("Double Stereographic");
045    }
046
047    @Override
048    public String getProj4Id() {
049        return "sterea";
050    }
051
052    @Override
053    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
054        if (params.lat0 == null)
055            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
056        ellps = params.ellps;
057        e = ellps.e;
058        initialize(params.lat0);
059    }
060
061    private void initialize(double lat_0) {
062        double phi0 = toRadians(lat_0);
063        double e2 = ellps.e2;
064        R = sqrt(1-e2) / (1 - e2*pow(sin(phi0), 2));
065        n = sqrt(1 + ellps.eb2 * pow(cos(phi0), 4));
066        double S1 = (1 + sin(phi0)) / (1 - sin(phi0));
067        double S2 = (1 - e * sin(phi0)) / (1 + e * sin(phi0));
068        double w1 = pow(S1 * pow(S2, e), n);
069        double sinchi00 = (w1 - 1) / (w1 + 1);
070        c = (n + sin(phi0)) * (1 - sinchi00) / ((n - sin(phi0)) * (1 + sinchi00));
071        double w2 = c * w1;
072        chi0 = asin((w2 - 1) / (w2 + 1));
073    }
074
075    @Override
076    public double[] project(double phi, double lambda) {
077        double Lambda = n * lambda;
078        double Sa = (1 + sin(phi)) / (1 - sin(phi));
079        double Sb = (1 - e * sin(phi)) / (1 + e * sin(phi));
080        double w = c * pow(Sa * pow(Sb, e), n);
081        double chi = asin((w - 1) / (w + 1));
082        double B = 1 + sin(chi) * sin(chi0) + cos(chi) * cos(chi0) * cos(Lambda);
083        double x = 2 * R * cos(chi) * sin(Lambda) / B;
084        double y = 2 * R * (sin(chi) * cos(chi0) - cos(chi) * sin(chi0) * cos(Lambda)) / B;
085        return new double[] {x, y};
086    }
087
088    @Override
089    public double[] invproject(double x, double y) {
090        double e2 = ellps.e2;
091        double g = 2 * R * tan(PI/4 - chi0/2);
092        double h = 4 * R * tan(chi0) + g;
093        double i = atan(x/(h + y));
094        double j = atan(x/(g - y)) - i;
095        double chi = chi0 + 2 * atan((y - x * tan(j/2)) / (2 * R));
096        double Lambda = j + 2*i;
097        double lambda = Lambda / n;
098        double psi = 0.5 * log((1 + sin(chi)) / (c*(1 - sin(chi)))) / n;
099        double phiprev = -1000;
100        int iteration = 0;
101        double phi = 2 * atan(exp(psi)) - PI/2;
102        while (abs(phi - phiprev) > EPSILON) {
103            if (++iteration > 10)
104                throw new RuntimeException("Too many iterations");
105            phiprev = phi;
106            double psii = log(tan(phi/2 + PI/4) * pow((1 - e * sin(phi)) / (1 + e * sin(phi)), e/2));
107            phi = phi - (psii - psi) * cos(phi) * (1 - e2 * pow(sin(phi), 2)) / (1 - e2);
108        }
109        return new double[] {phi, lambda};
110    }
111
112    @Override
113    public Bounds getAlgorithmBounds() {
114        return new Bounds(-89, -87, 89, 87, false);
115    }
116}