1 | /*
|
---|
2 | * Import from fr.geo.convert package, a geographic coordinates converter.
|
---|
3 | * (http://www.i3s.unice.fr/~johan/gps/)
|
---|
4 | * License: GPL. For details, see LICENSE file.
|
---|
5 | * Copyright (C) 2002 Johan Montagnat (johan@creatis.insa-lyon.fr)
|
---|
6 | */
|
---|
7 |
|
---|
8 | package org.openstreetmap.josm.data.projection;
|
---|
9 |
|
---|
10 | /**
|
---|
11 | * the reference ellipsoids
|
---|
12 | */
|
---|
13 | class Ellipsoid {
|
---|
14 | /**
|
---|
15 | * Clarke's ellipsoid (NTF system)
|
---|
16 | */
|
---|
17 | public static final Ellipsoid clarke = new Ellipsoid(6378249.2, 6356515.0);
|
---|
18 | /**
|
---|
19 | * Hayford's ellipsoid (ED50 system)
|
---|
20 | */
|
---|
21 | public static final Ellipsoid hayford =
|
---|
22 | new Ellipsoid(6378388.0, 6356911.9461);
|
---|
23 | /**
|
---|
24 | * WGS84 ellipsoid
|
---|
25 | */
|
---|
26 | public static final Ellipsoid GRS80 = new Ellipsoid(6378137.0, 6356752.314);
|
---|
27 |
|
---|
28 | /**
|
---|
29 | * half long axis
|
---|
30 | */
|
---|
31 | public final double a;
|
---|
32 | /**
|
---|
33 | * half short axis
|
---|
34 | */
|
---|
35 | public final double b;
|
---|
36 | /**
|
---|
37 | * first eccentricity
|
---|
38 | */
|
---|
39 | public final double e;
|
---|
40 | /**
|
---|
41 | * first eccentricity squared
|
---|
42 | */
|
---|
43 | public final double e2;
|
---|
44 |
|
---|
45 | /**
|
---|
46 | * square of the second eccentricity
|
---|
47 | */
|
---|
48 | public final double eb2;
|
---|
49 |
|
---|
50 | /**
|
---|
51 | * create a new ellipsoid and precompute its parameters
|
---|
52 | *
|
---|
53 | * @param a ellipsoid long axis (in meters)
|
---|
54 | * @param b ellipsoid short axis (in meters)
|
---|
55 | */
|
---|
56 | public Ellipsoid(double a, double b) {
|
---|
57 | this.a = a;
|
---|
58 | this.b = b;
|
---|
59 | e2 = (a*a - b*b) / (a*a);
|
---|
60 | e = Math.sqrt(e2);
|
---|
61 | eb2 = e2 / (1.0 - e2);
|
---|
62 | }
|
---|
63 |
|
---|
64 | /**
|
---|
65 | * Returns the <i>radius of curvature in the prime vertical</i>
|
---|
66 | * for this reference ellipsoid at the specified latitude.
|
---|
67 | *
|
---|
68 | * @param phi The local latitude (radians).
|
---|
69 | * @return The radius of curvature in the prime vertical (meters).
|
---|
70 | */
|
---|
71 | public double verticalRadiusOfCurvature(final double phi) {
|
---|
72 | return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
|
---|
73 | }
|
---|
74 |
|
---|
75 | private static double sqr(final double x) {
|
---|
76 | return x * x;
|
---|
77 | }
|
---|
78 |
|
---|
79 | /**
|
---|
80 | * Returns the meridional arc, the true meridional distance on the
|
---|
81 | * ellipsoid from the equator to the specified latitude, in meters.
|
---|
82 | *
|
---|
83 | * @param phi The local latitude (in radians).
|
---|
84 | * @return The meridional arc (in meters).
|
---|
85 | */
|
---|
86 | public double meridionalArc(final double phi) {
|
---|
87 | final double sin2Phi = Math.sin(2.0 * phi);
|
---|
88 | final double sin4Phi = Math.sin(4.0 * phi);
|
---|
89 | final double sin6Phi = Math.sin(6.0 * phi);
|
---|
90 | final double sin8Phi = Math.sin(8.0 * phi);
|
---|
91 | // TODO . calculate 'f'
|
---|
92 | //double f = 1.0 / 298.257222101; // GRS80
|
---|
93 | double f = 1.0 / 298.257223563; // WGS84
|
---|
94 | final double n = f / (2.0 - f);
|
---|
95 | final double n2 = n * n;
|
---|
96 | final double n3 = n2 * n;
|
---|
97 | final double n4 = n3 * n;
|
---|
98 | final double n5 = n4 * n;
|
---|
99 | final double n1n2 = n - n2;
|
---|
100 | final double n2n3 = n2 - n3;
|
---|
101 | final double n3n4 = n3 - n4;
|
---|
102 | final double n4n5 = n4 - n5;
|
---|
103 | final double ap = a * (1.0 - n + (5.0 / 4.0) * (n2n3) + (81.0 / 64.0) * (n4n5));
|
---|
104 | final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * (n3n4) + (55.0 / 64.0) * n5);
|
---|
105 | final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * (n4n5));
|
---|
106 | final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
|
---|
107 | final double ep = (315.0 / 512.0) * a * (n4n5);
|
---|
108 | return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
|
---|
109 | }
|
---|
110 |
|
---|
111 | /**
|
---|
112 | * Returns the <i>radius of curvature in the meridian<i>
|
---|
113 | * for this reference ellipsoid at the specified latitude.
|
---|
114 | *
|
---|
115 | * @param phi The local latitude (in radians).
|
---|
116 | * @return The radius of curvature in the meridian (in meters).
|
---|
117 | */
|
---|
118 | public double meridionalRadiusOfCurvature(final double phi) {
|
---|
119 | return verticalRadiusOfCurvature(phi)
|
---|
120 | / (1.0 + eb2 * sqr(Math.cos(phi)));
|
---|
121 | }
|
---|
122 |
|
---|
123 | /**
|
---|
124 | * Returns isometric latitude of phi on given first eccentricity (e)
|
---|
125 | * @param phi The local latitude (radians).
|
---|
126 | * @return isometric latitude of phi on first eccentricity (e)
|
---|
127 | */
|
---|
128 | public double latitudeIsometric(double phi, double e) {
|
---|
129 | double v1 = 1-e*Math.sin(phi);
|
---|
130 | double v2 = 1+e*Math.sin(phi);
|
---|
131 | return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
|
---|
132 | }
|
---|
133 |
|
---|
134 | /**
|
---|
135 | * Returns isometric latitude of phi on first eccentricity (e)
|
---|
136 | * @param phi The local latitude (radians).
|
---|
137 | * @return isometric latitude of phi on first eccentricity (e)
|
---|
138 | */
|
---|
139 | public double latitudeIsometric(double phi) {
|
---|
140 | double v1 = 1-e*Math.sin(phi);
|
---|
141 | double v2 = 1+e*Math.sin(phi);
|
---|
142 | return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2,e/2));
|
---|
143 | }
|
---|
144 |
|
---|
145 | /*
|
---|
146 | * Returns geographic latitude of isometric latitude of first eccentricity (e)
|
---|
147 | * and epsilon precision
|
---|
148 | */
|
---|
149 | public double latitude(double latIso, double e, double epsilon) {
|
---|
150 | double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
|
---|
151 | double lati = lat0;
|
---|
152 | double lati1 = 1.0; // random value to start the iterative processus
|
---|
153 | while(Math.abs(lati1-lati)>=epsilon) {
|
---|
154 | lati = lati1;
|
---|
155 | double v1 = 1+e*Math.sin(lati);
|
---|
156 | double v2 = 1-e*Math.sin(lati);
|
---|
157 | lati1 = 2*Math.atan(Math.pow(v1/v2,e/2)*Math.exp(latIso))-Math.PI/2;
|
---|
158 | }
|
---|
159 | return lati1;
|
---|
160 | }
|
---|
161 | }
|
---|
162 |
|
---|
163 |
|
---|
164 |
|
---|