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 < 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²×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²×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}