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