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