001 /*
002 * Zmanim Java API
003 * Copyright (C) 2004-2008 Eliyahu Hershfeld
004 *
005 * This program is free software; you can redistribute it and/or modify it under the terms of the
006 * GNU General Public License as published by the Free Software Foundation; either version 2 of the
007 * License, or (at your option) any later version.
008 *
009 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
010 * even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
011 * General Public License for more details.
012 *
013 * You should have received a copy of the GNU General Public License along with this program; if
014 * not, write to the Free Software Foundation, Inc. 59 Temple Place - Suite 330, Boston, MA
015 * 02111-1307, USA or connect to: http://www.fsf.org/copyleft/gpl.html
016 */
017 package net.sourceforge.zmanim.util;
018
019 import java.util.TimeZone;
020
021 /**
022 * A class for various location calculations
023 * Most of the code in this class is ported from <a href="http://www.movable-type.co.uk/">Chris Veness'</a>
024 * <a href="http://www.fsf.org/licensing/licenses/lgpl.html">LGPL</a> Javascript Implementation
025 *
026 * @author © Eliyahu Hershfeld 2008
027 * @version 0.1
028 */
029 public class GeoLocationUtils {
030 private static int DISTANCE = 0;
031 private static int INITIAL_BEARING = 1;
032 private static int FINAL_BEARING = 2;
033
034 /**
035 * Calculate the initial <a
036 * href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing
037 * between this Object and a second Object passed to this method using <a
038 * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
039 * inverse formula See T Vincenty, "<a
040 * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
041 * Solutions of Geodesics on the Ellipsoid with application of nested
042 * equations</a>", Survey Review, vol XXII no 176, 1975.
043 *
044 * @param location
045 * the destination location
046 */
047 public static double getGeodesicInitialBearing(GeoLocation location, GeoLocation destination) {
048 return vincentyFormula(location, destination, INITIAL_BEARING);
049 }
050
051 /**
052 * Calculate the final <a
053 * href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing
054 * between this Object and a second Object passed to this method using <a
055 * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
056 * inverse formula See T Vincenty, "<a
057 * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
058 * Solutions of Geodesics on the Ellipsoid with application of nested
059 * equations</a>", Survey Review, vol XXII no 176, 1975.
060 *
061 * @param location
062 * the destination location
063 */
064 public static double getGeodesicFinalBearing(GeoLocation location, GeoLocation destination) {
065 return vincentyFormula(location, destination, FINAL_BEARING);
066 }
067
068 /**
069 * Calculate <a
070 * href="http://en.wikipedia.org/wiki/Great-circle_distance">geodesic
071 * distance</a> in Meters between this Object and a second Object passed to
072 * this method using <a
073 * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
074 * inverse formula See T Vincenty, "<a
075 * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
076 * Solutions of Geodesics on the Ellipsoid with application of nested
077 * equations</a>", Survey Review, vol XXII no 176, 1975.
078 *
079 * @param location
080 * the destination location
081 */
082 public static double getGeodesicDistance(GeoLocation location, GeoLocation destination) {
083 return vincentyFormula(location, destination, DISTANCE);
084 }
085
086 /**
087 * Calculate <a
088 * href="http://en.wikipedia.org/wiki/Great-circle_distance">geodesic
089 * distance</a> in Meters between this Object and a second Object passed to
090 * this method using <a
091 * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a>
092 * inverse formula See T Vincenty, "<a
093 * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse
094 * Solutions of Geodesics on the Ellipsoid with application of nested
095 * equations</a>", Survey Review, vol XXII no 176, 1975.
096 *
097 * @param location
098 * the destination location
099 * @param formula
100 * This formula calculates initial bearing ({@link #INITIAL_BEARING}),
101 * final bearing ({@link #FINAL_BEARING}) and distance ({@link #DISTANCE}).
102 */
103 private static double vincentyFormula(GeoLocation location, GeoLocation destination, int formula) {
104 double a = 6378137;
105 double b = 6356752.3142;
106 double f = 1 / 298.257223563; // WGS-84 ellipsiod
107 double L = Math.toRadians(destination.getLongitude() - location.getLongitude());
108 double U1 = Math
109 .atan((1 - f) * Math.tan(Math.toRadians(location.getLatitude())));
110 double U2 = Math.atan((1 - f)
111 * Math.tan(Math.toRadians(destination.getLatitude())));
112 double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
113 double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
114
115 double lambda = L;
116 double lambdaP = 2 * Math.PI;
117 double iterLimit = 20;
118 double sinLambda = 0;
119 double cosLambda = 0;
120 double sinSigma = 0;
121 double cosSigma = 0;
122 double sigma = 0;
123 double sinAlpha = 0;
124 double cosSqAlpha = 0;
125 double cos2SigmaM = 0;
126 double C;
127 while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) {
128 sinLambda = Math.sin(lambda);
129 cosLambda = Math.cos(lambda);
130 sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda)
131 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)
132 * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
133 if (sinSigma == 0)
134 return 0; // co-incident points
135 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
136 sigma = Math.atan2(sinSigma, cosSigma);
137 sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
138 cosSqAlpha = 1 - sinAlpha * sinAlpha;
139 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
140 if (Double.isNaN(cos2SigmaM))
141 cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6)
142 C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
143 lambdaP = lambda;
144 lambda = L
145 + (1 - C)
146 * f
147 * sinAlpha
148 * (sigma + C
149 * sinSigma
150 * (cos2SigmaM + C * cosSigma
151 * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
152 }
153 if (iterLimit == 0)
154 return Double.NaN; // formula failed to converge
155
156 double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
157 double A = 1 + uSq / 16384
158 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
159 double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
160 double deltaSigma = B
161 * sinSigma
162 * (cos2SigmaM + B
163 / 4
164 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B
165 / 6 * cos2SigmaM
166 * (-3 + 4 * sinSigma * sinSigma)
167 * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
168 double distance = b * A * (sigma - deltaSigma);
169
170 // initial bearing
171 double fwdAz = Math.toDegrees(Math.atan2(cosU2 * sinLambda, cosU1
172 * sinU2 - sinU1 * cosU2 * cosLambda));
173 // final bearing
174 double revAz = Math.toDegrees(Math.atan2(cosU1 * sinLambda, -sinU1
175 * cosU2 + cosU1 * sinU2 * cosLambda));
176 if (formula == DISTANCE) {
177 return distance;
178 } else if (formula == INITIAL_BEARING) {
179 return fwdAz;
180 } else if (formula == FINAL_BEARING) {
181 return revAz;
182 } else { // should never happpen
183 return Double.NaN;
184 }
185 }
186
187 /**
188 * Returns the <a href="http://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a>
189 * bearing from the current location to the GeoLocation passed in.
190 *
191 * @param location
192 * destination location
193 * @return the bearing in degrees
194 */
195 public static double getRhumbLineBearing(GeoLocation location, GeoLocation destination) {
196 double dLon = Math.toRadians(destination.getLongitude() - location.getLongitude());
197 double dPhi = Math.log(Math.tan(Math.toRadians(destination.getLatitude())
198 / 2 + Math.PI / 4)
199 / Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4));
200 if (Math.abs(dLon) > Math.PI)
201 dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon);
202 return Math.toDegrees(Math.atan2(dLon, dPhi));
203 }
204
205 /**
206 * Returns the <a href="http://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a>
207 * distance from the current location to the GeoLocation passed in.
208 * Ported from <a href="http://www.movable-type.co.uk/">Chris Veness'</a> Javascript Implementation
209 *
210 * @param location
211 * the destination location
212 * @return the distance in Meters
213 */
214 public static double getRhumbLineDistance(GeoLocation location, GeoLocation destination) {
215 double R = 6371; // earth's mean radius in km
216 double dLat = Math.toRadians(destination.getLatitude() - location.getLatitude());
217 double dLon = Math.toRadians(Math.abs(destination.getLongitude()
218 - location.getLongitude()));
219 double dPhi = Math.log(Math.tan(Math.toRadians(destination.getLongitude())
220 / 2 + Math.PI / 4)
221 / Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4));
222 double q = (Math.abs(dLat) > 1e-10) ? dLat / dPhi : Math.cos(Math
223 .toRadians(location.getLatitude()));
224 // if dLon over 180° take shorter rhumb across 180° meridian:
225 if (dLon > Math.PI)
226 dLon = 2 * Math.PI - dLon;
227 double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon);
228 return d * R;
229 }
230
231 }