001// License: GPL. For details, see LICENSE file.
002package org.openstreetmap.josm.data.projection.proj;
003
004import org.openstreetmap.josm.data.projection.Ellipsoid;
005import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
006import org.openstreetmap.josm.tools.CheckParameterUtil;
007
008/**
009 * Abstract base class providing utilities for implementations of the Proj
010 * interface.
011 *
012 * This class has been derived from the implementation of the Geotools project;
013 * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection
014 * at the time of migration.
015 * <p>
016 *
017 * @author André Gosselin
018 * @author Martin Desruisseaux (PMO, IRD)
019 * @author Rueben Schulz
020*/
021public abstract class AbstractProj implements Proj {
022
023    /**
024     * Maximum number of iterations for iterative computations.
025     */
026    private static final int MAXIMUM_ITERATIONS = 15;
027
028    /**
029     * Difference allowed in iterative computations.
030     */
031    private static final double ITERATION_TOLERANCE = 1E-10;
032
033    /**
034     * Relative iteration precision used in the <code>mlfn</code> method
035     */
036    private static final double MLFN_TOL = 1E-11;
037
038    /**
039     * Constants used to calculate {@link #en0}, {@link #en1},
040     * {@link #en2}, {@link #en3}, {@link #en4}.
041     */
042    private static final double C00 = 1.0;
043    private static final double C02 = 0.25;
044    private static final double C04 = 0.046875;
045    private static final double C06 = 0.01953125;
046    private static final double C08 = 0.01068115234375;
047    private static final double C22 = 0.75;
048    private static final double C44 = 0.46875;
049    private static final double C46 = 0.01302083333333333333;
050    private static final double C48 = 0.00712076822916666666;
051    private static final double C66 = 0.36458333333333333333;
052    private static final double C68 = 0.00569661458333333333;
053    private static final double C88 = 0.3076171875;
054
055    /**
056     * Constant needed for the <code>mlfn</code> method.
057     * Setup at construction time.
058     */
059    protected double en0, en1, en2, en3, en4;
060
061    /**
062     * Ellipsoid excentricity, equals to <code>sqrt({@link #e2 excentricity squared})</code>.
063     * Value 0 means that the ellipsoid is spherical.
064     *
065     * @see #e2
066     */
067    protected double e;
068
069    /**
070     * The square of excentricity: e² = (a²-b²)/a² where
071     * <var>e</var> is the excentricity,
072     * <var>a</var> is the semi major axis length and
073     * <var>b</var> is the semi minor axis length.
074     *
075     * @see #e
076     */
077    protected double e2;
078
079    /**
080     * is ellipsoid spherical?
081     * @see Ellipsoid#spherical
082     */
083    protected boolean spherical;
084
085    @Override
086    public void initialize(ProjParameters params) throws ProjectionConfigurationException {
087        CheckParameterUtil.ensureParameterNotNull(params, "params");
088        CheckParameterUtil.ensureParameterNotNull(params.ellps, "params.ellps");
089        e2 = params.ellps.e2;
090        e = params.ellps.e;
091        spherical = params.ellps.spherical;
092        //  Compute constants for the mlfn
093        double t;
094        en0 = C00 - e2  *  (C02 + e2  *
095             (C04 + e2  *  (C06 + e2  * C08)));
096        en1 =       e2  *  (C22 - e2  *
097             (C04 + e2  *  (C06 + e2  * C08)));
098        en2 =  (t = e2  *         e2) *
099             (C44 - e2  *  (C46 + e2  * C48));
100        en3 = (t *= e2) *  (C66 - e2  * C68);
101        en4 =   t * e2  *  C88;
102    }
103
104    @Override
105    public boolean isGeographic() {
106        return false;
107    }
108
109    /**
110     * Calculates the meridian distance. This is the distance along the central
111     * meridian from the equator to {@code phi}. Accurate to &lt; 1e-5 meters
112     * when used in conjuction with typical major axis values.
113     *
114     * @param phi latitude to calculate meridian distance for.
115     * @param sphi sin(phi).
116     * @param cphi cos(phi).
117     * @return meridian distance for the given latitude.
118     */
119    protected final double mlfn(final double phi, double sphi, double cphi) {
120        cphi *= sphi;
121        sphi *= sphi;
122        return en0 * phi - cphi *
123              (en1 + sphi *
124              (en2 + sphi *
125              (en3 + sphi *
126              (en4))));
127    }
128
129    /**
130     * Calculates the latitude ({@code phi}) from a meridian distance.
131     * Determines phi to TOL (1e-11) radians, about 1e-6 seconds.
132     *
133     * @param arg meridian distance to calulate latitude for.
134     * @return the latitude of the meridian distance.
135     * @throws RuntimeException if the itteration does not converge.
136     */
137    protected final double inv_mlfn(double arg) {
138        double s, t, phi, k = 1.0/(1.0 - e2);
139        int i;
140        phi = arg;
141        for (i = MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations
142            if (--i < 0) {
143                throw new IllegalStateException("Too many iterations");
144            }
145            s = Math.sin(phi);
146            t = 1.0 - e2 * s * s;
147            t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k;
148            phi -= t;
149            if (Math.abs(t) < MLFN_TOL) {
150                return phi;
151            }
152        }
153    }
154
155    // Iteratively solve equation (7-9) from Snyder.
156    final double cphi2(final double ts) {
157        final double eccnth = 0.5 * e;
158        double phi = (Math.PI/2) - 2.0 * Math.atan(ts);
159        for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
160            final double con  = e * Math.sin(phi);
161            final double dphi = (Math.PI/2) - 2.0*Math.atan(ts * Math.pow((1-con)/(1+con), eccnth)) - phi;
162            phi += dphi;
163            if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
164                return phi;
165            }
166        }
167        throw new IllegalStateException("no convergence for ts="+ts);
168    }
169
170    /**
171     * Computes function <code>f(s,c,e²) = c/sqrt(1 - s²&times;e²)</code> needed for the true scale
172     * latitude (Snyder 14-15), where <var>s</var> and <var>c</var> are the sine and cosine of
173     * the true scale latitude, and <var>e²</var> is the {@linkplain #e2 eccentricity squared}.
174     * @param s sine of the true scale latitude
175     * @param c cosine of the true scale latitude
176     * @return <code>c/sqrt(1 - s²&times;e²)</code>
177     */
178    final double msfn(final double s, final double c) {
179        return c / Math.sqrt(1.0 - (s*s) * e2);
180    }
181
182    /**
183     * Computes function (15-9) and (9-13) from Snyder.
184     * Equivalent to negative of function (7-7).
185     * @param lat the latitude
186     * @param sinlat sine of the latitude
187     * @return auxiliary value computed from <code>lat</code> and <code>sinlat</code>
188     */
189    final double tsfn(final double lat, double sinlat) {
190        sinlat *= e;
191        // NOTE: change sign to get the equivalent of Snyder (7-7).
192        return Math.tan(0.5 * (Math.PI/2 - lat)) / Math.pow((1 - sinlat) / (1 + sinlat), 0.5*e);
193    }
194}