001/* 002 * Zmanim Java API 003 * Copyright (C) 2004-2025 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: https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html 015 */ 016package com.kosherjava.zmanim.util; 017 018import java.util.Objects; 019import java.util.TimeZone; 020 021/** 022 * A class that contains location information such as latitude and longitude required for astronomical calculations. The 023 * elevation field may not be used by some calculation engines and would be ignored if set. Check the documentation for 024 * specific implementations of the {@link AstronomicalCalculator} to see if elevation is calculated as part of the 025 * algorithm. 026 * 027 * @author © Eliyahu Hershfeld 2004 - 2025 028 */ 029public class GeoLocation implements Cloneable { 030 /** 031 * The latitude, for example 40.096 for Lakewood, NJ. 032 * @see #getLatitude() 033 * @see #setLatitude(double) 034 * @see #setLatitude(int, int, double, String) 035 */ 036 private double latitude; 037 038 /** 039 * The longitude, for example -74.222 for Lakewood, NJ. 040 * @see #getLongitude() 041 * @see #setLongitude(double) 042 * @see #setLongitude(int, int, double, String) 043 */ 044 private double longitude; 045 046 /** 047 * The location name used for display, for example "Lakewood, NJ". 048 * @see #getLocationName() 049 * @see #setLocationName(String) 050 */ 051 private String locationName; 052 053 /** 054 * The location's time zone. 055 * @see #getTimeZone() 056 * @see #setTimeZone(TimeZone) 057 */ 058 private TimeZone timeZone; 059 060 /** 061 * The elevation in Meters <b>above</b> sea level. 062 * @see #getElevation() 063 * @see #setElevation(double) 064 */ 065 private double elevation; 066 067 /** 068 * Constant for a distance type calculation. 069 * @see #getGeodesicDistance(GeoLocation) 070 */ 071 private static final int DISTANCE = 0; 072 073 /** 074 * Constant for an initial bearing type calculation. 075 * @see #getGeodesicInitialBearing(GeoLocation) 076 */ 077 private static final int INITIAL_BEARING = 1; 078 079 /** 080 * Constant for a final bearing type calculation. 081 * @see #getGeodesicFinalBearing(GeoLocation) 082 */ 083 private static final int FINAL_BEARING = 2; 084 085 /** constant for milliseconds in a minute (60,000) */ 086 private static final long MINUTE_MILLIS = 60 * 1000; 087 088 /** constant for milliseconds in an hour (3,600,000) */ 089 private static final long HOUR_MILLIS = MINUTE_MILLIS * 60; 090 091 /** 092 * Method to return the elevation in Meters <b>above</b> sea level. 093 * 094 * @return Returns the elevation in Meters. 095 * @see #setElevation(double) 096 */ 097 public double getElevation() { 098 return elevation; 099 } 100 101 /** 102 * Method to set the elevation in Meters <b>above</b> sea level. 103 * 104 * @param elevation 105 * The elevation to set in Meters. An IllegalArgumentException will be thrown if the value is a negative, 106 * {@link java.lang.Double#isNaN(double) NaN} or {@link java.lang.Double#isInfinite(double) infinite}. 107 */ 108 public void setElevation(double elevation) { 109 if (elevation < 0) { 110 throw new IllegalArgumentException("Elevation cannot be negative"); 111 } 112 if (Double.isNaN(elevation) || Double.isInfinite(elevation)) { 113 throw new IllegalArgumentException("Elevation cannot be NaN or infinite"); 114 } 115 this.elevation = elevation; 116 } 117 118 /** 119 * GeoLocation constructor with parameters for all required fields. 120 * 121 * @param name 122 * The location name for display, for example "Lakewood, NJ". 123 * @param latitude 124 * the latitude as a <code>double</code>, for example 40.096 for Lakewood, NJ. 125 * <b>Note:</b> For latitudes south of the equator, a negative value should be used. 126 * @param longitude 127 * the longitude as a <code>double</code>, for example -74.222 for Lakewood, NJ. <b>Note:</b> For longitudes 128 * east of the <a href="https://en.wikipedia.org/wiki/Prime_Meridian">Prime Meridian</a> (Greenwich), 129 * a negative value should be used. 130 * @param timeZone 131 * the <code>TimeZone</code> for the location. 132 */ 133 public GeoLocation(String name, double latitude, double longitude, TimeZone timeZone) { 134 this(name, latitude, longitude, 0, timeZone); 135 } 136 137 /** 138 * GeoLocation constructor with parameters for all required fields. 139 * 140 * @param name 141 * The location name for display, for example "Lakewood, NJ". 142 * @param latitude 143 * the latitude as a <code>double</code>, for example 40.096 for Lakewood, NJ. 144 * <b>Note:</b> For latitudes south of the equator, a negative value should be used. 145 * @param longitude 146 * double the longitude as a <code>double</code>, for example -74.222 for Lakewood, NJ. 147 * <b>Note:</b> For longitudes east of the <a href="https://en.wikipedia.org/wiki/Prime_Meridian">Prime 148 * Meridian</a> (Greenwich), a negative value should be used. 149 * @param elevation 150 * the elevation above sea level in Meters. 151 * @param timeZone 152 * the <code>TimeZone</code> for the location. 153 */ 154 public GeoLocation(String name, double latitude, double longitude, double elevation, TimeZone timeZone) { 155 setLocationName(name); 156 setLatitude(latitude); 157 setLongitude(longitude); 158 setElevation(elevation); 159 setTimeZone(timeZone); 160 } 161 162 /** 163 * Default GeoLocation constructor will set location to the Prime Meridian at Greenwich, England and a TimeZone of 164 * GMT. The longitude will be set to 0 and the latitude will be 51.4772 to match the location of the <a 165 * href="https://www.rmg.co.uk/royal-observatory">Royal Observatory, Greenwich</a>. No daylight savings time will be used. 166 */ 167 public GeoLocation() { 168 setLocationName("Greenwich, England"); 169 setLongitude(0); // added for clarity 170 setLatitude(51.4772); 171 setTimeZone(TimeZone.getTimeZone("GMT")); 172 } 173 174 /** 175 * Method to set the latitude as a <code>double</code>, for example 40.096 for Lakewood, NJ. 176 * 177 * @param latitude 178 * The degrees of latitude to set. The values should be between -90° and 90°. For example 40.096 179 * would be used for Lakewood, NJ. <b>Note:</b> For latitudes south of the equator, a negative value 180 * should be used. An IllegalArgumentException will be thrown if the value exceeds the limits. 181 */ 182 public void setLatitude(double latitude) { 183 if (latitude > 90 || latitude < -90 || Double.isNaN(latitude)) { 184 throw new IllegalArgumentException("Latitude must be between -90 and 90"); 185 } 186 this.latitude = latitude; 187 } 188 189 /** 190 * Method to set the latitude in degrees, minutes and seconds. For example, set the degrees to 40, minutes to 5, 191 * seconds to 45.6 and direction to "N" for Lakewood, NJ. 192 * 193 * @param degrees 194 * The degrees of latitude to set between 0° and 90°, for example 40 would be used for Lakewood, NJ. 195 * An IllegalArgumentException will be thrown if the value exceeds the limit. 196 * @param minutes 197 * <a href="https://en.wikipedia.org/wiki/Minute_of_arc#Cartography">minutes of arc</a>, for example 5 198 * would be used for Lakewood, NJ. 199 * @param seconds 200 * <a href="https://en.wikipedia.org/wiki/Minute_of_arc#Cartography">seconds of arc</a>, for example 45.6 201 * would be used for Lakewood, NJ. 202 * @param direction 203 * "N" for north and "S" for south, for example "N" would be used for Lakewood, NJ. An 204 * IllegalArgumentException will be thrown if the value is not "S" or "N". 205 */ 206 public void setLatitude(int degrees, int minutes, double seconds, String direction) { 207 double tempLat = degrees + ((minutes + (seconds / 60.0)) / 60.0); 208 if (tempLat > 90 || tempLat < 0 || Double.isNaN(tempLat)) { //FIXME An exception should be thrown if degrees, minutes or seconds are negative 209 throw new IllegalArgumentException( 210 "Latitude must be between 0 and 90. Use direction of S instead of negative."); 211 } 212 if (direction.equals("S")) { 213 tempLat *= -1; 214 } else if (!direction.equals("N")) { 215 throw new IllegalArgumentException("Latitude direction must be N or S"); 216 } 217 this.latitude = tempLat; 218 } 219 220 /** 221 * Method to return the latitude as a <code>double</code>, for example 40.096 for Lakewood, NJ. 222 * @return Returns the latitude. 223 */ 224 public double getLatitude() { 225 return latitude; 226 } 227 228 /** 229 * Method to set the longitude as a <code>double</code>, for example -74.222 for Lakewood, NJ 230 * 231 * @param longitude 232 * The degrees of longitude to set as a <code>double</code> between -180.0° and 180.0°. For example 233 * -74.222 would be used for Lakewood, NJ. Note: for longitudes east of the <a href= 234 * "https://en.wikipedia.org/wiki/Prime_Meridian">Prime Meridian</a> (Greenwich) a negative value 235 * should be used.An IllegalArgumentException will be thrown if the value exceeds the limit. 236 */ 237 public void setLongitude(double longitude) { 238 if (longitude > 180 || longitude < -180 || Double.isNaN(longitude)) { 239 throw new IllegalArgumentException("Longitude must be between -180 and 180"); 240 } 241 this.longitude = longitude; 242 } 243 244 /** 245 * Method to set the longitude in degrees, minutes and seconds. For example, set the degrees to 74, minutes to 13, 246 * seconds to 19.2 and direction to "W" for Lakewood, NJ. 247 * 248 * @param degrees 249 * The degrees of longitude to set between 0° and 180°. For example 74 would be set for Lakewood, NJ. 250 * An IllegalArgumentException will be thrown if the value exceeds the limits. 251 * @param minutes 252 * <a href="https://en.wikipedia.org/wiki/Minute_of_arc#Cartography">minutes of arc</a>. For example 13 253 * would be set for Lakewood, NJ. 254 * @param seconds 255 * <a href="https://en.wikipedia.org/wiki/Minute_of_arc#Cartography">seconds of arc</a>. For example 19.2 256 * would be set for Lakewood, NJ. 257 * @param direction 258 * "E" for east of the <a href="https://en.wikipedia.org/wiki/Prime_Meridian">Prime Meridian</a> 259 * or "W"for west of it. For example, "W" would be set for Lakewood, NJ. 260 * An IllegalArgumentException will be thrown if the value is not E or W. 261 */ 262 public void setLongitude(int degrees, int minutes, double seconds, String direction) { 263 double longTemp = degrees + ((minutes + (seconds / 60.0)) / 60.0); 264 if (longTemp > 180 || this.longitude < 0 || Double.isNaN(longTemp)) { //FIXME An exception should be thrown if degrees, minutes or seconds are negative 265 throw new IllegalArgumentException("Longitude must be between 0 and 180. Use a direction of W instead of negative."); 266 } 267 if (direction.equals("W")) { 268 longTemp *= -1; 269 } else if (!direction.equals("E")) { 270 throw new IllegalArgumentException("Longitude direction must be E or W"); 271 } 272 this.longitude = longTemp; 273 } 274 275 /** 276 * Method to return the longitude as a <code>double</code>, for example -74.222 for Lakewood, NJ. 277 * @return Returns the longitude. 278 */ 279 public double getLongitude() { 280 return longitude; 281 } 282 283 /** 284 * Method to return the location name (for display), for example "Lakewood, NJ". 285 * @return Returns the location name. 286 */ 287 public String getLocationName() { 288 return locationName; 289 } 290 291 /** 292 * Setter for the location name (for display), for example "Lakewood, NJ". 293 * @param name 294 * The setter method for the display name. 295 */ 296 public void setLocationName(String name) { 297 this.locationName = name; 298 } 299 300 /** 301 * Method to return the time zone. 302 * @return Returns the timeZone. 303 */ 304 public TimeZone getTimeZone() { 305 return timeZone; 306 } 307 308 /** 309 * Method to set the TimeZone. If this is ever set after the GeoLocation is set in the 310 * {@link com.kosherjava.zmanim.AstronomicalCalendar}, it is critical that 311 * {@link com.kosherjava.zmanim.AstronomicalCalendar#getCalendar()}. 312 * {@link java.util.Calendar#setTimeZone(TimeZone) setTimeZone(TimeZone)} be called in order for the 313 * AstronomicalCalendar to output times in the expected offset. This situation will arise if the 314 * AstronomicalCalendar is ever {@link com.kosherjava.zmanim.AstronomicalCalendar#clone() cloned}. 315 * 316 * @param timeZone 317 * The timeZone to set. 318 */ 319 public void setTimeZone(TimeZone timeZone) { 320 this.timeZone = timeZone; 321 } 322 323 /** 324 * A method that will return the location's local mean time offset in milliseconds from local <a 325 * href="https://en.wikipedia.org/wiki/Standard_time">standard time</a>. The globe is split into 360°, with 326 * 15° per hour of the day. For a local that is at a longitude that is evenly divisible by 15 (longitude % 15 == 327 * 0), at solar {@link com.kosherjava.zmanim.AstronomicalCalendar#getSunTransit() noon} (with adjustment for the <a 328 * href="https://en.wikipedia.org/wiki/Equation_of_time">equation of time</a>) the sun should be directly overhead, 329 * so a user who is 1° west of this will have noon at 4 minutes after standard time noon, and conversely, a user 330 * who is 1° east of the 15° longitude will have noon at 11:56 AM. Lakewood, N.J., whose longitude is 331 * -74.222, is 0.778 away from the closest multiple of 15 at -75°. This is multiplied by 4 to yield 3 minutes 332 * and 10 seconds earlier than standard time. The offset returned does not account for the <a 333 * href="https://en.wikipedia.org/wiki/Daylight_saving_time">Daylight saving time</a> offset since this class is 334 * unaware of dates. 335 * 336 * @return the offset in milliseconds not accounting for Daylight saving time. A positive value will be returned 337 * East of the 15° timezone line, and a negative value West of it. 338 */ 339 public long getLocalMeanTimeOffset() { 340 return (long) (getLongitude() * 4 * MINUTE_MILLIS - getTimeZone().getRawOffset()); 341 } 342 343 /** 344 * Adjust the date for <a href="https://en.wikipedia.org/wiki/180th_meridian">antimeridian</a> crossover. This is 345 * needed to deal with edge cases such as Samoa that use a different calendar date than expected based on their 346 * geographic location. 347 * 348 * The actual Time Zone offset may deviate from the expected offset based on the longitude. Since the 'absolute time' 349 * calculations are always based on longitudinal offset from UTC for a given date, the date is presumed to only 350 * increase East of the Prime Meridian, and to only decrease West of it. For Time Zones that cross the antimeridian, 351 * the date will be artificially adjusted before calculation to conform with this presumption. 352 * 353 * For example, Apia, Samoa with a longitude of -171.75 uses a local offset of +14:00. When calculating sunrise for 354 * 2018-02-03, the calculator should operate using 2018-02-02 since the expected zone is -11. After determining the 355 * UTC time, the local DST offset of <a href="https://en.wikipedia.org/wiki/UTC%2B14:00">UTC+14:00</a> should be applied 356 * to bring the date back to 2018-02-03. 357 * 358 * @return the number of days to adjust the date This will typically be 0 unless the date crosses the antimeridian 359 */ 360 public int getAntimeridianAdjustment() { 361 double localHoursOffset = getLocalMeanTimeOffset() / (double)HOUR_MILLIS; 362 363 if (localHoursOffset >= 20){// if the offset is 20 hours or more in the future (never expected anywhere other 364 // than a location using a timezone across the antimeridian to the east such as Samoa) 365 return 1; // roll the date forward a day 366 } else if (localHoursOffset <= -20) { // if the offset is 20 hours or more in the past (no current location is known 367 //that crosses the antimeridian to the west, but better safe than sorry) 368 return -1; // roll the date back a day 369 } 370 return 0; //99.999% of the world will have no adjustment 371 } 372 373 /** 374 * Calculate the initial <a href="https://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing between this 375 * Object and a second Object passed to this method using <a 376 * href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, "<a 377 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 378 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975 379 * 380 * @param location 381 * the destination location 382 * @return the initial bearing 383 */ 384 public double getGeodesicInitialBearing(GeoLocation location) { 385 return vincentyInverseFormula(location, INITIAL_BEARING); 386 } 387 388 /** 389 * Calculate the final <a href="https://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing between this Object 390 * and a second Object passed to this method using <a href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus 391 * Vincenty's</a> inverse formula See T Vincenty, "<a href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and 392 * Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations</a>", Survey Review, vol 393 * XXII no 176, 1975 394 * 395 * @param location 396 * the destination location 397 * @return the final bearing 398 */ 399 public double getGeodesicFinalBearing(GeoLocation location) { 400 return vincentyInverseFormula(location, FINAL_BEARING); 401 } 402 403 /** 404 * Calculate <a href="https://en.wikipedia.org/wiki/Great-circle_distance">geodesic distance</a> in Meters between 405 * this Object and a second Object passed to this method using <a 406 * href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, "<a 407 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 408 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975 409 * 410 * @see #vincentyInverseFormula(GeoLocation, int) 411 * @param location 412 * the destination location 413 * @return the geodesic distance in Meters 414 */ 415 public double getGeodesicDistance(GeoLocation location) { 416 return vincentyInverseFormula(location, DISTANCE); 417 } 418 419 /** 420 * Calculate the destination point based on an initial bearing and distance in meters from the current location using 421 * <a href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> direct formula. See T Vincenty, "<a 422 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 423 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975. 424 * 425 * @param initialBearing the initialBearing 426 * @param distance the distance in meters. 427 * @return the GeoLocation containing the destination point. The Timezone is set to the origin point (current object). 428 */ 429 private GeoLocation vincentyDirectFormulaDestination(double initialBearing, double distance) { 430 double major_semi_axis = 6378137; 431 double minor_semi_axis = 6356752.3142; 432 double flattening = 1 / 298.257223563; // WGS-84 ellipsoid 433 double initial_bearing_radians = Math.toRadians(initialBearing); 434 double sinAzimuth1 = Math.sin(initial_bearing_radians); 435 double cosAzimuth1 = Math.cos(initial_bearing_radians); 436 double tanU1 = (1 - flattening) * Math.tan(Math.toRadians(getLatitude())); 437 double cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)); 438 double sinU1 = tanU1 * cosU1; 439 double eq_p1_ang_dist = Math.atan2(tanU1, cosAzimuth1); // eq_p1_ang_dist = angular distance on the sphere from the equator to P1 440 double sinAzimuth = cosU1 * sinAzimuth1; //azimuth of the geodesic at the equator 441 double cosSqAzimuth = 1 - sinAzimuth*sinAzimuth; 442 double uSq = cosSqAzimuth * (Math.pow(major_semi_axis, 2) - Math.pow(minor_semi_axis, 2) / Math.pow(minor_semi_axis, 2)); 443 double a = 1 + uSq/16384*(4096 + uSq *(-768 + uSq *(320 - 175 * uSq))); 444 double b = uSq / 1024 * (256 + uSq *(-128 + uSq * (74-47 * uSq))); 445 double p1_p2_ang_dist = distance / (minor_semi_axis * a); //p1_p2_ang_dist = angular distance P1 P2 on the sphere 446 double sinSigma = Double.NaN; 447 double cosSigma = Double.NaN; 448 double cos2_eq_mid_ang_distance = Double.NaN; // # eq_mid_ang_distance = angular distance on the sphere from the equator to the midpoint of the line 449 double a_prime = Double.NaN; 450 int iterations = 0; 451 452 do { 453 cos2_eq_mid_ang_distance = Math.cos(2*eq_p1_ang_dist + p1_p2_ang_dist); 454 sinSigma = Math.sin(p1_p2_ang_dist); 455 cosSigma = Math.cos(p1_p2_ang_dist); 456 double delta_ang_distance = b * sinSigma * (cos2_eq_mid_ang_distance + b / 4 * 457 (cosSigma * (-1 + 2 * cos2_eq_mid_ang_distance * cos2_eq_mid_ang_distance) - b / 6 * cos2_eq_mid_ang_distance * 458 (-3+4*sinSigma*sinSigma)*(-3+4*cos2_eq_mid_ang_distance*cos2_eq_mid_ang_distance))); 459 a_prime = p1_p2_ang_dist; 460 p1_p2_ang_dist = distance / (minor_semi_axis * a) + delta_ang_distance; 461 } while (Math.abs(p1_p2_ang_dist-a_prime) > 1e-12 && ++iterations < 100); // iterate until negligible change in lambda (about 0.006mm) 462 463 double x = sinU1 * sinSigma - cosU1 * cosSigma * cosAzimuth1; 464 double other_latitude = Math.toDegrees(Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAzimuth1, (1 - flattening) * Math.sqrt(sinAzimuth*sinAzimuth + x * x))); 465 double lambda = Math.atan2(sinSigma*sinAzimuth1, cosU1*cosSigma - sinU1*sinSigma*cosAzimuth1); 466 double c = flattening/16*cosSqAzimuth*(4+flattening*(4-3*cosSqAzimuth)); 467 double l = lambda - (1-c) * flattening * sinAzimuth * 468 (p1_p2_ang_dist + c * sinSigma * (cos2_eq_mid_ang_distance + c * cosSigma * (-1 + 2 * cos2_eq_mid_ang_distance * cos2_eq_mid_ang_distance))); 469 double other_longitude = longitude + Math.toDegrees(l); 470 471 return new GeoLocation("Destination", other_latitude, other_longitude, getTimeZone()); //ToDo - we can easily return final_bearing, it just needs some minor refactoring 472 } 473 474 /** 475 * Calculate <a href="https://en.wikipedia.org/wiki/Great-circle_distance">geodesic distance</a> in Meters between 476 * this Object and a second Object passed to this method using <a 477 * href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, "<a 478 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 479 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975 480 * 481 * @param location 482 * the destination location 483 * @param formula 484 * This formula calculates initial bearing ({@link #INITIAL_BEARING}), final bearing ( 485 * {@link #FINAL_BEARING}) and distance ({@link #DISTANCE}). 486 * @return geodesic distance in Meters 487 */ 488 private double vincentyInverseFormula(GeoLocation location, int formula) { 489 double majorSemiAxis = 6378137; 490 double minorSemiAxis = 6356752.3142; 491 double f = 1 / 298.257223563; // WGS-84 ellipsoid 492 double L = Math.toRadians(location.getLongitude() - getLongitude()); 493 double U1 = Math.atan((1 - f) * Math.tan(Math.toRadians(getLatitude()))); 494 double U2 = Math.atan((1 - f) * Math.tan(Math.toRadians(location.getLatitude()))); 495 double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1); 496 double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2); 497 498 double lambda = L; 499 double lambdaP = 2 * Math.PI; 500 double iterLimit = 20; 501 double sinLambda = 0; 502 double cosLambda = 0; 503 double sinSigma = 0; 504 double cosSigma = 0; 505 double sigma = 0; 506 double sinAlpha = 0; 507 double cosSqAlpha = 0; 508 double cos2SigmaM = 0; 509 double C; 510 while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) { 511 sinLambda = Math.sin(lambda); 512 cosLambda = Math.cos(lambda); 513 sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) 514 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)); 515 if (sinSigma == 0) 516 return 0; // co-incident points 517 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda; 518 sigma = Math.atan2(sinSigma, cosSigma); 519 sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; 520 cosSqAlpha = 1 - sinAlpha * sinAlpha; 521 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha; 522 if (Double.isNaN(cos2SigmaM)) 523 cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6) 524 C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)); 525 lambdaP = lambda; 526 lambda = L + (1 - C) * f * sinAlpha 527 * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM))); 528 } 529 if (iterLimit == 0) 530 return Double.NaN; // formula failed to converge 531 532 double uSq = cosSqAlpha * (majorSemiAxis * majorSemiAxis - minorSemiAxis * minorSemiAxis) / (minorSemiAxis * minorSemiAxis); 533 double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))); 534 double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); 535 double deltaSigma = B 536 * sinSigma 537 * (cos2SigmaM + B 538 / 4 539 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM 540 * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM))); 541 double distance = minorSemiAxis * A * (sigma - deltaSigma); 542 543 // initial bearing 544 double fwdAz = Math.toDegrees(Math.atan2(cosU2 * sinLambda, cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)); 545 // final bearing 546 double revAz = Math.toDegrees(Math.atan2(cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda)); 547 if (formula == DISTANCE) { 548 return distance; 549 } else if (formula == INITIAL_BEARING) { 550 return fwdAz; 551 } else if (formula == FINAL_BEARING) { 552 return revAz; 553 } else { // should never happen 554 return Double.NaN; 555 } 556 } 557 558 /** 559 * Returns the <a href="https://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a> bearing from the current location to 560 * the GeoLocation passed in. 561 * 562 * @param location 563 * destination location 564 * @return the bearing in degrees 565 */ 566 public double getRhumbLineBearing(GeoLocation location) { 567 double dLon = Math.toRadians(location.getLongitude() - getLongitude()); 568 double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4) 569 / Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4)); 570 if (Math.abs(dLon) > Math.PI) 571 dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon); 572 return Math.toDegrees(Math.atan2(dLon, dPhi)); 573 } 574 575 /** 576 * Returns the <a href="https://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a> distance from the current location 577 * to the GeoLocation passed in. 578 * 579 * @param location 580 * the destination location 581 * @return the distance in Meters 582 */ 583 public double getRhumbLineDistance(GeoLocation location) { 584 double earthRadius = 6378137; // Earth's radius in meters (WGS-84) 585 double dLat = Math.toRadians(location.getLatitude()) - Math.toRadians(getLatitude()); 586 double dLon = Math.abs(Math.toRadians(location.getLongitude()) - Math.toRadians(getLongitude())); 587 double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4) 588 / Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4)); 589 double q = dLat / dPhi; 590 591 if (!(Math.abs(q) <= Double.MAX_VALUE)) { 592 q = Math.cos(Math.toRadians(getLatitude())); 593 } 594 // if dLon over 180° take shorter rhumb across 180° meridian: 595 if (dLon > Math.PI) { 596 dLon = 2 * Math.PI - dLon; 597 } 598 double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon); 599 return d * earthRadius; 600 } 601 602 /** 603 * A method that returns an XML formatted <code>String</code> representing the serialized <code>Object</code>. Very 604 * similar to the toString method but the return value is in an xml format. The format currently used (subject to 605 * change) is: 606 * 607 * <pre> 608 * <GeoLocation> 609 * <LocationName>Lakewood, NJ</LocationName> 610 * <Latitude>40.096&deg</Latitude> 611 * <Longitude>-74.222amp;deg</Longitude> 612 * <Elevation>0 Meters</Elevation> 613 * <TimezoneName>America/New_York</TimezoneName> 614 * <TimeZoneDisplayName>Eastern Standard Time</TimeZoneDisplayName> 615 * <TimezoneGMTOffset>-5</TimezoneGMTOffset> 616 * <TimezoneDSTOffset>1</TimezoneDSTOffset> 617 * </GeoLocation> 618 * </pre> 619 * 620 * @return The XML formatted <code>String</code>. 621 */ 622 public String toXML() { 623 return "<GeoLocation>\n" + 624 "\t<LocationName>" + getLocationName() + "</LocationName>\n" + 625 "\t<Latitude>" + getLatitude() + "</Latitude>\n" + 626 "\t<Longitude>" + getLongitude() + "</Longitude>\n" + 627 "\t<Elevation>" + getElevation() + " Meters" + "</Elevation>\n" + 628 "\t<TimezoneName>" + getTimeZone().getID() + "</TimezoneName>\n" + 629 "\t<TimeZoneDisplayName>" + getTimeZone().getDisplayName() + "</TimeZoneDisplayName>\n" + 630 "\t<TimezoneGMTOffset>" + getTimeZone().getRawOffset() / HOUR_MILLIS + 631 "</TimezoneGMTOffset>\n" + 632 "\t<TimezoneDSTOffset>" + getTimeZone().getDSTSavings() / HOUR_MILLIS + 633 "</TimezoneDSTOffset>\n" + 634 "</GeoLocation>"; 635 } 636 637 /** 638 * @see java.lang.Object#equals(Object) 639 */ 640 public boolean equals(Object object) { 641 if (this == object) 642 return true; 643 if (!(object instanceof GeoLocation)) 644 return false; 645 GeoLocation geo = (GeoLocation) object; 646 return Double.doubleToLongBits(this.latitude) == Double.doubleToLongBits(geo.latitude) 647 && Double.doubleToLongBits(this.longitude) == Double.doubleToLongBits(geo.longitude) 648 && this.elevation == geo.elevation 649 && (Objects.equals(this.locationName, geo.locationName)) 650 && (Objects.equals(this.timeZone, geo.timeZone)); 651 } 652 653 /** 654 * @see java.lang.Object#hashCode() 655 */ 656 public int hashCode() { 657 658 int result = 17; 659 long latLong = Double.doubleToLongBits(this.latitude); 660 long lonLong = Double.doubleToLongBits(this.longitude); 661 long elevLong = Double.doubleToLongBits(this.elevation); 662 int latInt = (int) (latLong ^ (latLong >>> 32)); 663 int lonInt = (int) (lonLong ^ (lonLong >>> 32)); 664 int elevInt = (int) (elevLong ^ (elevLong >>> 32)); 665 result = 37 * result + getClass().hashCode(); 666 result += 37 * result + latInt; 667 result += 37 * result + lonInt; 668 result += 37 * result + elevInt; 669 result += 37 * result + (this.locationName == null ? 0 : this.locationName.hashCode()); 670 result += 37 * result + (this.timeZone == null ? 0 : this.timeZone.hashCode()); 671 return result; 672 } 673 674 /** 675 * @see java.lang.Object#toString() 676 */ 677 public String toString() { 678 return "\nLocation Name:\t\t\t" + getLocationName() + 679 "\nLatitude:\t\t\t" + getLatitude() + "\u00B0" + 680 "\nLongitude:\t\t\t" + getLongitude() + "\u00B0" + 681 "\nElevation:\t\t\t" + getElevation() + " Meters" + 682 "\nTimezone ID:\t\t\t" + getTimeZone().getID() + 683 "\nTimezone Display Name:\t\t" + getTimeZone().getDisplayName() + 684 " (" + getTimeZone().getDisplayName(false, TimeZone.SHORT) + ")" + 685 "\nTimezone GMT Offset:\t\t" + getTimeZone().getRawOffset() / HOUR_MILLIS + 686 "\nTimezone DST Offset:\t\t" + getTimeZone().getDSTSavings() / HOUR_MILLIS; 687 } 688 689 /** 690 * An implementation of the {@link java.lang.Object#clone()} method that creates a <a 691 * href="https://en.wikipedia.org/wiki/Object_copy#Deep_copy">deep copy</a> of the object. 692 * <b>Note:</b> If the {@link java.util.TimeZone} in the clone will be changed from the original, it is critical 693 * that {@link com.kosherjava.zmanim.AstronomicalCalendar#getCalendar()}. 694 * {@link java.util.Calendar#setTimeZone(TimeZone) setTimeZone(TimeZone)} is called after cloning in order for the 695 * AstronomicalCalendar to output times in the expected offset. 696 * 697 * @see java.lang.Object#clone() 698 */ 699 public Object clone() { 700 GeoLocation clone = null; 701 try { 702 clone = (GeoLocation) super.clone(); 703 } catch (CloneNotSupportedException cnse) { 704 //Required by the compiler. Should never be reached since we implement clone() 705 } 706 if (clone != null) { 707 clone.timeZone = (TimeZone) getTimeZone().clone(); 708 clone.locationName = getLocationName(); 709 } 710 return clone; 711 } 712}