source: josm/trunk/src/org/openstreetmap/josm/data/projection/proj/AbstractProj.java@ 9558

Last change on this file since 9558 was 9558, checked in by bastiK, 8 years ago

always normalize longitude before projection and after inverse projection (see #12186)
  • Property svn:eol-style set to native
File size: 6.2 KB
Line 
1// License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.data.projection.proj;
3
4import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
5
6/**
7 * Abstract base class providing utilities for implementations of the Proj
8 * interface.
9 *
10 * This class has been derived from the implementation of the Geotools project;
11 * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection
12 * at the time of migration.
13 * <p>
14 *
15 * @author André Gosselin
16 * @author Martin Desruisseaux (PMO, IRD)
17 * @author Rueben Schulz
18*/
19public abstract class AbstractProj implements Proj {
20
21 /**
22 * Maximum number of iterations for iterative computations.
23 */
24 private static final int MAXIMUM_ITERATIONS = 15;
25
26 /**
27 * Difference allowed in iterative computations.
28 */
29 private static final double ITERATION_TOLERANCE = 1E-10;
30
31 /**
32 * Relative iteration precision used in the <code>mlfn</code> method
33 */
34 private static final double MLFN_TOL = 1E-11;
35
36 /**
37 * Constants used to calculate {@link #en0}, {@link #en1},
38 * {@link #en2}, {@link #en3}, {@link #en4}.
39 */
40 private static final double C00 = 1.0,
41 C02 = 0.25,
42 C04 = 0.046875,
43 C06 = 0.01953125,
44 C08 = 0.01068115234375,
45 C22 = 0.75,
46 C44 = 0.46875,
47 C46 = 0.01302083333333333333,
48 C48 = 0.00712076822916666666,
49 C66 = 0.36458333333333333333,
50 C68 = 0.00569661458333333333,
51 C88 = 0.3076171875;
52
53 /**
54 * Constant needed for the <code>mlfn</code> method.
55 * Setup at construction time.
56 */
57 protected double en0, en1, en2, en3, en4;
58
59 /**
60 * Ellipsoid excentricity, equals to <code>sqrt({@link #e2 excentricity squared})</code>.
61 * Value 0 means that the ellipsoid is spherical.
62 *
63 * @see #e2
64 */
65 protected double e;
66
67 /**
68 * The square of excentricity: e² = (a²-b²)/a² where
69 * <var>e</var> is the excentricity,
70 * <var>a</var> is the semi major axis length and
71 * <var>b</var> is the semi minor axis length.
72 *
73 * @see #e
74 */
75 protected double e2;
76
77 @Override
78 public void initialize(ProjParameters params) throws ProjectionConfigurationException {
79 e2 = params.ellps.e2;
80 e = params.ellps.e;
81 // Compute constants for the mlfn
82 double t;
83 en0 = C00 - e2 * (C02 + e2 *
84 (C04 + e2 * (C06 + e2 * C08)));
85 en1 = e2 * (C22 - e2 *
86 (C04 + e2 * (C06 + e2 * C08)));
87 en2 = (t = e2 * e2) *
88 (C44 - e2 * (C46 + e2 * C48));
89 en3 = (t *= e2) * (C66 - e2 * C68);
90 en4 = t * e2 * C88;
91 }
92
93 /**
94 * Calculates the meridian distance. This is the distance along the central
95 * meridian from the equator to {@code phi}. Accurate to &lt; 1e-5 meters
96 * when used in conjuction with typical major axis values.
97 *
98 * @param phi latitude to calculate meridian distance for.
99 * @param sphi sin(phi).
100 * @param cphi cos(phi).
101 * @return meridian distance for the given latitude.
102 */
103 protected final double mlfn(final double phi, double sphi, double cphi) {
104 cphi *= sphi;
105 sphi *= sphi;
106 return en0 * phi - cphi *
107 (en1 + sphi *
108 (en2 + sphi *
109 (en3 + sphi *
110 (en4))));
111 }
112
113 /**
114 * Calculates the latitude ({@code phi}) from a meridian distance.
115 * Determines phi to TOL (1e-11) radians, about 1e-6 seconds.
116 *
117 * @param arg meridian distance to calulate latitude for.
118 * @return the latitude of the meridian distance.
119 * @throws RuntimeException if the itteration does not converge.
120 */
121 protected final double inv_mlfn(double arg) {
122 double s, t, phi, k = 1.0/(1.0 - e2);
123 int i;
124 phi = arg;
125 for (i = MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations
126 if (--i < 0) {
127 throw new RuntimeException("Too many iterations");
128 }
129 s = Math.sin(phi);
130 t = 1.0 - e2 * s * s;
131 t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k;
132 phi -= t;
133 if (Math.abs(t) < MLFN_TOL) {
134 return phi;
135 }
136 }
137 }
138
139 // Iteratively solve equation (7-9) from Snyder.
140 final double cphi2(final double ts) {
141 final double eccnth = 0.5 * e;
142 double phi = (Math.PI/2) - 2.0 * Math.atan(ts);
143 for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
144 final double con = e * Math.sin(phi);
145 final double dphi = (Math.PI/2) - 2.0*Math.atan(ts * Math.pow((1-con)/(1+con), eccnth)) - phi;
146 phi += dphi;
147 if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
148 return phi;
149 }
150 }
151 throw new RuntimeException("no convergence");
152 }
153
154 /**
155 * Computes function <code>f(s,c,e²) = c/sqrt(1 - s²&times;e²)</code> needed for the true scale
156 * latitude (Snyder 14-15), where <var>s</var> and <var>c</var> are the sine and cosine of
157 * the true scale latitude, and <var>e²</var> is the {@linkplain #e2 eccentricity squared}.
158 * @param s sine of the true scale latitude
159 * @param c cosine of the true scale latitude
160 * @return <code>c/sqrt(1 - s²&times;e²)</code>
161 */
162 final double msfn(final double s, final double c) {
163 return c / Math.sqrt(1.0 - (s*s) * e2);
164 }
165
166 /**
167 * Computes function (15-9) and (9-13) from Snyder.
168 * Equivalent to negative of function (7-7).
169 * @param lat the latitude
170 * @param sinlat sine of the latitude
171 * @return auxiliary value computed from <code>lat</code> and <code>sinlat</code>
172 */
173 final double tsfn(final double lat, double sinlat) {
174 sinlat *= e;
175 // NOTE: change sign to get the equivalent of Snyder (7-7).
176 return Math.tan(0.5 * (Math.PI/2 - lat)) / Math.pow((1 - sinlat) / (1 + sinlat), 0.5*e);
177 }
178}
Note: See TracBrowser for help on using the repository browser.