001// License: GPL. For details, see LICENSE file.
002package org.openstreetmap.josm.data.projection.proj;
003
004import static org.openstreetmap.josm.tools.I18n.tr;
005
006import org.openstreetmap.josm.data.Bounds;
007import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
008import org.openstreetmap.josm.tools.Utils;
009
010/**
011 * Cassini-Soldner Projection (EPSG code 9806).
012 * The Cassini-Soldner Projection is the ellipsoidal version of the Cassini
013 * projection for the sphere. It is not conformal but as it is relatively simple
014 * to construct it was extensively used in the last century and is still useful
015 * for mapping areas with limited longitudinal extent. It has now largely
016 * been replaced by the conformal Transverse Mercator which it resembles. Like this,
017 * it has a straight central meridian along which the scale is true, all other
018 * meridians and parallels are curved, and the scale distortion increases
019 * rapidly with increasing distance from the central meridian.
020 * <p>
021 *
022 * This class has been derived from the implementation of the Geotools project;
023 * git 8cbf52d, org.geotools.referencing.operation.projection.CassiniSoldner
024 * at the time of migration.
025 */
026public class CassiniSoldner extends AbstractProj {
027
028    /**
029     * Meridian distance at the {@code latitudeOfOrigin}.
030     * Used for calculations for the ellipsoid.
031     */
032    private double ml0;
033
034    /**
035     * Latitude of origin.
036     */
037    private double phi0;
038
039    /**
040     * Constants used for the forward and inverse transform for the elliptical
041     * case of the Cassini-Soldner.
042     */
043    private static final double C1 = 0.16666666666666666666;
044    private static final double C2 = 0.008333333333333333333;
045    private static final double C3 = 0.041666666666666666666;
046    private static final double C4 = 0.33333333333333333333;
047    private static final double C5 = 0.066666666666666666666;
048
049    @Override
050    public String getName() {
051        return tr("Cassini-Soldner");
052    }
053
054    @Override
055    public String getProj4Id() {
056        return "cass";
057    }
058
059    @Override
060    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
061        super.initialize(params);
062        if (params.lat0 == null)
063            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
064        phi0 = Utils.toRadians(params.lat0);
065        ml0 = mlfn(phi0, Math.sin(phi0), Math.cos(phi0));
066    }
067
068    @Override
069    public double[] project(double phi, double lam) {
070        if (spherical) {
071            double x = aasin(Math.cos(phi) * Math.sin(lam));
072            double y = Math.atan2(Math.tan(phi), Math.cos(lam));
073            return new double[] {x, y};
074        } else {
075            double sinphi = Math.sin(phi);
076            double cosphi = Math.cos(phi);
077
078            double n = 1.0 / (Math.sqrt(1.0 - e2 * sinphi * sinphi));
079            double tn = Math.tan(phi);
080            double t = tn * tn;
081            double a1 = lam * cosphi;
082            double c = cosphi * cosphi * e2 / (1 - e2);
083            double a2 = a1 * a1;
084
085            double x = n * a1 * (1.0 - a2 * t * (C1 - (8.0 - t + 8.0 * c) * a2 * C2));
086            double y = mlfn(phi, sinphi, cosphi) - ml0 + n * tn * a2 * (0.5 + (5.0 - t + 6.0 * c) * a2 * C3);
087            return new double[] {x, y};
088        }
089    }
090
091    @Override
092    public double[] invproject(double x, double y) {
093        if (spherical) {
094            double dd = y + phi0;
095            double phi = aasin(Math.sin(dd * Math.cos(x)));
096            double lam = Math.atan2(Math.tan(x), Math.cos(dd));
097            return new double[] {phi, lam};
098        } else {
099            double ph1 = invMlfn(ml0 + y);
100            double tn = Math.tan(ph1);
101            double t = tn * tn;
102            double n = Math.sin(ph1);
103            double r = 1.0 / (1.0 - e2 * n * n);
104            n = Math.sqrt(r);
105            r *= (1.0 - e2) * n;
106            double dd = x / n;
107            double d2 = dd * dd;
108            double phi = ph1 - (n * tn / r) * d2 * (0.5 - (1.0 + 3.0 * t) * d2 * C3);
109            double lam = dd * (1.0 + t * d2 * (-C4 + (1.0 + 3.0 * t) * d2 * C5)) / Math.cos(ph1);
110            return new double[] {phi, lam};
111        }
112    }
113
114    @Override
115    public Bounds getAlgorithmBounds() {
116        return new Bounds(-89, -1.0, 89, 1.0, false);
117    }
118}