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: https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html 015 */ 016package com.kosherjava.zmanim.util; 017 018import java.util.TimeZone; 019 020/** 021 * A class that contains location information such as latitude and longitude required for astronomical calculations. The 022 * elevation field may not be used by some calculation engines and would be ignored if set. Check the documentation for 023 * specific implementations of the {@link AstronomicalCalculator} to see if elevation is calculated as part of the 024 * algorithm. 025 * 026 * @author © Eliyahu Hershfeld 2004 - 2022 027 */ 028public class GeoLocation implements Cloneable { 029 /** 030 * @see #getLatitude() 031 * @see #setLatitude(double) 032 * @see #setLatitude(int, int, double, String) 033 */ 034 private double latitude; 035 036 /** 037 * @see #getLongitude() 038 * @see #setLongitude(double) 039 * @see #setLongitude(int, int, double, String) 040 */ 041 private double longitude; 042 043 /** 044 * @see #getLocationName() 045 * @see #setLocationName(String) 046 */ 047 private String locationName; 048 049 /** 050 * @see #getTimeZone() 051 * @see #setTimeZone(TimeZone) 052 */ 053 private TimeZone timeZone; 054 055 /** 056 * @see #getElevation() 057 * @see #setElevation(double) 058 */ 059 private double elevation; 060 061 /** 062 * Constant for a distance type calculation. 063 * @see #getGeodesicDistance(GeoLocation) 064 */ 065 private static final int DISTANCE = 0; 066 067 /** 068 * Constant for a initial bearing type calculation. 069 * @see #getGeodesicInitialBearing(GeoLocation) 070 */ 071 private static final int INITIAL_BEARING = 1; 072 073 /** 074 * Constant for a final bearing type calculation. 075 * @see #getGeodesicFinalBearing(GeoLocation) 076 */ 077 private static final int FINAL_BEARING = 2; 078 079 080 /** 081 * Constant for a final bearing type calculation. 082 * @see #getGeodesicFinalBearing(GeoLocation) 083 */ 084 private static final int POINT = 3; 085 086 /** constant for milliseconds in a minute (60,000) */ 087 private static final long MINUTE_MILLIS = 60 * 1000; 088 089 /** constant for milliseconds in an hour (3,600,000) */ 090 private static final long HOUR_MILLIS = MINUTE_MILLIS * 60; 091 092 /** 093 * Method to get the elevation in Meters. 094 * 095 * @return Returns the elevation in Meters. 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, NaN or infinite. 106 */ 107 public void setElevation(double elevation) { 108 if (elevation < 0) { 109 throw new IllegalArgumentException("Elevation cannot be negative"); 110 } 111 if (Double.isNaN(elevation) || Double.isInfinite(elevation)) { 112 throw new IllegalArgumentException("Elevation must not be NaN or infinite"); 113 } 114 this.elevation = elevation; 115 } 116 117 /** 118 * GeoLocation constructor with parameters for all required fields. 119 * 120 * @param name 121 * The location name for display use such as "Lakewood, NJ" 122 * @param latitude 123 * the latitude in a double format such as 40.095965 for Lakewood, NJ. 124 * <b>Note:</b> For latitudes south of the equator, a negative value should be used. 125 * @param longitude 126 * double the longitude in a double format such as -74.222130 for Lakewood, NJ. 127 * <b>Note:</b> For longitudes east of the <a href="https://en.wikipedia.org/wiki/Prime_Meridian">Prime 128 * Meridian</a> (Greenwich), a negative value should be used. 129 * @param timeZone 130 * the <code>TimeZone</code> for the location. 131 */ 132 public GeoLocation(String name, double latitude, double longitude, TimeZone timeZone) { 133 this(name, latitude, longitude, 0, timeZone); 134 } 135 136 /** 137 * GeoLocation constructor with parameters for all required fields. 138 * 139 * @param name 140 * The location name for display use such as "Lakewood, NJ" 141 * @param latitude 142 * the latitude in a double format such as 40.095965 for Lakewood, NJ. 143 * <b>Note:</b> For latitudes south of the equator, a negative value should be used. 144 * @param longitude 145 * double the longitude in a double format such as -74.222130 for Lakewood, NJ. 146 * <b>Note:</b> For longitudes east of the <a href="https://en.wikipedia.org/wiki/Prime_Meridian">Prime 147 * Meridian</a> (Greenwich), a negative value should be used. 148 * @param elevation 149 * the elevation above sea level in Meters. Elevation is not used in most algorithms used for calculating 150 * sunrise and set. 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. 176 * 177 * @param latitude 178 * The degrees of latitude to set. The values should be between -90° and 90°. An 179 * IllegalArgumentException will be thrown if the value exceeds the limit. For example 40.095965 would be 180 * used for Lakewood, NJ. <b>Note:</b> For latitudes south of the equator, a negative value should be 181 * used. 182 */ 183 public void setLatitude(double latitude) { 184 if (latitude > 90 || latitude < -90) { 185 throw new IllegalArgumentException("Latitude must be between -90 and 90"); 186 } 187 this.latitude = latitude; 188 } 189 190 /** 191 * Method to set the latitude in degrees, minutes and seconds. 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> 198 * @param seconds 199 * <a href="https://en.wikipedia.org/wiki/Minute_of_arc#Cartography">seconds of arc</a> 200 * @param direction 201 * N for north and S for south. An IllegalArgumentException will be thrown if the value is not S or N. 202 */ 203 public void setLatitude(int degrees, int minutes, double seconds, String direction) { 204 double tempLat = degrees + ((minutes + (seconds / 60.0)) / 60.0); 205 if (tempLat > 90 || tempLat < 0) { //FIXME An exception should be thrown if degrees, minutes or seconds are negative 206 throw new IllegalArgumentException( 207 "Latitude must be between 0 and 90. Use direction of S instead of negative."); 208 } 209 if (direction.equals("S")) { 210 tempLat *= -1; 211 } else if (!direction.equals("N")) { 212 throw new IllegalArgumentException("Latitude direction must be N or S"); 213 } 214 this.latitude = tempLat; 215 } 216 217 /** 218 * @return Returns the latitude. 219 */ 220 public double getLatitude() { 221 return latitude; 222 } 223 224 /** 225 * Method to set the longitude in a double format. 226 * 227 * @param longitude 228 * The degrees of longitude to set in a double format between -180° and 180°. An 229 * IllegalArgumentException will be thrown if the value exceeds the limit. For example -74.2094 would be 230 * used for Lakewood, NJ. Note: for longitudes east of the <a 231 * href="https://en.wikipedia.org/wiki/Prime_Meridian">Prime Meridian</a> (Greenwich) a negative value 232 * should be used. 233 */ 234 public void setLongitude(double longitude) { 235 if (longitude > 180 || longitude < -180) { 236 throw new IllegalArgumentException("Longitude must be between -180 and 180"); 237 } 238 this.longitude = longitude; 239 } 240 241 /** 242 * Method to set the longitude in degrees, minutes and seconds. 243 * 244 * @param degrees 245 * The degrees of longitude to set between 0° and 180°. As an example 74 would be set for Lakewood, NJ. 246 * An IllegalArgumentException will be thrown if the value exceeds the limits. 247 * @param minutes 248 * <a href="https://en.wikipedia.org/wiki/Minute_of_arc#Cartography">minutes of arc</a> 249 * @param seconds 250 * <a href="https://en.wikipedia.org/wiki/Minute_of_arc#Cartography">seconds of arc</a> 251 * @param direction 252 * E for east of the <a href="https://en.wikipedia.org/wiki/Prime_Meridian">Prime Meridian</a> or W for west of it. 253 * An IllegalArgumentException will be thrown if 254 * the value is not E or W. 255 */ 256 public void setLongitude(int degrees, int minutes, double seconds, String direction) { 257 double longTemp = degrees + ((minutes + (seconds / 60.0)) / 60.0); 258 if (longTemp > 180 || this.longitude < 0) { //FIXME An exception should be thrown if degrees, minutes or seconds are negative 259 throw new IllegalArgumentException("Longitude must be between 0 and 180. Use a direction of W instead of negative."); 260 } 261 if (direction.equals("W")) { 262 longTemp *= -1; 263 } else if (!direction.equals("E")) { 264 throw new IllegalArgumentException("Longitude direction must be E or W"); 265 } 266 this.longitude = longTemp; 267 } 268 269 /** 270 * @return Returns the longitude. 271 */ 272 public double getLongitude() { 273 return longitude; 274 } 275 276 /** 277 * @return Returns the location name. 278 */ 279 public String getLocationName() { 280 return locationName; 281 } 282 283 /** 284 * @param name 285 * The setter method for the display name. 286 */ 287 public void setLocationName(String name) { 288 this.locationName = name; 289 } 290 291 /** 292 * @return Returns the timeZone. 293 */ 294 public TimeZone getTimeZone() { 295 return timeZone; 296 } 297 298 /** 299 * Method to set the TimeZone. If this is ever set after the GeoLocation is set in the 300 * {@link com.kosherjava.zmanim.AstronomicalCalendar}, it is critical that 301 * {@link com.kosherjava.zmanim.AstronomicalCalendar#getCalendar()}. 302 * {@link java.util.Calendar#setTimeZone(TimeZone) setTimeZone(TimeZone)} be called in order for the 303 * AstronomicalCalendar to output times in the expected offset. This situation will arise if the 304 * AstronomicalCalendar is ever {@link com.kosherjava.zmanim.AstronomicalCalendar#clone() cloned}. 305 * 306 * @param timeZone 307 * The timeZone to set. 308 */ 309 public void setTimeZone(TimeZone timeZone) { 310 this.timeZone = timeZone; 311 } 312 313 /** 314 * A method that will return the location's local mean time offset in milliseconds from local <a 315 * href="https://en.wikipedia.org/wiki/Standard_time">standard time</a>. The globe is split into 360°, with 316 * 15° per hour of the day. For a local that is at a longitude that is evenly divisible by 15 (longitude % 15 == 317 * 0), at solar {@link com.kosherjava.zmanim.AstronomicalCalendar#getSunTransit() noon} (with adjustment for the <a 318 * href="https://en.wikipedia.org/wiki/Equation_of_time">equation of time</a>) the sun should be directly overhead, 319 * so a user who is 1° west of this will have noon at 4 minutes after standard time noon, and conversely, a user 320 * who is 1° east of the 15° longitude will have noon at 11:56 AM. Lakewood, N.J., whose longitude is 321 * -74.2094, is 0.7906 away from the closest multiple of 15 at -75°. This is multiplied by 4 to yield 3 minutes 322 * and 10 seconds earlier than standard time. The offset returned does not account for the <a 323 * href="https://en.wikipedia.org/wiki/Daylight_saving_time">Daylight saving time</a> offset since this class is 324 * unaware of dates. 325 * 326 * @return the offset in milliseconds not accounting for Daylight saving time. A positive value will be returned 327 * East of the 15° timezone line, and a negative value West of it. 328 * @since 1.1 329 */ 330 public long getLocalMeanTimeOffset() { 331 return (long) (getLongitude() * 4 * MINUTE_MILLIS - getTimeZone().getRawOffset()); 332 } 333 334 /** 335 * Adjust the date for <a href="https://en.wikipedia.org/wiki/180th_meridian">antimeridian</a> crossover. This is 336 * needed to deal with edge cases such as Samoa that use a different calendar date than expected based on their 337 * geographic location. 338 * 339 * The actual Time Zone offset may deviate from the expected offset based on the longitude. Since the 'absolute time' 340 * calculations are always based on longitudinal offset from UTC for a given date, the date is presumed to only 341 * increase East of the Prime Meridian, and to only decrease West of it. For Time Zones that cross the antimeridian, 342 * the date will be artificially adjusted before calculation to conform with this presumption. 343 * 344 * For example, Apia, Samoa with a longitude of -171.75 uses a local offset of +14:00. When calculating sunrise for 345 * 2018-02-03, the calculator should operate using 2018-02-02 since the expected zone is -11. After determining the 346 * UTC time, the local DST offset of <a href="https://en.wikipedia.org/wiki/UTC%2B14:00">UTC+14:00</a> should be applied 347 * to bring the date back to 2018-02-03. 348 * 349 * @return the number of days to adjust the date This will typically be 0 unless the date crosses the antimeridian 350 */ 351 public int getAntimeridianAdjustment() { 352 double localHoursOffset = getLocalMeanTimeOffset() / (double)HOUR_MILLIS; 353 354 if (localHoursOffset >= 20){// if the offset is 20 hours or more in the future (never expected anywhere other 355 // than a location using a timezone across the anti meridian to the east such as Samoa) 356 return 1; // roll the date forward a day 357 } else if (localHoursOffset <= -20) { // if the offset is 20 hours or more in the past (no current location is known 358 //that crosses the antimeridian to the west, but better safe than sorry) 359 return -1; // roll the date back a day 360 } 361 return 0; //99.999% of the world will have no adjustment 362 } 363 364 /** 365 * Calculate the initial <a href="https://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing between this 366 * Object and a second Object passed to this method using <a 367 * href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, "<a 368 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 369 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975 370 * 371 * @param location 372 * the destination location 373 * @return the initial bearing 374 */ 375 public double getGeodesicInitialBearing(GeoLocation location) { 376 return vincentyInverseFormula(location, INITIAL_BEARING); 377 } 378 379 /** 380 * Calculate the final <a href="https://en.wikipedia.org/wiki/Great_circle">geodesic</a> bearing between this Object 381 * and a second Object passed to this method using <a href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus 382 * Vincenty's</a> inverse formula See T Vincenty, "<a href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and 383 * Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations</a>", Survey Review, vol 384 * XXII no 176, 1975 385 * 386 * @param location 387 * the destination location 388 * @return the final bearing 389 */ 390 public double getGeodesicFinalBearing(GeoLocation location) { 391 return vincentyInverseFormula(location, FINAL_BEARING); 392 } 393 394 /** 395 * Calculate <a href="https://en.wikipedia.org/wiki/Great-circle_distance">geodesic distance</a> in Meters between 396 * this Object and a second Object passed to this method using <a 397 * href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, "<a 398 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 399 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975 400 * 401 * @see #vincentyInverseFormula(GeoLocation, int) 402 * @param location 403 * the destination location 404 * @return the geodesic distance in Meters 405 */ 406 public double getGeodesicDistance(GeoLocation location) { 407 return vincentyInverseFormula(location, DISTANCE); 408 } 409 410 /** 411 * Calculate the destination point based on an initial bearing and distance in meters from the current location using 412 * <a href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> direct formula. See T Vincenty, "<a 413 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 414 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975. 415 * 416 * @param initialBearing the initialBearing 417 * @param distance the distance in meters. 418 * @return the GeoLocation containing the destination point. The Timezone is set to the origin point (current object). 419 */ 420 private GeoLocation vincentyDirectFormulaDestination(double initialBearing, double distance) { 421 double major_semi_axis = 6378137; 422 double minor_semi_axis = 6356752.3142; 423 double flattening = 1 / 298.257223563; // WGS-84 ellipsiod 424 double initial_bearing_radians = Math.toRadians(initialBearing); 425 double sinAzimuth1 = Math.sin(initial_bearing_radians); 426 double cosAzimuth1 = Math.cos(initial_bearing_radians); 427 double tanU1 = (1 - flattening) * Math.tan(Math.toRadians(getLatitude())); 428 double cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)); 429 double sinU1 = tanU1 * cosU1; 430 double eq_p1_ang_dist = Math.atan2(tanU1, cosAzimuth1); // eq_p1_ang_dist = angular distance on the sphere from the equator to P1 431 double sinAzimuth = cosU1 * sinAzimuth1; //azimuth of the geodesic at the equator 432 double cosSqAzimuth = 1 - sinAzimuth*sinAzimuth; 433 double uSq = cosSqAzimuth * (Math.pow(major_semi_axis, 2) - Math.pow(minor_semi_axis, 2) / Math.pow(minor_semi_axis, 2)); 434 double a = 1 + uSq/16384*(4096 + uSq *(-768 + uSq *(320 - 175 * uSq))); 435 double b = uSq / 1024 * (256 + uSq *(-128 + uSq * (74-47 * uSq))); 436 double p1_p2_ang_dist = distance / (minor_semi_axis * a); //p1_p2_ang_dist = angular distance P1 P2 on the sphere 437 double sinSigma = Double.NaN; 438 double cosSigma = Double.NaN; 439 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 440 double a_prime = Double.NaN; 441 int iterations = 0; 442 443 do { 444 cos2_eq_mid_ang_distance = Math.cos(2*eq_p1_ang_dist + p1_p2_ang_dist); 445 sinSigma = Math.sin(p1_p2_ang_dist); 446 cosSigma = Math.cos(p1_p2_ang_dist); 447 double delta_ang_distance = b * sinSigma * (cos2_eq_mid_ang_distance + b / 4 * 448 (cosSigma * (-1 + 2 * cos2_eq_mid_ang_distance * cos2_eq_mid_ang_distance) - b / 6 * cos2_eq_mid_ang_distance * 449 (-3+4*sinSigma*sinSigma)*(-3+4*cos2_eq_mid_ang_distance*cos2_eq_mid_ang_distance))); 450 a_prime = p1_p2_ang_dist; 451 p1_p2_ang_dist = distance / (minor_semi_axis * a) + delta_ang_distance; 452 } while (Math.abs(p1_p2_ang_dist-a_prime) > 1e-12 && ++iterations < 100); // iterate until negligible change in lambda (about 0.006mm) 453 454 double x = sinU1 * sinSigma - cosU1 * cosSigma * cosAzimuth1; 455 double other_latitude = Math.toDegrees(Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAzimuth1, (1 - flattening) * Math.sqrt(sinAzimuth*sinAzimuth + x * x))); 456 double lambda = Math.atan2(sinSigma*sinAzimuth1, cosU1*cosSigma - sinU1*sinSigma*cosAzimuth1); 457 double c = flattening/16*cosSqAzimuth*(4+flattening*(4-3*cosSqAzimuth)); 458 double l = lambda - (1-c) * flattening * sinAzimuth * 459 (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))); 460 double other_longitude = longitude + Math.toDegrees(l); 461 462 return new GeoLocation("Destination", other_latitude, other_longitude, getTimeZone()); //ToDo - we can easily return final_bearing, it just needs some minor refactoring 463 } 464 465 public GeoLocation findLocationAtTime(GeoLocation destination, long start_time, long end_time, long target_time){ 466 double total_distance = vincentyInverseFormula(destination, DISTANCE); 467 double initial_bearing = vincentyInverseFormula(destination, INITIAL_BEARING); 468 double total_elapsed_time = end_time - start_time; 469 double time_difference = target_time - start_time; 470 double distance_travelled = total_distance * time_difference / total_elapsed_time; 471 return vincentyDirectFormulaDestination(initial_bearing, distance_travelled); 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 ellipsiod 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.0828&deg</Latitude> 611 * <Longitude>-74.2094&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 StringBuilder sb = new StringBuilder(); 624 sb.append("<GeoLocation>\n"); 625 sb.append("\t<LocationName>").append(getLocationName()).append("</LocationName>\n"); 626 sb.append("\t<Latitude>").append(getLatitude()).append("</Latitude>\n"); 627 sb.append("\t<Longitude>").append(getLongitude()).append("</Longitude>\n"); 628 sb.append("\t<Elevation>").append(getElevation()).append(" Meters").append("</Elevation>\n"); 629 sb.append("\t<TimezoneName>").append(getTimeZone().getID()).append("</TimezoneName>\n"); 630 sb.append("\t<TimeZoneDisplayName>").append(getTimeZone().getDisplayName()).append("</TimeZoneDisplayName>\n"); 631 sb.append("\t<TimezoneGMTOffset>").append(getTimeZone().getRawOffset() / HOUR_MILLIS) 632 .append("</TimezoneGMTOffset>\n"); 633 sb.append("\t<TimezoneDSTOffset>").append(getTimeZone().getDSTSavings() / HOUR_MILLIS) 634 .append("</TimezoneDSTOffset>\n"); 635 sb.append("</GeoLocation>"); 636 return sb.toString(); 637 } 638 639 /** 640 * @see java.lang.Object#equals(Object) 641 */ 642 public boolean equals(Object object) { 643 if (this == object) 644 return true; 645 if (!(object instanceof GeoLocation)) 646 return false; 647 GeoLocation geo = (GeoLocation) object; 648 return Double.doubleToLongBits(this.latitude) == Double.doubleToLongBits(geo.latitude) 649 && Double.doubleToLongBits(this.longitude) == Double.doubleToLongBits(geo.longitude) 650 && this.elevation == geo.elevation 651 && (this.locationName == null ? geo.locationName == null : this.locationName.equals(geo.locationName)) 652 && (this.timeZone == null ? geo.timeZone == null : this.timeZone.equals(geo.timeZone)); 653 } 654 655 /** 656 * @see java.lang.Object#hashCode() 657 */ 658 public int hashCode() { 659 660 int result = 17; 661 long latLong = Double.doubleToLongBits(this.latitude); 662 long lonLong = Double.doubleToLongBits(this.longitude); 663 long elevLong = Double.doubleToLongBits(this.elevation); 664 int latInt = (int) (latLong ^ (latLong >>> 32)); 665 int lonInt = (int) (lonLong ^ (lonLong >>> 32)); 666 int elevInt = (int) (elevLong ^ (elevLong >>> 32)); 667 result = 37 * result + getClass().hashCode(); 668 result += 37 * result + latInt; 669 result += 37 * result + lonInt; 670 result += 37 * result + elevInt; 671 result += 37 * result + (this.locationName == null ? 0 : this.locationName.hashCode()); 672 result += 37 * result + (this.timeZone == null ? 0 : this.timeZone.hashCode()); 673 return result; 674 } 675 676 /** 677 * @see java.lang.Object#toString() 678 */ 679 public String toString() { 680 StringBuilder sb = new StringBuilder(); 681 sb.append("\nLocation Name:\t\t\t").append(getLocationName()); 682 sb.append("\nLatitude:\t\t\t").append(getLatitude()).append("\u00B0"); 683 sb.append("\nLongitude:\t\t\t").append(getLongitude()).append("\u00B0"); 684 sb.append("\nElevation:\t\t\t").append(getElevation()).append(" Meters"); 685 sb.append("\nTimezone ID:\t\t\t").append(getTimeZone().getID()); 686 sb.append("\nTimezone Display Name:\t\t").append(getTimeZone().getDisplayName()) 687 .append(" (").append(getTimeZone().getDisplayName(false, TimeZone.SHORT)).append(")"); 688 sb.append("\nTimezone GMT Offset:\t\t").append(getTimeZone().getRawOffset() / HOUR_MILLIS); 689 sb.append("\nTimezone DST Offset:\t\t").append(getTimeZone().getDSTSavings() / HOUR_MILLIS); 690 return sb.toString(); 691 } 692 693 /** 694 * An implementation of the {@link java.lang.Object#clone()} method that creates a <a 695 * href="https://en.wikipedia.org/wiki/Object_copy#Deep_copy">deep copy</a> of the object. 696 * <b>Note:</b> If the {@link java.util.TimeZone} in the clone will be changed from the original, it is critical 697 * that {@link com.kosherjava.zmanim.AstronomicalCalendar#getCalendar()}. 698 * {@link java.util.Calendar#setTimeZone(TimeZone) setTimeZone(TimeZone)} is called after cloning in order for the 699 * AstronomicalCalendar to output times in the expected offset. 700 * 701 * @see java.lang.Object#clone() 702 * @since 1.1 703 */ 704 public Object clone() { 705 GeoLocation clone = null; 706 try { 707 clone = (GeoLocation) super.clone(); 708 } catch (CloneNotSupportedException cnse) { 709 //Required by the compiler. Should never be reached since we implement clone() 710 } 711 clone.timeZone = (TimeZone) getTimeZone().clone(); 712 clone.locationName = getLocationName(); 713 return clone; 714 } 715}