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    }