001 /*
002 * Zmanim Java API
003 * Copyright (C) 2004-2007 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 - 2007
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 /**
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 "US National Oceanic and Atmospheric Administration Algorithm";
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 //if(timeMins/60 < 4) {System.out.println("zenith="+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 //if(timeMins/60 > 28) {System.out.println("zenith="+zenith);}
110 return timeMins / 60;
111 }
112
113 /**
114 * Calculate the UTC of a morning phenomenon for the given day at the given
115 * latitude and longitude on Earth
116 *
117 * @param julian
118 * Julian day
119 * @param latitude
120 * latitude of observer in degrees
121 * @param longitude
122 * longitude of observer in degrees
123 * @param zenithDistance
124 * one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE,
125 * Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE,
126 * Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE,
127 * Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE.
128 * @return time in minutes from zero Z
129 */
130 private static double morningPhenomenon(double julian, double latitude,
131 double longitude, double zenithDistance) {
132 double t = julianDayToJulianCenturies(julian);
133 double eqtime = equationOfTime(t);
134 double solarDec = sunDeclination(t);
135 double hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
136 double delta = longitude - Math.toDegrees(hourangle);
137 double timeDiff = 4 * delta;
138 double timeUTC = 720 + timeDiff - eqtime;
139
140 // Second pass includes fractional julian day in gamma calc
141 double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t)
142 + timeUTC / 1440);
143 eqtime = equationOfTime(newt);
144 solarDec = sunDeclination(newt);
145 hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
146 delta = longitude - Math.toDegrees(hourangle);
147 timeDiff = 4 * delta;
148
149 return 720 + timeDiff - eqtime;
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 return 720 + timeDiff - eqtime;
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 return (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 }
211
212 /**
213 * Convert Julian Day to centuries since J2000.0
214 *
215 * @param julian
216 * The Julian Day to convert
217 * @return the value corresponding to the Julian Day
218 */
219 private static double julianDayToJulianCenturies(double julian) {
220 return (julian - 2451545) / 36525;
221 }
222
223 /**
224 * Convert centuries since J2000.0 to Julian Day
225 *
226 * @param t
227 * Number of Julian centuries since J2000.0
228 * @return The Julian Day corresponding to the value of t
229 */
230 private static double julianCenturiesToJulianDay(double t) {
231 return (t * 36525) + 2451545;
232 }
233
234 /**
235 * Calculate the difference between true solar time and mean solar time
236 *
237 * @param t
238 * Number of Julian centuries since J2000.0
239 * @return
240 */
241 private static double equationOfTime(double t) {
242 double epsilon = obliquityCorrection(t);
243 double l0 = geomMeanLongSun(t);
244 double e = eccentricityOfEarthsOrbit(t);
245 double m = geometricMeanAnomalyOfSun(t);
246 double y = Math.pow((Math.tan(Math.toRadians(epsilon) / 2)), 2);
247
248 double eTime = y * Math.sin(2 * Math.toRadians(l0)) - 2 * e
249 * Math.sin(Math.toRadians(m)) + 4 * e * y
250 * Math.sin(Math.toRadians(m))
251 * Math.cos(2 * Math.toRadians(l0)) - 0.5 * y * y
252 * Math.sin(4 * Math.toRadians(l0)) - 1.25 * e * e
253 * Math.sin(2 * Math.toRadians(m));
254
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
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 }