1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import 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 | */
|
---|
19 | public 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 | * Relative iteration precision used in the <code>mlfn</code> method
|
---|
28 | */
|
---|
29 | private static final double MLFN_TOL = 1E-11;
|
---|
30 |
|
---|
31 | /**
|
---|
32 | * Constants used to calculate {@link #en0}, {@link #en1},
|
---|
33 | * {@link #en2}, {@link #en3}, {@link #en4}.
|
---|
34 | */
|
---|
35 | private static final double C00 = 1.0,
|
---|
36 | C02 = 0.25,
|
---|
37 | C04 = 0.046875,
|
---|
38 | C06 = 0.01953125,
|
---|
39 | C08 = 0.01068115234375,
|
---|
40 | C22 = 0.75,
|
---|
41 | C44 = 0.46875,
|
---|
42 | C46 = 0.01302083333333333333,
|
---|
43 | C48 = 0.00712076822916666666,
|
---|
44 | C66 = 0.36458333333333333333,
|
---|
45 | C68 = 0.00569661458333333333,
|
---|
46 | C88 = 0.3076171875;
|
---|
47 |
|
---|
48 | /**
|
---|
49 | * Constant needed for the <code>mlfn</code> method.
|
---|
50 | * Setup at construction time.
|
---|
51 | */
|
---|
52 | protected double en0, en1, en2, en3, en4;
|
---|
53 |
|
---|
54 | /**
|
---|
55 | * The square of excentricity: e² = (a²-b²)/a² where
|
---|
56 | * <var>e</var> is the excentricity,
|
---|
57 | * <var>a</var> is the semi major axis length and
|
---|
58 | * <var>b</var> is the semi minor axis length.
|
---|
59 | */
|
---|
60 | protected double e2;
|
---|
61 |
|
---|
62 | @Override
|
---|
63 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
64 | e2 = params.ellps.e2;
|
---|
65 | // Compute constants for the mlfn
|
---|
66 | double t;
|
---|
67 | en0 = C00 - e2 * (C02 + e2 *
|
---|
68 | (C04 + e2 * (C06 + e2 * C08)));
|
---|
69 | en1 = e2 * (C22 - e2 *
|
---|
70 | (C04 + e2 * (C06 + e2 * C08)));
|
---|
71 | en2 = (t = e2 * e2) *
|
---|
72 | (C44 - e2 * (C46 + e2 * C48));
|
---|
73 | en3 = (t *= e2) * (C66 - e2 * C68);
|
---|
74 | en4 = t * e2 * C88;
|
---|
75 | }
|
---|
76 |
|
---|
77 | /**
|
---|
78 | * Calculates the meridian distance. This is the distance along the central
|
---|
79 | * meridian from the equator to {@code phi}. Accurate to < 1e-5 meters
|
---|
80 | * when used in conjuction with typical major axis values.
|
---|
81 | *
|
---|
82 | * @param phi latitude to calculate meridian distance for.
|
---|
83 | * @param sphi sin(phi).
|
---|
84 | * @param cphi cos(phi).
|
---|
85 | * @return meridian distance for the given latitude.
|
---|
86 | */
|
---|
87 | protected final double mlfn(final double phi, double sphi, double cphi) {
|
---|
88 | cphi *= sphi;
|
---|
89 | sphi *= sphi;
|
---|
90 | return en0 * phi - cphi *
|
---|
91 | (en1 + sphi *
|
---|
92 | (en2 + sphi *
|
---|
93 | (en3 + sphi *
|
---|
94 | (en4))));
|
---|
95 | }
|
---|
96 |
|
---|
97 | /**
|
---|
98 | * Calculates the latitude ({@code phi}) from a meridian distance.
|
---|
99 | * 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 | }
|
---|