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 | import org.openstreetmap.josm.tools.CheckParameterUtil;
|
---|
7 |
|
---|
8 | /**
|
---|
9 | * Abstract base class providing utilities for implementations of the Proj
|
---|
10 | * interface.
|
---|
11 | *
|
---|
12 | * This class has been derived from the implementation of the Geotools project;
|
---|
13 | * git 8cbf52d, org.geotools.referencing.operation.projection.MapProjection
|
---|
14 | * at the time of migration.
|
---|
15 | * <p>
|
---|
16 | *
|
---|
17 | * @author André Gosselin
|
---|
18 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
19 | * @author Rueben Schulz
|
---|
20 | */
|
---|
21 | public abstract class AbstractProj implements Proj {
|
---|
22 |
|
---|
23 | /**
|
---|
24 | * Maximum number of iterations for iterative computations.
|
---|
25 | */
|
---|
26 | private static final int MAXIMUM_ITERATIONS = 15;
|
---|
27 |
|
---|
28 | /**
|
---|
29 | * Difference allowed in iterative computations.
|
---|
30 | */
|
---|
31 | private static final double ITERATION_TOLERANCE = 1E-10;
|
---|
32 |
|
---|
33 | /**
|
---|
34 | * Relative iteration precision used in the <code>mlfn</code> method
|
---|
35 | */
|
---|
36 | private static final double MLFN_TOL = 1E-11;
|
---|
37 |
|
---|
38 | /**
|
---|
39 | * Constants used to calculate {@link #en0}, {@link #en1},
|
---|
40 | * {@link #en2}, {@link #en3}, {@link #en4}.
|
---|
41 | */
|
---|
42 | private static final double C00 = 1.0;
|
---|
43 | private static final double C02 = 0.25;
|
---|
44 | private static final double C04 = 0.046875;
|
---|
45 | private static final double C06 = 0.01953125;
|
---|
46 | private static final double C08 = 0.01068115234375;
|
---|
47 | private static final double C22 = 0.75;
|
---|
48 | private static final double C44 = 0.46875;
|
---|
49 | private static final double C46 = 0.01302083333333333333;
|
---|
50 | private static final double C48 = 0.00712076822916666666;
|
---|
51 | private static final double C66 = 0.36458333333333333333;
|
---|
52 | private static final double C68 = 0.00569661458333333333;
|
---|
53 | private static final double C88 = 0.3076171875;
|
---|
54 |
|
---|
55 | /**
|
---|
56 | * Constant needed for the <code>mlfn</code> method.
|
---|
57 | * Setup at construction time.
|
---|
58 | */
|
---|
59 | protected double en0, en1, en2, en3, en4;
|
---|
60 |
|
---|
61 | /**
|
---|
62 | * Ellipsoid excentricity, equals to <code>sqrt({@link #e2 excentricity squared})</code>.
|
---|
63 | * Value 0 means that the ellipsoid is spherical.
|
---|
64 | *
|
---|
65 | * @see #e2
|
---|
66 | */
|
---|
67 | protected double e;
|
---|
68 |
|
---|
69 | /**
|
---|
70 | * The square of excentricity: e² = (a²-b²)/a² where
|
---|
71 | * <var>e</var> is the excentricity,
|
---|
72 | * <var>a</var> is the semi major axis length and
|
---|
73 | * <var>b</var> is the semi minor axis length.
|
---|
74 | *
|
---|
75 | * @see #e
|
---|
76 | */
|
---|
77 | protected double e2;
|
---|
78 |
|
---|
79 | /**
|
---|
80 | * is ellipsoid spherical?
|
---|
81 | * @see Ellipsoid#spherical
|
---|
82 | */
|
---|
83 | protected boolean spherical;
|
---|
84 |
|
---|
85 | @Override
|
---|
86 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
87 | CheckParameterUtil.ensureParameterNotNull(params, "params");
|
---|
88 | CheckParameterUtil.ensureParameterNotNull(params.ellps, "params.ellps");
|
---|
89 | e2 = params.ellps.e2;
|
---|
90 | e = params.ellps.e;
|
---|
91 | spherical = params.ellps.spherical;
|
---|
92 | // Compute constants for the mlfn
|
---|
93 | double t;
|
---|
94 | // CHECKSTYLE.OFF: SingleSpaceSeparator
|
---|
95 | en0 = C00 - e2 * (C02 + e2 *
|
---|
96 | (C04 + e2 * (C06 + e2 * C08)));
|
---|
97 | en1 = e2 * (C22 - e2 *
|
---|
98 | (C04 + e2 * (C06 + e2 * C08)));
|
---|
99 | 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 conjuction 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 calulate 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 | // Iteratively solve equation (7-9) from Snyder.
|
---|
158 | final double cphi2(final double ts) {
|
---|
159 | final double eccnth = 0.5 * e;
|
---|
160 | double phi = (Math.PI/2) - 2.0 * Math.atan(ts);
|
---|
161 | for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
|
---|
162 | final double con = e * Math.sin(phi);
|
---|
163 | final double dphi = (Math.PI/2) - 2.0*Math.atan(ts * Math.pow((1-con)/(1+con), eccnth)) - phi;
|
---|
164 | phi += dphi;
|
---|
165 | if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
|
---|
166 | return phi;
|
---|
167 | }
|
---|
168 | }
|
---|
169 | throw new IllegalStateException("no convergence for ts="+ts);
|
---|
170 | }
|
---|
171 |
|
---|
172 | /**
|
---|
173 | * Computes function <code>f(s,c,e²) = c/sqrt(1 - s²×e²)</code> needed for the true scale
|
---|
174 | * latitude (Snyder 14-15), where <var>s</var> and <var>c</var> are the sine and cosine of
|
---|
175 | * the true scale latitude, and <var>e²</var> is the {@linkplain #e2 eccentricity squared}.
|
---|
176 | * @param s sine of the true scale latitude
|
---|
177 | * @param c cosine of the true scale latitude
|
---|
178 | * @return <code>c/sqrt(1 - s²×e²)</code>
|
---|
179 | */
|
---|
180 | final double msfn(final double s, final double c) {
|
---|
181 | return c / Math.sqrt(1.0 - (s*s) * e2);
|
---|
182 | }
|
---|
183 |
|
---|
184 | /**
|
---|
185 | * Computes function (15-9) and (9-13) from Snyder.
|
---|
186 | * Equivalent to negative of function (7-7).
|
---|
187 | * @param lat the latitude
|
---|
188 | * @param sinlat sine of the latitude
|
---|
189 | * @return auxiliary value computed from <code>lat</code> and <code>sinlat</code>
|
---|
190 | */
|
---|
191 | final double tsfn(final double lat, double sinlat) {
|
---|
192 | sinlat *= e;
|
---|
193 | // NOTE: change sign to get the equivalent of Snyder (7-7).
|
---|
194 | return Math.tan(0.5 * (Math.PI/2 - lat)) / Math.pow((1 - sinlat) / (1 + sinlat), 0.5*e);
|
---|
195 | }
|
---|
196 | }
|
---|