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        // CHECKSTYLE.OFF: SingleSpaceSeparator
095        en0 = C00 - e2  *  (C02 + e2  *
096             (C04 + e2  *  (C06 + e2  * C08)));
097        en1 =       e2  *  (C22 - e2  *
098             (C04 + e2  *  (C06 + e2  * C08)));
099        en2 =  (t = e2  *         e2) *
100             (C44 - e2  *  (C46 + e2  * C48));
101        en3 = (t *= e2) *  (C66 - e2  * C68);
102        en4 =   t * e2  *  C88;
103        // CHECKSTYLE.ON: SingleSpaceSeparator
104    }
105
106    @Override
107    public boolean isGeographic() {
108        return false;
109    }
110
111    /**
112     * Calculates the meridian distance. This is the distance along the central
113     * meridian from the equator to {@code phi}. Accurate to &lt; 1e-5 meters
114     * when used in conjunction with typical major axis values.
115     *
116     * @param phi latitude to calculate meridian distance for.
117     * @param sphi sin(phi).
118     * @param cphi cos(phi).
119     * @return meridian distance for the given latitude.
120     */
121    protected final double mlfn(final double phi, double sphi, double cphi) {
122        cphi *= sphi;
123        sphi *= sphi;
124        return en0 * phi - cphi *
125              (en1 + sphi *
126              (en2 + sphi *
127              (en3 + sphi *
128              (en4))));
129    }
130
131    /**
132     * Calculates the latitude ({@code phi}) from a meridian distance.
133     * Determines phi to TOL (1e-11) radians, about 1e-6 seconds.
134     *
135     * @param arg meridian distance to calculate latitude for.
136     * @return the latitude of the meridian distance.
137     * @throws RuntimeException if the itteration does not converge.
138     */
139    protected final double invMlfn(double arg) {
140        double s, t, phi, k = 1.0/(1.0 - e2);
141        int i;
142        phi = arg;
143        for (i = MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations
144            if (--i < 0) {
145                throw new IllegalStateException("Too many iterations");
146            }
147            s = Math.sin(phi);
148            t = 1.0 - e2 * s * s;
149            t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k;
150            phi -= t;
151            if (Math.abs(t) < MLFN_TOL) {
152                return phi;
153            }
154        }
155    }
156
157    /**
158     * Tolerant asin that will just return the limits of its output range if the input is out of range
159     * @param v the value whose arc sine is to be returned.
160     * @return the arc sine of the argument.
161     */
162    protected final double aasin(double v) {
163        double av = Math.abs(v);
164        if (av >= 1.) {
165            return (v < 0. ? -Math.PI / 2 : Math.PI / 2);
166        }
167        return Math.asin(v);
168    }
169
170    // Iteratively solve equation (7-9) from Snyder.
171    final double cphi2(final double ts) {
172        final double eccnth = 0.5 * e;
173        double phi = (Math.PI/2) - 2.0 * Math.atan(ts);
174        for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
175            final double con = e * Math.sin(phi);
176            final double dphi = (Math.PI/2) - 2.0*Math.atan(ts * Math.pow((1-con)/(1+con), eccnth)) - phi;
177            phi += dphi;
178            if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
179                return phi;
180            }
181        }
182        throw new IllegalStateException("no convergence for ts="+ts);
183    }
184
185    /**
186     * Computes function <code>f(s,c,e²) = c/sqrt(1 - s²&times;e²)</code> needed for the true scale
187     * latitude (Snyder 14-15), where <var>s</var> and <var>c</var> are the sine and cosine of
188     * the true scale latitude, and <var>e²</var> is the {@linkplain #e2 eccentricity squared}.
189     * @param s sine of the true scale latitude
190     * @param c cosine of the true scale latitude
191     * @return <code>c/sqrt(1 - s²&times;e²)</code>
192     */
193    final double msfn(final double s, final double c) {
194        return c / Math.sqrt(1.0 - (s*s) * e2);
195    }
196
197    /**
198     * Computes function (15-9) and (9-13) from Snyder.
199     * Equivalent to negative of function (7-7).
200     * @param lat the latitude
201     * @param sinlat sine of the latitude
202     * @return auxiliary value computed from <code>lat</code> and <code>sinlat</code>
203     */
204    final double tsfn(final double lat, double sinlat) {
205        sinlat *= e;
206        // NOTE: change sign to get the equivalent of Snyder (7-7).
207        return Math.tan(0.5 * (Math.PI/2 - lat)) / Math.pow((1 - sinlat) / (1 + sinlat), 0.5*e);
208    }
209}