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