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