001// License: GPL. For details, see LICENSE file.
002package org.openstreetmap.josm.data.projection.proj;
003
004import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
005
006/**
007 * Abstract base class providing utilities for implementations of the Proj
008 * interface.
009 *
010 * This class has been derived from the implementation of the Geotools project;
011 * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection
012 * at the time of migration.
013 * <p>
014 *
015 * @author André Gosselin
016 * @author Martin Desruisseaux (PMO, IRD)
017 * @author Rueben Schulz
018*/
019public abstract class AbstractProj implements Proj {
020
021    /**
022     * Maximum number of iterations for iterative computations.
023     */
024    private static final int MAXIMUM_ITERATIONS = 15;
025
026    /**
027     * Relative iteration precision used in the <code>mlfn</code> method
028     */
029    private static final double MLFN_TOL = 1E-11;
030
031    /**
032     * Constants used to calculate {@link #en0}, {@link #en1},
033     * {@link #en2}, {@link #en3}, {@link #en4}.
034     */
035    private static final double C00 = 1.0,
036                                C02 = 0.25,
037                                C04 = 0.046875,
038                                C06 = 0.01953125,
039                                C08 = 0.01068115234375,
040                                C22 = 0.75,
041                                C44 = 0.46875,
042                                C46 = 0.01302083333333333333,
043                                C48 = 0.00712076822916666666,
044                                C66 = 0.36458333333333333333,
045                                C68 = 0.00569661458333333333,
046                                C88 = 0.3076171875;
047
048    /**
049     * Constant needed for the <code>mlfn</code> method.
050     * Setup at construction time.
051     */
052    protected double en0, en1, en2, en3, en4;
053
054    /**
055     * The square of excentricity: e² = (a²-b²)/a² where
056     * <var>e</var> is the excentricity,
057     * <var>a</var> is the semi major axis length and
058     * <var>b</var> is the semi minor axis length.
059     */
060    protected double e2;
061
062    @Override
063    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
064        e2 = params.ellps.e2;
065        //  Compute constants for the mlfn
066        double t;
067        en0 = C00 - e2  *  (C02 + e2  *
068             (C04 + e2  *  (C06 + e2  * C08)));
069        en1 =       e2  *  (C22 - e2  *
070             (C04 + e2  *  (C06 + e2  * C08)));
071        en2 =  (t = e2  *         e2) *
072             (C44 - e2  *  (C46 + e2  * C48));
073        en3 = (t *= e2) *  (C66 - e2  * C68);
074        en4 =   t * e2  *  C88;
075    }
076
077    /**
078     * Calculates the meridian distance. This is the distance along the central
079     * meridian from the equator to {@code phi}. Accurate to &lt; 1e-5 meters
080     * when used in conjuction with typical major axis values.
081     *
082     * @param phi latitude to calculate meridian distance for.
083     * @param sphi sin(phi).
084     * @param cphi cos(phi).
085     * @return meridian distance for the given latitude.
086     */
087    protected final double mlfn(final double phi, double sphi, double cphi) {
088        cphi *= sphi;
089        sphi *= sphi;
090        return en0 * phi - cphi *
091              (en1 + sphi *
092              (en2 + sphi *
093              (en3 + sphi *
094              (en4))));
095    }
096
097    /**
098     * Calculates the latitude ({@code phi}) from a meridian distance.
099     * Determines phi to TOL (1e-11) radians, about 1e-6 seconds.
100     *
101     * @param arg meridian distance to calulate latitude for.
102     * @return the latitude of the meridian distance.
103     * @throws RuntimeException if the itteration does not converge.
104     */
105    protected final double inv_mlfn(double arg) {
106        double s, t, phi, k = 1.0/(1.0 - e2);
107        int i;
108        phi = arg;
109        for (i = MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations
110            if (--i < 0) {
111                throw new RuntimeException("Too many iterations");
112            }
113            s = Math.sin(phi);
114            t = 1.0 - e2 * s * s;
115            t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k;
116            phi -= t;
117            if (Math.abs(t) < MLFN_TOL) {
118                return phi;
119            }
120        }
121    }
122
123    public static double normalizeLon(double lon) {
124        if (lon >= -Math.PI && lon <= Math.PI)
125            return lon;
126        else {
127            lon = lon % (2 * Math.PI);
128            if (lon > Math.PI) {
129                return lon - 2 * Math.PI;
130            } else if (lon < -Math.PI) {
131                return lon + 2 * Math.PI;
132            }
133            return lon;
134        }
135    }
136}