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