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