001/* 002 * Zmanim Java API 003 * Copyright (C) 2004-2011 Eliyahu Hershfeld 004 * 005 * This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General 006 * Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) 007 * any later version. 008 * 009 * This library is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; without even the implied 010 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 011 * details. 012 * You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to 013 * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA, 014 * or connect to: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html 015 */ 016package net.sourceforge.zmanim.util; 017 018import java.util.Calendar; 019 020/** 021 * Implementation of sunrise and sunset methods to calculate astronomical times. This calculator uses the Java algorithm 022 * written by <a href="http://www.jstot.me.uk/jsuntimes/">Jonathan Stott</a> that is based on the implementation by <a 023 * href=""http://noaa.gov">NOAA - National Oceanic and Atmospheric Administration</a>'s <a href = 024 * "http://www.srrb.noaa.gov/highlights/sunrise/sunrisehtml">Surface Radiation Research Branch</a>. NOAA's <a 025 * href="http://www.srrb.noaa.gov/highlights/sunrise/solareqns.PDF">implementation</a> is based on equations from <a 026 * href="http://www.willbell.com/math/mc1.htm">Astronomical Algorithms</a> by <a 027 * href="http://en.wikipedia.org/wiki/Jean_Meeus">Jean Meeus</a>. Jonathan's implementation was released under the GPL. 028 * Added to the algorithm is an adjustment of the zenith to account for elevation. 029 * 030 * @author Jonathan Stott 2000 - 2004 031 * @author © Eliyahu Hershfeld 2004 - 2011 032 * @version 1.1 033 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA algorithm 034 * JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the 035 * NOAACalculator class. 036 * @see net.sourceforge.zmanim.util.NOAACalculator 037 */ 038public class JSuntimeCalculator extends AstronomicalCalculator { 039 private String calculatorName = "US National Oceanic and Atmospheric Administration Algorithm"; 040 041 /** 042 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA 043 * algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been 044 * replaced by the NOAACalculator class. 045 * @see net.sourceforge.zmanim.util.NOAACalculator#getCalculatorName() 046 */ 047 public String getCalculatorName() { 048 return this.calculatorName; 049 } 050 051 /** 052 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA 053 * algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been 054 * replaced by the NOAACalculator class. 055 * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean) 056 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean) 057 * @throws ZmanimException 058 * if the year entered == 2000. This calculator can't properly deal with the year 2000. It can properly 059 * calculate times for years <> 2000. 060 */ 061 public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 062 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 063 double adjustedZenith = adjustZenith(zenith, elevation); 064 double timeMins = morningPhenomenon(dateToJulian(calendar), geoLocation.getLatitude(), 065 -geoLocation.getLongitude(), adjustedZenith); 066 return timeMins / 60; 067 } 068 069 /** 070 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAAA 071 * algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been 072 * replaced by the NOAACalculator class. 073 * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunset(Calendar, GeoLocation, double, boolean) 074 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean) 075 * @throws ZmanimException 076 * if the year entered == 2000. This calculator can't properly deal with the year 2000. It can properly 077 * calculate times for years <> 2000. 078 */ 079 public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 080 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 081 double adjustedZenith = adjustZenith(zenith, elevation); 082 double timeMins = eveningPhenomenon(dateToJulian(calendar), geoLocation.getLatitude(), 083 -geoLocation.getLongitude(), adjustedZenith); 084 return timeMins / 60; 085 } 086 087 /** 088 * Calculate the UTC of a morning phenomenon for the given day at the given latitude and longitude on Earth 089 * 090 * @param julian 091 * Julian day 092 * @param latitude 093 * latitude of observer in degrees 094 * @param longitude 095 * longitude of observer in degrees 096 * @param zenithDistance 097 * one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE, Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE, 098 * Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE, Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE. 099 * @return time in minutes from zero Z 100 */ 101 private static double morningPhenomenon(double julian, double latitude, double longitude, double zenithDistance) { 102 double t = julianDayToJulianCenturies(julian); 103 double eqtime = equationOfTime(t); 104 double solarDec = sunDeclination(t); 105 double hourangle = hourAngleMorning(latitude, solarDec, zenithDistance); 106 double delta = longitude - Math.toDegrees(hourangle); 107 double timeDiff = 4 * delta; 108 double timeUTC = 720 + timeDiff - eqtime; 109 110 // Second pass includes fractional julian day in gamma calc 111 double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t) + timeUTC / 1440); 112 eqtime = equationOfTime(newt); 113 solarDec = sunDeclination(newt); 114 hourangle = hourAngleMorning(latitude, solarDec, zenithDistance); 115 delta = longitude - Math.toDegrees(hourangle); 116 timeDiff = 4 * delta; 117 118 double morning = 720 + timeDiff - eqtime; 119 return morning; 120 } 121 122 /** 123 * Calculate the UTC of an evening phenomenon for the given day at the given latitude and longitude on Earth 124 * 125 * @param julian 126 * Julian day 127 * @param latitude 128 * latitude of observer in degrees 129 * @param longitude 130 * longitude of observer in degrees 131 * @param zenithDistance 132 * one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE, Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE, 133 * Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE, Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE. 134 * @return time in minutes from zero Z 135 */ 136 private static double eveningPhenomenon(double julian, double latitude, double longitude, double zenithDistance) { 137 double t = julianDayToJulianCenturies(julian); 138 139 // First calculates sunrise and approx length of day 140 double eqtime = equationOfTime(t); 141 double solarDec = sunDeclination(t); 142 double hourangle = hourAngleEvening(latitude, solarDec, zenithDistance); 143 144 double delta = longitude - Math.toDegrees(hourangle); 145 double timeDiff = 4 * delta; 146 double timeUTC = 720 + timeDiff - eqtime; 147 148 // first pass used to include fractional day in gamma calc 149 double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t) + timeUTC / 1440); 150 eqtime = equationOfTime(newt); 151 solarDec = sunDeclination(newt); 152 hourangle = hourAngleEvening(latitude, solarDec, zenithDistance); 153 154 delta = longitude - Math.toDegrees(hourangle); 155 timeDiff = 4 * delta; 156 157 double evening = 720 + timeDiff - eqtime; 158 return evening; 159 } 160 161 private static double dateToJulian(Calendar date) { 162 int year = date.get(Calendar.YEAR); 163 int month = date.get(Calendar.MONTH) + 1; 164 int day = date.get(Calendar.DAY_OF_MONTH); 165 int hour = date.get(Calendar.HOUR_OF_DAY); 166 int minute = date.get(Calendar.MINUTE); 167 int second = date.get(Calendar.SECOND); 168 169 double extra = (100.0 * year) + month - 190002.5; 170 double JD = (367.0 * year) - (Math.floor(7.0 * (year + Math.floor((month + 9.0) / 12.0)) / 4.0)) 171 + Math.floor((275.0 * month) / 9.0) + day + ((hour + ((minute + (second / 60.0)) / 60.0)) / 24.0) 172 + 1721013.5 - ((0.5 * extra) / Math.abs(extra)) + 0.5; 173 return JD; 174 } 175 176 /** 177 * Convert Julian Day to centuries since J2000.0 178 * 179 * @param julian 180 * The Julian Day to convert 181 * @return the value corresponding to the Julian Day 182 */ 183 private static double julianDayToJulianCenturies(double julian) { 184 return (julian - 2451545) / 36525; 185 } 186 187 /** 188 * Convert centuries since J2000.0 to Julian Day 189 * 190 * @param t 191 * Number of Julian centuries since J2000.0 192 * @return The Julian Day corresponding to the value of t 193 */ 194 private static double julianCenturiesToJulianDay(double t) { 195 return (t * 36525) + 2451545; 196 } 197 198 /** 199 * Calculate the difference between true solar time and mean solar time 200 * 201 * @param t 202 * Number of Julian centuries since J2000.0 203 * @return 204 */ 205 private static double equationOfTime(double t) { 206 double epsilon = obliquityCorrection(t); 207 double l0 = geomMeanLongSun(t); 208 double e = eccentricityOfEarthsOrbit(t); 209 double m = geometricMeanAnomalyOfSun(t); 210 double y = Math.pow((Math.tan(Math.toRadians(epsilon) / 2)), 2); 211 212 double eTime = y * Math.sin(2 * Math.toRadians(l0)) - 2 * e * Math.sin(Math.toRadians(m)) + 4 * e * y 213 * Math.sin(Math.toRadians(m)) * Math.cos(2 * Math.toRadians(l0)) - 0.5 * y * y 214 * Math.sin(4 * Math.toRadians(l0)) - 1.25 * e * e * Math.sin(2 * Math.toRadians(m)); 215 return Math.toDegrees(eTime) * 4; 216 } 217 218 /** 219 * Calculate the declination of the sun 220 * 221 * @param t 222 * Number of Julian centuries since J2000.0 223 * @return The Sun's declination in degrees 224 */ 225 private static double sunDeclination(double t) { 226 double e = obliquityCorrection(t); 227 double lambda = sunsApparentLongitude(t); 228 229 double sint = Math.sin(Math.toRadians(e)) * Math.sin(Math.toRadians(lambda)); 230 return Math.toDegrees(Math.asin(sint)); 231 } 232 233 /** 234 * calculate the hour angle of the sun for a morning phenomenon for the given latitude 235 * 236 * @param lat 237 * Latitude of the observer in degrees 238 * @param solarDec 239 * declination of the sun in degrees 240 * @param zenithDistance 241 * zenith distance of the sun in degrees 242 * @return hour angle of sunrise in radians 243 */ 244 private static double hourAngleMorning(double lat, double solarDec, double zenithDistance) { 245 return (Math.acos(Math.cos(Math.toRadians(zenithDistance)) 246 / (Math.cos(Math.toRadians(lat)) * Math.cos(Math.toRadians(solarDec))) - Math.tan(Math.toRadians(lat)) 247 * Math.tan(Math.toRadians(solarDec)))); 248 } 249 250 /** 251 * Calculate the hour angle of the sun for an evening phenomenon for the given latitude 252 * 253 * @param lat 254 * Latitude of the observer in degrees 255 * @param solarDec 256 * declination of the Sun in degrees 257 * @param zenithDistance 258 * zenith distance of the sun in degrees 259 * @return hour angle of sunset in radians 260 */ 261 private static double hourAngleEvening(double lat, double solarDec, double zenithDistance) { 262 return -hourAngleMorning(lat, solarDec, zenithDistance); 263 } 264 265 /** 266 * Calculate the corrected obliquity of the ecliptic 267 * 268 * @param t 269 * Number of Julian centuries since J2000.0 270 * @return Corrected obliquity in degrees 271 */ 272 private static double obliquityCorrection(double t) { 273 return meanObliquityOfEcliptic(t) + 0.00256 * Math.cos(Math.toRadians(125.04 - 1934.136 * t)); 274 } 275 276 /** 277 * Calculate the mean obliquity of the ecliptic 278 * 279 * @param t 280 * Number of Julian centuries since J2000.0 281 * @return Mean obliquity in degrees 282 */ 283 private static double meanObliquityOfEcliptic(double t) { 284 return 23 + (26 + (21.448 - t * (46.815 + t * (0.00059 - t * (0.001813))) / 60)) / 60; 285 } 286 287 /** 288 * Calculate the geometric mean longitude of the sun 289 * 290 * @param t 291 * number of Julian centuries since J2000.0 292 * @return the geometric mean longitude of the sun in degrees 293 */ 294 private static double geomMeanLongSun(double t) { 295 double l0 = 280.46646 + t * (36000.76983 + 0.0003032 * t); 296 297 while ((l0 >= 0) && (l0 <= 360)) { 298 if (l0 > 360) { 299 l0 = l0 - 360; 300 } 301 302 if (l0 < 0) { 303 l0 = l0 + 360; 304 } 305 } 306 return l0; 307 } 308 309 /** 310 * Calculate the eccentricity of Earth's orbit 311 * 312 * @param t 313 * Number of Julian centuries since J2000.0 314 * @return the eccentricity 315 */ 316 private static double eccentricityOfEarthsOrbit(double t) { 317 return 0.016708634 - t * (0.000042037 + 0.0000001267 * t); 318 } 319 320 /** 321 * Calculate the geometric mean anomaly of the Sun 322 * 323 * @param t 324 * Number of Julian centuries since J2000.0 325 * @return the geometric mean anomaly of the Sun in degrees 326 */ 327 private static double geometricMeanAnomalyOfSun(double t) { 328 return 357.52911 + t * (35999.05029 - 0.0001537 * t); 329 } 330 331 /** 332 * Calculate the apparent longitude of the sun 333 * 334 * @param t 335 * Number of Julian centuries since J2000.0 336 * @return The apparent longitude of the Sun in degrees 337 */ 338 private static double sunsApparentLongitude(double t) { 339 return sunsTrueLongitude(t) - 0.00569 - 0.00478 * Math.sin(Math.toRadians(125.04 - 1934.136 * t)); 340 } 341 342 /** 343 * Calculate the true longitude of the sun 344 * 345 * @param t 346 * Number of Julian centuries since J2000.0 347 * @return The Sun's true longitude in degrees 348 */ 349 private static double sunsTrueLongitude(double t) { 350 return geomMeanLongSun(t) + equationOfCentreForSun(t); 351 } 352 353 /** 354 * Calculate the equation of centre for the Sun 355 * 356 * @param centuries 357 * Number of Julian centuries since J2000.0 358 * @return The equation of centre for the Sun in degrees 359 */ 360 private static double equationOfCentreForSun(double t) { 361 double m = geometricMeanAnomalyOfSun(t); 362 363 return Math.sin(Math.toRadians(m)) * (1.914602 - t * (0.004817 + 0.000014 * t)) 364 + Math.sin(2 * Math.toRadians(m)) * (0.019993 - 0.000101 * t) + Math.sin(3 * Math.toRadians(m)) 365 * 0.000289; 366 } 367}