001/* 002 * Zmanim Java API 003 * Copyright (C) 2004-2022 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 com.kosherjava.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 - 2022 024 * @deprecated All methods in this call are available in the {@link GeoLocation} class, and this class that duplicates that 025 * code will be removed in release 3.0. 026 */ 027@Deprecated //(since="2.0.3", forRemoval=true) // add back once Java 9 is the minimum supported version 028public class GeoLocationUtils { 029 /** 030 * Constant for a distance type calculation. 031 * @see #getGeodesicDistance(GeoLocation, GeoLocation) 032 */ 033 private static int DISTANCE = 0; 034 035 /** 036 * Constant for a initial bearing type calculation. 037 * @see #getGeodesicInitialBearing(GeoLocation, GeoLocation) 038 */ 039 private static int INITIAL_BEARING = 1; 040 041 /** 042 * Constant for a final bearing type calculation. 043 * @see #getGeodesicFinalBearing(GeoLocation, GeoLocation) 044 */ 045 private static int FINAL_BEARING = 2; 046 047 /** 048 * Calculate the <a href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> initial bearing between this Object and 049 * a second Object passed to this method using <a href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus 050 * Vincenty's</a> inverse formula See T Vincenty, "<a href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and 051 * Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations</a>", Survey Review, vol XXII 052 * no 176, 1975. 053 * 054 * @param location 055 * the initial location 056 * @param destination 057 * the destination location 058 * @return the geodesic bearing 059 */ 060 public static double getGeodesicInitialBearing(GeoLocation location, GeoLocation destination) { 061 return vincentyFormula(location, destination, INITIAL_BEARING); 062 } 063 064 /** 065 * Calculate the <a href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> final bearing between this Object 066 * and a second Object passed to this method using <a href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> 067 * inverse formula See T Vincenty, "<a href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics 068 * on the Ellipsoid with application of nested equations</a>", Survey Review, vol XXII no 176, 1975. 069 * 070 * @param location 071 * the initial location 072 * @param destination 073 * the destination location 074 * @return the geodesic bearing 075 */ 076 public static double getGeodesicFinalBearing(GeoLocation location, GeoLocation destination) { 077 return vincentyFormula(location, destination, FINAL_BEARING); 078 } 079 080 /** 081 * Calculate <a href="http://en.wikipedia.org/wiki/Great-circle_distance">geodesic distance</a> in Meters 082 * between this Object and a second Object passed to this method using <a 083 * href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, 084 * "<a href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the 085 * Ellipsoid with application of nested equations</a>", Survey Review, vol XXII no 176, 1975. This uses the 086 * WGS-84 geodetic model. 087 * 088 * @param location 089 * the initial location 090 * @param destination 091 * the destination location 092 * @return the geodesic distance in Meters 093 */ 094 public static double getGeodesicDistance(GeoLocation location, GeoLocation destination) { 095 return vincentyFormula(location, destination, DISTANCE); 096 } 097 098 /** 099 * Calculates the initial <a href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing, final bearing or 100 * <a href="http://en.wikipedia.org/wiki/Great-circle_distance">geodesic distance</a> using <a href= 101 * "http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, "<a 102 * href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 103 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975. 104 * 105 * @param location 106 * the initial location 107 * @param destination 108 * the destination location 109 * @param formula 110 * This formula calculates initial bearing ({@link #INITIAL_BEARING}), 111 * final bearing ({@link #FINAL_BEARING}) and distance ({@link #DISTANCE}). 112 * @return 113 * the geodesic distance, initial or final bearing (based on the formula passed in) between the location 114 * and destination in Meters 115 * @see #getGeodesicDistance(GeoLocation, GeoLocation) 116 * @see #getGeodesicInitialBearing(GeoLocation, GeoLocation) 117 * @see #getGeodesicFinalBearing(GeoLocation, GeoLocation) 118 */ 119 private static double vincentyFormula(GeoLocation location, GeoLocation destination, int formula) { 120 double a = 6378137; // length of semi-major axis of the ellipsoid (radius at equator) in metres based on WGS-84 121 double b = 6356752.3142; // length of semi-minor axis of the ellipsoid (radius at the poles) in meters based on WGS-84 122 double f = 1 / 298.257223563; // flattening of the ellipsoid based on WGS-84 123 double L = Math.toRadians(destination.getLongitude() - location.getLongitude()); //difference in longitude of two points; 124 double U1 = Math.atan((1 - f) * Math.tan(Math.toRadians(location.getLatitude()))); // reduced latitude (latitude on the auxiliary sphere) 125 double U2 = Math.atan((1 - f) * Math.tan(Math.toRadians(destination.getLatitude()))); // reduced latitude (latitude on the auxiliary sphere) 126 127 double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1); 128 double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2); 129 130 double lambda = L; 131 double lambdaP = 2 * Math.PI; 132 double iterLimit = 20; 133 double sinLambda = 0; 134 double cosLambda = 0; 135 double sinSigma = 0; 136 double cosSigma = 0; 137 double sigma = 0; 138 double sinAlpha = 0; 139 double cosSqAlpha = 0; 140 double cos2SigmaM = 0; 141 double C; 142 while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) { 143 sinLambda = Math.sin(lambda); 144 cosLambda = Math.cos(lambda); 145 sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) 146 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) 147 * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)); 148 if (sinSigma == 0) 149 return 0; // co-incident points 150 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda; 151 sigma = Math.atan2(sinSigma, cosSigma); 152 sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; 153 cosSqAlpha = 1 - sinAlpha * sinAlpha; 154 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha; 155 if (Double.isNaN(cos2SigmaM)) 156 cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6) 157 C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)); 158 lambdaP = lambda; 159 lambda = L 160 + (1 - C) 161 * f 162 * sinAlpha 163 * (sigma + C 164 * sinSigma 165 * (cos2SigmaM + C * cosSigma 166 * (-1 + 2 * cos2SigmaM * cos2SigmaM))); 167 } 168 if (iterLimit == 0) 169 return Double.NaN; // formula failed to converge 170 171 double uSq = cosSqAlpha * (a * a - b * b) / (b * b); 172 double A = 1 + uSq / 16384 173 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))); 174 double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); 175 double deltaSigma = B 176 * sinSigma 177 * (cos2SigmaM + B 178 / 4 179 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B 180 / 6 * cos2SigmaM 181 * (-3 + 4 * sinSigma * sinSigma) 182 * (-3 + 4 * cos2SigmaM * cos2SigmaM))); 183 double distance = b * A * (sigma - deltaSigma); 184 185 // initial bearing 186 double fwdAz = Math.toDegrees(Math.atan2(cosU2 * sinLambda, cosU1 187 * sinU2 - sinU1 * cosU2 * cosLambda)); 188 // final bearing 189 double revAz = Math.toDegrees(Math.atan2(cosU1 * sinLambda, -sinU1 190 * cosU2 + cosU1 * sinU2 * cosLambda)); 191 if (formula == DISTANCE) { 192 return distance; 193 } else if (formula == INITIAL_BEARING) { 194 return fwdAz; 195 } else if (formula == FINAL_BEARING) { 196 return revAz; 197 } else { // should never happen 198 return Double.NaN; 199 } 200 } 201 202 /** 203 * Returns the <a href="http://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a> 204 * bearing from the current location to the GeoLocation passed in. 205 * 206 * @param location 207 * the initial location 208 * @param destination 209 * the destination location 210 * @return the bearing in degrees 211 */ 212 public static double getRhumbLineBearing(GeoLocation location, GeoLocation destination) { 213 double dLon = Math.toRadians(destination.getLongitude() - location.getLongitude()); 214 double dPhi = Math.log(Math.tan(Math.toRadians(destination.getLatitude()) 215 / 2 + Math.PI / 4) 216 / Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4)); 217 if (Math.abs(dLon) > Math.PI) 218 dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon); 219 return Math.toDegrees(Math.atan2(dLon, dPhi)); 220 } 221 222 /** 223 * Returns the <a href="http://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a> distance between two GeoLocations 224 * passed in. Ported from <a href="http://www.movable-type.co.uk/">Chris Veness'</a> Javascript Implementation. 225 * 226 * @param location 227 * the initial location 228 * @param destination 229 * the destination location 230 * @return the distance in Meters 231 */ 232 public static double getRhumbLineDistance(GeoLocation location, GeoLocation destination) { 233 double earthRadius = 6378137; // Earth's radius in meters (WGS-84) 234 double dLat = Math.toRadians(location.getLatitude()) - Math.toRadians(destination.getLatitude()); 235 double dLon = Math.abs(Math.toRadians(location.getLongitude()) - Math.toRadians(destination.getLongitude())); 236 double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4) 237 / Math.tan(Math.toRadians(destination.getLatitude()) / 2 + Math.PI / 4)); 238 double q = dLat / dPhi; 239 240 if (!(Math.abs(q) <= Double.MAX_VALUE)) { 241 q = Math.cos(Math.toRadians(destination.getLatitude())); 242 } 243 // if dLon over 180° take shorter rhumb across 180° meridian: 244 if (dLon > Math.PI) { 245 dLon = 2 * Math.PI - dLon; 246 } 247 double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon); 248 return d * earthRadius; 249 } 250}