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 }