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 }