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