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.atan;
007import static java.lang.Math.cos;
008import static java.lang.Math.exp;
009import static java.lang.Math.log;
010import static java.lang.Math.pow;
011import static java.lang.Math.sin;
012import static java.lang.Math.sqrt;
013import static java.lang.Math.tan;
014import static java.lang.Math.toRadians;
015import static org.openstreetmap.josm.tools.I18n.tr;
016
017import org.openstreetmap.josm.data.projection.Ellipsoid;
018import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
019
020/**
021 * Implementation of the Lambert Conformal Conic projection.
022 *
023 * @author Pieren
024 */
025public class LambertConformalConic implements Proj {
026
027    protected Ellipsoid ellps;
028    protected double e;
029
030    public abstract static class Parameters {
031        public final double latitudeOrigin;
032        public Parameters(double latitudeOrigin) {
033            this.latitudeOrigin = latitudeOrigin;
034        }
035    }
036
037    public static class Parameters1SP extends Parameters {
038        public Parameters1SP(double latitudeOrigin) {
039            super(latitudeOrigin);
040        }
041    }
042
043    public static class Parameters2SP extends Parameters {
044        public final double standardParallel1;
045        public final double standardParallel2;
046        public Parameters2SP(double latitudeOrigin, double standardParallel1, double standardParallel2) {
047            super(latitudeOrigin);
048            this.standardParallel1 = standardParallel1;
049            this.standardParallel2 = standardParallel2;
050        }
051    }
052
053    private Parameters params;
054
055    /**
056     * projection exponent
057     */
058    protected double n;
059    /**
060     * projection factor
061     */
062    protected double F;
063    /**
064     * radius of the parallel of latitude of the false origin (2SP) or at
065     * natural origin (1SP)
066     */
067    protected double r0;
068
069    /**
070     * precision in iterative schema
071     */
072    protected static final double epsilon = 1e-12;
073
074    @Override
075    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
076        ellps = params.ellps;
077        e = ellps.e;
078        if (params.lat_0 == null)
079            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
080        if (params.lat_1 != null && params.lat_2 != null) {
081            initialize2SP(params.lat_0, params.lat_1, params.lat_2);
082        } else {
083            initialize1SP(params.lat_0);
084        }
085    }
086
087    /**
088     * Initialize for LCC with 2 standard parallels.
089     *
090     * @param lat_0 latitude of false origin (in degrees)
091     * @param lat_1 latitude of first standard parallel (in degrees)
092     * @param lat_2 latitude of second standard parallel (in degrees)
093     */
094    private void initialize2SP(double lat_0, double lat_1, double lat_2) {
095        this.params = new Parameters2SP(lat_0, lat_1, lat_2);
096
097        final double m1 = m(toRadians(lat_1));
098        final double m2 = m(toRadians(lat_2));
099
100        final double t1 = t(toRadians(lat_1));
101        final double t2 = t(toRadians(lat_2));
102        final double tf = t(toRadians(lat_0));
103
104        n  = (log(m1) - log(m2)) / (log(t1) - log(t2));
105        F  = m1 / (n * pow(t1, n));
106        r0 = F * pow(tf, n);
107    }
108
109    /**
110     * Initialize for LCC with 1 standard parallel.
111     *
112     * @param lat_0 latitude of natural origin (in degrees)
113     */
114    private void initialize1SP(double lat_0) {
115        this.params = new Parameters1SP(lat_0);
116        final double lat_0_rad = toRadians(lat_0);
117
118        final double m0 = m(lat_0_rad);
119        final double t0 = t(lat_0_rad);
120
121        n = sin(lat_0_rad);
122        F  = m0 / (n * pow(t0, n));
123        r0 = F * pow(t0, n);
124    }
125
126    /**
127     * auxiliary function t
128     */
129    protected double t(double lat_rad) {
130        return tan(PI/4 - lat_rad / 2.0)
131            / pow(( (1.0 - e * sin(lat_rad)) / (1.0 + e * sin(lat_rad))) , e/2);
132    }
133
134    /**
135     * auxiliary function m
136     */
137    protected double m(double lat_rad) {
138        return cos(lat_rad) / (sqrt(1 - e * e * pow(sin(lat_rad), 2)));
139    }
140
141    @Override
142    public String getName() {
143        return tr("Lambert Conformal Conic");
144    }
145
146    @Override
147    public String getProj4Id() {
148        return "lcc";
149    }
150
151    @Override
152    public double[] project(double phi, double lambda) {
153        double sinphi = sin(phi);
154        double L = (0.5*log((1+sinphi)/(1-sinphi))) - e/2*log((1+e*sinphi)/(1-e*sinphi));
155        double r = F*exp(-n*L);
156        double gamma = n*lambda;
157        double X = r*sin(gamma);
158        double Y = r0 - r*cos(gamma);
159        return new double[] { X, Y };
160    }
161
162    @Override
163    public double[] invproject(double east, double north) {
164        double r = sqrt(pow(east,2) + pow(north-r0, 2));
165        double gamma = atan(east / (r0-north));
166        double lambda = gamma/n;
167        double latIso = (-1/n) * log(abs(r/F));
168        double phi = ellps.latitude(latIso, e, epsilon);
169        return new double[] { phi, lambda };
170    }
171
172    public final Parameters getParameters() {
173        return params;
174    }
175}