1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28 package edu.sc.seis.TauP;
29
30
31 /*** Utility class for spherical coordinate (lat-lon) transformations. Given
32 * lat, lon, lat, lon you can find the distance or azimuth and
33 * given lat, lon, distance, azimuth you can find the lat lon of the resultant
34 * point. Just uses spherical relations, no ellpticity correction is applied.
35 *
36 * See Appendix A of "Seismology and Plate Tectonics" by David Gubbins
37 * Cambridge University Press, 1990
38 *
39 * and
40 * Chapter 3 of "Plate Tectonics: How it Works"
41 * by Allan Cox and Robert Brian Hart
42 * Blackwell Scientific Publications, 1986
43 *
44 * @author H. Philip Crotwell
45 * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
46
47
48
49 */
50 public class SphericalCoords {
51
52 protected static final double dtor = Math.PI/180.0;
53 protected static final double rtod = 180.0/Math.PI;
54
55 /*** Calculates angular distance between two lat lon pairs. */
56 public static double distance(double latA, double lonA, double latB, double lonB) {
57 return rtod * Math.acos(Math.sin(latA*dtor)*Math.sin(latB*dtor) +
58 Math.cos(latA*dtor)*Math.cos(latB*dtor)*Math.cos((lonB-lonA)*dtor));
59 }
60
61 /*** Calculates azimuth between two lat lon pairs. */
62 public static double azimuth(double latA, double lonA,
63 double latB, double lonB) {
64 double cosAzimuth = (Math.cos(latA*dtor)*Math.sin(latB*dtor) -
65 Math.sin(latA*dtor)*Math.cos(latB*dtor)*Math.cos((lonB-lonA)*dtor)) /
66 Math.sin(distance(latA, lonA, latB, lonB)*dtor);
67
68 double sinAzimuth = Math.cos(latB*dtor)*Math.sin((lonB-lonA)*dtor) /
69 Math.sin(distance(latA, lonA, latB, lonB)*dtor);
70 return rtod*Math.atan2(sinAzimuth, cosAzimuth);
71 }
72
73 /*** Find the rotation pole required to rotate the first lat lon pair to
74 * the second. Just does a cross product.
75 * @returns a 3 element double array with the X, Y and Z components of
76 * the pole.
77 */
78 public static double[] rotationPole(double latA, double lonA,
79 double latB, double lonB) {
80 double[] pointA = new double[3];
81 double[] pointB = new double[3];
82 double[] pole = new double[3];
83
84 double dToR = Math.PI/180.0;
85
86 pointA[0] = Math.cos(latA*dToR) * Math.cos(lonA*dToR);
87 pointA[1] = Math.cos(latA*dToR) * Math.sin(lonA*dToR);
88 pointA[2] = Math.sin(latA*dToR);
89
90 pointB[0] = Math.cos(latB*dToR) * Math.cos(lonB*dToR);
91 pointB[1] = Math.cos(latB*dToR) * Math.sin(lonB*dToR);
92 pointB[2] = Math.sin(latB*dToR);
93
94 pole[0] = pointA[1]*pointB[2] - pointA[2]*pointB[1];
95 pole[1] = pointA[2]*pointB[0] - pointA[0]*pointB[2];
96 pole[2] = pointA[0]*pointB[1] - pointA[1]*pointB[0];
97
98 return pole;
99 }
100
101 /*** rotates a point about a pole by an angle.
102 * @param pole is a 3 element double array with X, Y and Z components of
103 * the pole.
104 * @returns [lat, lon] in array. */
105 public static double[] rotate(double latA, double lonA,
106 double[] pole, double angleDeg) {
107 double[][] R = new double[3][3];
108 double[] point = new double[3];
109 double[] newPoint = new double[3];
110
111 double rToDeg = 180.0/Math.PI;
112 double angle = angleDeg / rToDeg;
113
114 R[0][0] = pole[0]*pole[0]*(1-Math.cos(angle)) + Math.cos(angle);
115 R[0][1] = pole[0]*pole[1]*(1-Math.cos(angle)) - pole[2]*Math.sin(angle);
116 R[0][2] = pole[0]*pole[2]*(1-Math.cos(angle)) + pole[1]*Math.sin(angle);
117
118 R[1][0] = pole[1]*pole[0]*(1-Math.cos(angle)) + pole[2]*Math.sin(angle);
119 R[1][1] = pole[1]*pole[1]*(1-Math.cos(angle)) + Math.cos(angle);
120 R[1][2] = pole[1]*pole[2]*(1-Math.cos(angle)) - pole[0]*Math.sin(angle);
121
122 R[2][0] = pole[2]*pole[0]*(1-Math.cos(angle)) - pole[1]*Math.sin(angle);
123 R[2][1] = pole[2]*pole[1]*(1-Math.cos(angle)) + pole[0]*Math.sin(angle);
124 R[2][2] = pole[2]*pole[2]*(1-Math.cos(angle)) + Math.cos(angle);
125
126 point[0] = Math.cos(latA/rToDeg) * Math.cos(lonA/rToDeg);
127 point[1] = Math.cos(latA/rToDeg) * Math.sin(lonA/rToDeg);
128 point[2] = Math.sin(latA/rToDeg);
129
130 newPoint[0] = R[0][0]*point[0] +
131 R[0][1]*point[1] +
132 R[0][2]*point[2];
133
134 newPoint[1] = R[1][0]*point[0] +
135 R[1][1]*point[1] +
136 R[1][2]*point[2];
137
138 newPoint[2] = R[2][0]*point[0] +
139 R[2][1]*point[1] +
140 R[2][2]*point[2];
141
142 double newLat = Math.asin(newPoint[2])*180.0/Math.PI;
143 double newLon = Math.atan2(newPoint[1],newPoint[0])*180.0/Math.PI;
144
145 newPoint = new double[2];
146 newPoint[0] = newLat;
147 newPoint[1] = newLon;
148
149 return newPoint;
150 }
151
152 /*** Calculates the latitude of a point a given distance along a given
153 * azimuth from a starting lat lon. */
154 public static double latFor(double latA, double lonA,
155 double distance, double azimuth) {
156 return rtod*Math.asin(
157 Math.cos(azimuth*dtor)*Math.sin(distance*dtor)*Math.cos(latA*dtor) +
158 Math.cos(distance*dtor)* Math.sin(latA*dtor));
159 }
160
161 /*** Calculates the longitude of a point a given distance along a given
162 * azimuth from a starting lat lon. */
163 public static double lonFor(double latA, double lonA,
164 double distance, double azimuth) {
165 double tempLat = latFor(latA, lonA, distance, azimuth);
166
167 double sinLon = Math.sin(azimuth*dtor)*Math.sin(distance*dtor) /
168 Math.cos(tempLat*dtor);
169 double cosLon = (Math.cos(distance*dtor)-
170 Math.sin(latA*dtor)*Math.sin(tempLat*dtor)) /
171 (Math.cos(latA*dtor)*Math.cos(tempLat*dtor));
172
173
174 double lon = lonA+rtod*Math.atan2(sinLon, cosLon);
175 if (lon <= -180.0) lon += 360.0;
176 if (lon > 180.0) lon -= 360.0;
177 return lon;
178 }
179
180 public static void main(String args[]) {
181
182 System.out.println(distance(0,0,0,45)+" "+
183 azimuth(0,0,0,45)+" "+azimuth(0,45,0,0));
184
185 System.out.println(latFor(0,0,45,90)+" "+
186 lonFor(0,0,45,90));
187
188 System.out.println("(35,42,36,43) "+distance(35,42,36,43)+" "+
189 azimuth(35,42,36,43)+" "+azimuth(36,43,35,42));
190
191 System.out.println(
192 latFor(35,42,distance(35,42,36,43), azimuth(35,42,36,43))+" "+
193 lonFor(35,42,distance(35,42,36,43), azimuth(35,42,36,43))
194 );
195 }
196 }