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