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 /** 466 * Returns the location along a great circle between the current GeoLocation instance and a second one assuming a start, end 467 * and target time. 468 * @param destination the destination location 469 * @param startTime the start time 470 * @param endTime the end time 471 * @param targetTime that you are looking for a location for 472 * @return the GeoLocation at the target time 473 */ 474 @Deprecated // (forRemoval=true) // add back once Java 9 is the minimum supported version 475 public GeoLocation findLocationAtTime(GeoLocation destination, long startTime, long endTime, long targetTime){ 476 double total_distance = vincentyInverseFormula(destination, DISTANCE); 477 double initial_bearing = vincentyInverseFormula(destination, INITIAL_BEARING); 478 double total_elapsed_time = endTime - startTime; 479 double time_difference = targetTime - startTime; 480 double distance_travelled = total_distance * time_difference / total_elapsed_time; 481 return vincentyDirectFormulaDestination(initial_bearing, distance_travelled); 482 } 483 484 /** 485 * Calculate <a href="https://en.wikipedia.org/wiki/Great-circle_distance">geodesic distance</a> in Meters between 486 * this Object and a second Object passed to this method using <a 487 * href="https://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus Vincenty's</a> inverse formula See T Vincenty, "<a 488 * href="https://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf">Direct and Inverse Solutions of Geodesics on the Ellipsoid 489 * with application of nested equations</a>", Survey Review, vol XXII no 176, 1975 490 * 491 * @param location 492 * the destination location 493 * @param formula 494 * This formula calculates initial bearing ({@link #INITIAL_BEARING}), final bearing ( 495 * {@link #FINAL_BEARING}) and distance ({@link #DISTANCE}). 496 * @return geodesic distance in Meters 497 */ 498 private double vincentyInverseFormula(GeoLocation location, int formula) { 499 double majorSemiAxis = 6378137; 500 double minorSemiAxis = 6356752.3142; 501 double f = 1 / 298.257223563; // WGS-84 ellipsiod 502 double L = Math.toRadians(location.getLongitude() - getLongitude()); 503 double U1 = Math.atan((1 - f) * Math.tan(Math.toRadians(getLatitude()))); 504 double U2 = Math.atan((1 - f) * Math.tan(Math.toRadians(location.getLatitude()))); 505 double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1); 506 double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2); 507 508 double lambda = L; 509 double lambdaP = 2 * Math.PI; 510 double iterLimit = 20; 511 double sinLambda = 0; 512 double cosLambda = 0; 513 double sinSigma = 0; 514 double cosSigma = 0; 515 double sigma = 0; 516 double sinAlpha = 0; 517 double cosSqAlpha = 0; 518 double cos2SigmaM = 0; 519 double C; 520 while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) { 521 sinLambda = Math.sin(lambda); 522 cosLambda = Math.cos(lambda); 523 sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) 524 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)); 525 if (sinSigma == 0) 526 return 0; // co-incident points 527 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda; 528 sigma = Math.atan2(sinSigma, cosSigma); 529 sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; 530 cosSqAlpha = 1 - sinAlpha * sinAlpha; 531 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha; 532 if (Double.isNaN(cos2SigmaM)) 533 cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6) 534 C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)); 535 lambdaP = lambda; 536 lambda = L + (1 - C) * f * sinAlpha 537 * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM))); 538 } 539 if (iterLimit == 0) 540 return Double.NaN; // formula failed to converge 541 542 double uSq = cosSqAlpha * (majorSemiAxis * majorSemiAxis - minorSemiAxis * minorSemiAxis) / (minorSemiAxis * minorSemiAxis); 543 double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))); 544 double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); 545 double deltaSigma = B 546 * sinSigma 547 * (cos2SigmaM + B 548 / 4 549 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM 550 * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM))); 551 double distance = minorSemiAxis * A * (sigma - deltaSigma); 552 553 // initial bearing 554 double fwdAz = Math.toDegrees(Math.atan2(cosU2 * sinLambda, cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)); 555 // final bearing 556 double revAz = Math.toDegrees(Math.atan2(cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda)); 557 if (formula == DISTANCE) { 558 return distance; 559 } else if (formula == INITIAL_BEARING) { 560 return fwdAz; 561 } else if (formula == FINAL_BEARING) { 562 return revAz; 563 } else { // should never happen 564 return Double.NaN; 565 } 566 } 567 568 /** 569 * Returns the <a href="https://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a> bearing from the current location to 570 * the GeoLocation passed in. 571 * 572 * @param location 573 * destination location 574 * @return the bearing in degrees 575 */ 576 public double getRhumbLineBearing(GeoLocation location) { 577 double dLon = Math.toRadians(location.getLongitude() - getLongitude()); 578 double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4) 579 / Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4)); 580 if (Math.abs(dLon) > Math.PI) 581 dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon); 582 return Math.toDegrees(Math.atan2(dLon, dPhi)); 583 } 584 585 /** 586 * Returns the <a href="https://en.wikipedia.org/wiki/Rhumb_line">rhumb line</a> distance from the current location 587 * to the GeoLocation passed in. 588 * 589 * @param location 590 * the destination location 591 * @return the distance in Meters 592 */ 593 public double getRhumbLineDistance(GeoLocation location) { 594 double earthRadius = 6378137; // Earth's radius in meters (WGS-84) 595 double dLat = Math.toRadians(location.getLatitude()) - Math.toRadians(getLatitude()); 596 double dLon = Math.abs(Math.toRadians(location.getLongitude()) - Math.toRadians(getLongitude())); 597 double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4) 598 / Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4)); 599 double q = dLat / dPhi; 600 601 if (!(Math.abs(q) <= Double.MAX_VALUE)) { 602 q = Math.cos(Math.toRadians(getLatitude())); 603 } 604 // if dLon over 180° take shorter rhumb across 180° meridian: 605 if (dLon > Math.PI) { 606 dLon = 2 * Math.PI - dLon; 607 } 608 double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon); 609 return d * earthRadius; 610 } 611 612 /** 613 * A method that returns an XML formatted <code>String</code> representing the serialized <code>Object</code>. Very 614 * similar to the toString method but the return value is in an xml format. The format currently used (subject to 615 * change) is: 616 * 617 * <pre> 618 * <GeoLocation> 619 * <LocationName>Lakewood, NJ</LocationName> 620 * <Latitude>40.0828&deg</Latitude> 621 * <Longitude>-74.2094&deg</Longitude> 622 * <Elevation>0 Meters</Elevation> 623 * <TimezoneName>America/New_York</TimezoneName> 624 * <TimeZoneDisplayName>Eastern Standard Time</TimeZoneDisplayName> 625 * <TimezoneGMTOffset>-5</TimezoneGMTOffset> 626 * <TimezoneDSTOffset>1</TimezoneDSTOffset> 627 * </GeoLocation> 628 * </pre> 629 * 630 * @return The XML formatted <code>String</code>. 631 */ 632 public String toXML() { 633 StringBuilder sb = new StringBuilder(); 634 sb.append("<GeoLocation>\n"); 635 sb.append("\t<LocationName>").append(getLocationName()).append("</LocationName>\n"); 636 sb.append("\t<Latitude>").append(getLatitude()).append("</Latitude>\n"); 637 sb.append("\t<Longitude>").append(getLongitude()).append("</Longitude>\n"); 638 sb.append("\t<Elevation>").append(getElevation()).append(" Meters").append("</Elevation>\n"); 639 sb.append("\t<TimezoneName>").append(getTimeZone().getID()).append("</TimezoneName>\n"); 640 sb.append("\t<TimeZoneDisplayName>").append(getTimeZone().getDisplayName()).append("</TimeZoneDisplayName>\n"); 641 sb.append("\t<TimezoneGMTOffset>").append(getTimeZone().getRawOffset() / HOUR_MILLIS) 642 .append("</TimezoneGMTOffset>\n"); 643 sb.append("\t<TimezoneDSTOffset>").append(getTimeZone().getDSTSavings() / HOUR_MILLIS) 644 .append("</TimezoneDSTOffset>\n"); 645 sb.append("</GeoLocation>"); 646 return sb.toString(); 647 } 648 649 /** 650 * @see java.lang.Object#equals(Object) 651 */ 652 public boolean equals(Object object) { 653 if (this == object) 654 return true; 655 if (!(object instanceof GeoLocation)) 656 return false; 657 GeoLocation geo = (GeoLocation) object; 658 return Double.doubleToLongBits(this.latitude) == Double.doubleToLongBits(geo.latitude) 659 && Double.doubleToLongBits(this.longitude) == Double.doubleToLongBits(geo.longitude) 660 && this.elevation == geo.elevation 661 && (this.locationName == null ? geo.locationName == null : this.locationName.equals(geo.locationName)) 662 && (this.timeZone == null ? geo.timeZone == null : this.timeZone.equals(geo.timeZone)); 663 } 664 665 /** 666 * @see java.lang.Object#hashCode() 667 */ 668 public int hashCode() { 669 670 int result = 17; 671 long latLong = Double.doubleToLongBits(this.latitude); 672 long lonLong = Double.doubleToLongBits(this.longitude); 673 long elevLong = Double.doubleToLongBits(this.elevation); 674 int latInt = (int) (latLong ^ (latLong >>> 32)); 675 int lonInt = (int) (lonLong ^ (lonLong >>> 32)); 676 int elevInt = (int) (elevLong ^ (elevLong >>> 32)); 677 result = 37 * result + getClass().hashCode(); 678 result += 37 * result + latInt; 679 result += 37 * result + lonInt; 680 result += 37 * result + elevInt; 681 result += 37 * result + (this.locationName == null ? 0 : this.locationName.hashCode()); 682 result += 37 * result + (this.timeZone == null ? 0 : this.timeZone.hashCode()); 683 return result; 684 } 685 686 /** 687 * @see java.lang.Object#toString() 688 */ 689 public String toString() { 690 StringBuilder sb = new StringBuilder(); 691 sb.append("\nLocation Name:\t\t\t").append(getLocationName()); 692 sb.append("\nLatitude:\t\t\t").append(getLatitude()).append("\u00B0"); 693 sb.append("\nLongitude:\t\t\t").append(getLongitude()).append("\u00B0"); 694 sb.append("\nElevation:\t\t\t").append(getElevation()).append(" Meters"); 695 sb.append("\nTimezone ID:\t\t\t").append(getTimeZone().getID()); 696 sb.append("\nTimezone Display Name:\t\t").append(getTimeZone().getDisplayName()) 697 .append(" (").append(getTimeZone().getDisplayName(false, TimeZone.SHORT)).append(")"); 698 sb.append("\nTimezone GMT Offset:\t\t").append(getTimeZone().getRawOffset() / HOUR_MILLIS); 699 sb.append("\nTimezone DST Offset:\t\t").append(getTimeZone().getDSTSavings() / HOUR_MILLIS); 700 return sb.toString(); 701 } 702 703 /** 704 * An implementation of the {@link java.lang.Object#clone()} method that creates a <a 705 * href="https://en.wikipedia.org/wiki/Object_copy#Deep_copy">deep copy</a> of the object. 706 * <b>Note:</b> If the {@link java.util.TimeZone} in the clone will be changed from the original, it is critical 707 * that {@link com.kosherjava.zmanim.AstronomicalCalendar#getCalendar()}. 708 * {@link java.util.Calendar#setTimeZone(TimeZone) setTimeZone(TimeZone)} is called after cloning in order for the 709 * AstronomicalCalendar to output times in the expected offset. 710 * 711 * @see java.lang.Object#clone() 712 * @since 1.1 713 */ 714 public Object clone() { 715 GeoLocation clone = null; 716 try { 717 clone = (GeoLocation) super.clone(); 718 } catch (CloneNotSupportedException cnse) { 719 //Required by the compiler. Should never be reached since we implement clone() 720 } 721 clone.timeZone = (TimeZone) getTimeZone().clone(); 722 clone.locationName = getLocationName(); 723 return clone; 724 } 725}