1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import org.openstreetmap.josm.data.projection.Ellipsoid;
|
---|
5 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
6 |
|
---|
7 | /**
|
---|
8 | * Abstract base class providing utilities for implementations of the Proj
|
---|
9 | * interface.
|
---|
10 | *
|
---|
11 | * This class has been derived from the implementation of the Geotools project;
|
---|
12 | * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection
|
---|
13 | * at the time of migration.
|
---|
14 | * <p>
|
---|
15 | *
|
---|
16 | * @author André Gosselin
|
---|
17 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
18 | * @author Rueben Schulz
|
---|
19 | */
|
---|
20 | public abstract class AbstractProj implements Proj {
|
---|
21 |
|
---|
22 | /**
|
---|
23 | * Maximum number of iterations for iterative computations.
|
---|
24 | */
|
---|
25 | private static final int MAXIMUM_ITERATIONS = 15;
|
---|
26 |
|
---|
27 | /**
|
---|
28 | * Difference allowed in iterative computations.
|
---|
29 | */
|
---|
30 | private static final double ITERATION_TOLERANCE = 1E-10;
|
---|
31 |
|
---|
32 | /**
|
---|
33 | * Relative iteration precision used in the <code>mlfn</code> method
|
---|
34 | */
|
---|
35 | private static final double MLFN_TOL = 1E-11;
|
---|
36 |
|
---|
37 | /**
|
---|
38 | * Constants used to calculate {@link #en0}, {@link #en1},
|
---|
39 | * {@link #en2}, {@link #en3}, {@link #en4}.
|
---|
40 | */
|
---|
41 | private static final double C00 = 1.0;
|
---|
42 | private static final double C02 = 0.25;
|
---|
43 | private static final double C04 = 0.046875;
|
---|
44 | private static final double C06 = 0.01953125;
|
---|
45 | private static final double C08 = 0.01068115234375;
|
---|
46 | private static final double C22 = 0.75;
|
---|
47 | private static final double C44 = 0.46875;
|
---|
48 | private static final double C46 = 0.01302083333333333333;
|
---|
49 | private static final double C48 = 0.00712076822916666666;
|
---|
50 | private static final double C66 = 0.36458333333333333333;
|
---|
51 | private static final double C68 = 0.00569661458333333333;
|
---|
52 | private static final double C88 = 0.3076171875;
|
---|
53 |
|
---|
54 | /**
|
---|
55 | * Constant needed for the <code>mlfn</code> method.
|
---|
56 | * Setup at construction time.
|
---|
57 | */
|
---|
58 | protected double en0, en1, en2, en3, en4;
|
---|
59 |
|
---|
60 | /**
|
---|
61 | * Ellipsoid excentricity, equals to <code>sqrt({@link #e2 excentricity squared})</code>.
|
---|
62 | * Value 0 means that the ellipsoid is spherical.
|
---|
63 | *
|
---|
64 | * @see #e2
|
---|
65 | */
|
---|
66 | protected double e;
|
---|
67 |
|
---|
68 | /**
|
---|
69 | * The square of excentricity: e² = (a²-b²)/a² where
|
---|
70 | * <var>e</var> is the excentricity,
|
---|
71 | * <var>a</var> is the semi major axis length and
|
---|
72 | * <var>b</var> is the semi minor axis length.
|
---|
73 | *
|
---|
74 | * @see #e
|
---|
75 | */
|
---|
76 | protected double e2;
|
---|
77 |
|
---|
78 | /**
|
---|
79 | * is ellipsoid spherical?
|
---|
80 | * @see Ellipsoid#spherical
|
---|
81 | */
|
---|
82 | protected boolean spherical;
|
---|
83 |
|
---|
84 | @Override
|
---|
85 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
86 | e2 = params.ellps.e2;
|
---|
87 | e = params.ellps.e;
|
---|
88 | spherical = params.ellps.spherical;
|
---|
89 | // Compute constants for the mlfn
|
---|
90 | double t;
|
---|
91 | en0 = C00 - e2 * (C02 + e2 *
|
---|
92 | (C04 + e2 * (C06 + e2 * C08)));
|
---|
93 | en1 = e2 * (C22 - e2 *
|
---|
94 | (C04 + e2 * (C06 + e2 * C08)));
|
---|
95 | en2 = (t = e2 * e2) *
|
---|
96 | (C44 - e2 * (C46 + e2 * C48));
|
---|
97 | en3 = (t *= e2) * (C66 - e2 * C68);
|
---|
98 | en4 = t * e2 * C88;
|
---|
99 | }
|
---|
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 IllegalStateException("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 IllegalStateException("no convergence for ts="+ts);
|
---|
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 | }
|
---|