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