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