001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.data.projection.proj; 003 004import org.openstreetmap.josm.data.projection.ProjectionConfigurationException; 005 006/** 007 * Abstract base class providing utilities for implementations of the Proj 008 * interface. 009 * 010 * This class has been derived from the implementation of the Geotools project; 011 * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection 012 * at the time of migration. 013 * <p> 014 * 015 * @author André Gosselin 016 * @author Martin Desruisseaux (PMO, IRD) 017 * @author Rueben Schulz 018*/ 019public abstract class AbstractProj implements Proj { 020 021 /** 022 * Maximum number of iterations for iterative computations. 023 */ 024 private static final int MAXIMUM_ITERATIONS = 15; 025 026 /** 027 * Relative iteration precision used in the <code>mlfn</code> method 028 */ 029 private static final double MLFN_TOL = 1E-11; 030 031 /** 032 * Constants used to calculate {@link #en0}, {@link #en1}, 033 * {@link #en2}, {@link #en3}, {@link #en4}. 034 */ 035 private static final double C00 = 1.0, 036 C02 = 0.25, 037 C04 = 0.046875, 038 C06 = 0.01953125, 039 C08 = 0.01068115234375, 040 C22 = 0.75, 041 C44 = 0.46875, 042 C46 = 0.01302083333333333333, 043 C48 = 0.00712076822916666666, 044 C66 = 0.36458333333333333333, 045 C68 = 0.00569661458333333333, 046 C88 = 0.3076171875; 047 048 /** 049 * Constant needed for the <code>mlfn</code> method. 050 * Setup at construction time. 051 */ 052 protected double en0, en1, en2, en3, en4; 053 054 /** 055 * The square of excentricity: e² = (a²-b²)/a² where 056 * <var>e</var> is the excentricity, 057 * <var>a</var> is the semi major axis length and 058 * <var>b</var> is the semi minor axis length. 059 */ 060 protected double e2; 061 062 @Override 063 public void initialize(ProjParameters params) throws ProjectionConfigurationException { 064 e2 = params.ellps.e2; 065 // Compute constants for the mlfn 066 double t; 067 en0 = C00 - e2 * (C02 + e2 * 068 (C04 + e2 * (C06 + e2 * C08))); 069 en1 = e2 * (C22 - e2 * 070 (C04 + e2 * (C06 + e2 * C08))); 071 en2 = (t = e2 * e2) * 072 (C44 - e2 * (C46 + e2 * C48)); 073 en3 = (t *= e2) * (C66 - e2 * C68); 074 en4 = t * e2 * C88; 075 } 076 077 /** 078 * Calculates the meridian distance. This is the distance along the central 079 * meridian from the equator to {@code phi}. Accurate to < 1e-5 meters 080 * when used in conjuction with typical major axis values. 081 * 082 * @param phi latitude to calculate meridian distance for. 083 * @param sphi sin(phi). 084 * @param cphi cos(phi). 085 * @return meridian distance for the given latitude. 086 */ 087 protected final double mlfn(final double phi, double sphi, double cphi) { 088 cphi *= sphi; 089 sphi *= sphi; 090 return en0 * phi - cphi * 091 (en1 + sphi * 092 (en2 + sphi * 093 (en3 + sphi * 094 (en4)))); 095 } 096 097 /** 098 * Calculates the latitude ({@code phi}) from a meridian distance. 099 * Determines phi to TOL (1e-11) radians, about 1e-6 seconds. 100 * 101 * @param arg meridian distance to calulate latitude for. 102 * @return the latitude of the meridian distance. 103 * @throws RuntimeException if the itteration does not converge. 104 */ 105 protected final double inv_mlfn(double arg) { 106 double s, t, phi, k = 1.0/(1.0 - e2); 107 int i; 108 phi = arg; 109 for (i = MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations 110 if (--i < 0) { 111 throw new RuntimeException("Too many iterations"); 112 } 113 s = Math.sin(phi); 114 t = 1.0 - e2 * s * s; 115 t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k; 116 phi -= t; 117 if (Math.abs(t) < MLFN_TOL) { 118 return phi; 119 } 120 } 121 } 122 123 public static double normalizeLon(double lon) { 124 if (lon >= -Math.PI && lon <= Math.PI) 125 return lon; 126 else { 127 lon = lon % (2 * Math.PI); 128 if (lon > Math.PI) { 129 return lon - 2 * Math.PI; 130 } else if (lon < -Math.PI) { 131 return lon + 2 * Math.PI; 132 } 133 return lon; 134 } 135 } 136}