001 /*
002 * Zmanim Java API
003 * Copyright (C) 2004-2005 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 java.util.Calendar;
020 import net.sourceforge.zmanim.AstronomicalCalendar;
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.kevinboone.com/suntimes.html">Kevin Boone</a> that is based
026 * on the <a href = "http://aa.usno.navy.mil/">US Naval Observatory's</a><a
027 * href="http://aa.usno.navy.mil/publications/docs/NewASA.htm">Almanac</a> for
028 * Computer algorithm ( <a
029 * href="http://www.amazon.com/exec/obidos/tg/detail/-/0160515106/">Amazon</a>,
030 * <a
031 * href="http://search.barnesandnoble.com/booksearch/isbnInquiry.asp?isbn=0160515106">Barnes
032 * & Noble</a>) and is used with his permission. Added to Kevin's code is
033 * adjustment of the zenith to account for elevation.
034 *
035 * @author © Eliyahu Hershfeld 2004 - 2007
036 * @author © Kevin Boone 2000
037 * @version 1.1
038 */
039
040 public class SunTimesCalculator extends AstronomicalCalculator {
041
042 public String getCalculatorName() {
043 return "US Naval Almanac Algorithm";
044 }
045
046 /**
047 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(AstronomicalCalendar,
048 * double, boolean)
049 */
050 public double getUTCSunrise(AstronomicalCalendar astronomicalCalendar,
051 double zenith, boolean adjustForElevation) {
052 double doubleTime = Double.NaN;
053
054 if (adjustForElevation) {
055 zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation()
056 .getElevation());
057 } else {
058 zenith = adjustZenith(zenith, 0);
059 }
060 doubleTime = getTimeUTC(astronomicalCalendar.getCalendar().get(
061 Calendar.YEAR), astronomicalCalendar.getCalendar().get(
062 Calendar.MONTH) + 1, astronomicalCalendar.getCalendar().get(
063 Calendar.DAY_OF_MONTH), astronomicalCalendar.getGeoLocation()
064 .getLongitude(), astronomicalCalendar.getGeoLocation()
065 .getLatitude(), zenith, TYPE_SUNRISE);
066 return doubleTime;
067 }
068
069 /**
070 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar,
071 * double, boolean)
072 */
073 public double getUTCSunset(AstronomicalCalendar astronomicalCalendar,
074 double zenith, boolean adjustForElevation) {
075 double doubleTime = Double.NaN;
076
077 if (adjustForElevation) {
078 zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation()
079 .getElevation());
080 } else {
081 zenith = adjustZenith(zenith, 0);
082 }
083 doubleTime = getTimeUTC(astronomicalCalendar.getCalendar().get(
084 Calendar.YEAR), astronomicalCalendar.getCalendar().get(
085 Calendar.MONTH) + 1, astronomicalCalendar.getCalendar().get(
086 Calendar.DAY_OF_MONTH), astronomicalCalendar.getGeoLocation()
087 .getLongitude(), astronomicalCalendar.getGeoLocation()
088 .getLatitude(), zenith, TYPE_SUNSET);
089 return doubleTime;
090 }
091
092 /**
093 * Default value for Sun's zenith and true rise/set
094 */
095 public static final double ZENITH = 90 + 50.0 / 60.0;
096
097 private static final int TYPE_SUNRISE = 0;
098
099 private static final int TYPE_SUNSET = 1;
100
101 // DEG_PER_HOUR is the number of degrees of longitude
102 // that corresponds to one hour time difference.
103 private static final double DEG_PER_HOUR = 360.0 / 24.0;
104
105 /**
106 * sin of an angle in degrees
107 */
108 private static double sinDeg(double deg) {
109 return Math.sin(deg * 2.0 * Math.PI / 360.0);
110 }
111
112 /**
113 * acos of an angle, result in degrees
114 */
115 private static double acosDeg(double x) {
116 return Math.acos(x) * 360.0 / (2 * Math.PI);
117 }
118
119 /**
120 * asin of an angle, result in degrees
121 */
122 private static double asinDeg(double x) {
123 return Math.asin(x) * 360.0 / (2 * Math.PI);
124 }
125
126 /**
127 * tan of an angle in degrees
128 */
129 private static double tanDeg(double deg) {
130 return Math.tan(deg * 2.0 * Math.PI / 360.0);
131 }
132
133 /**
134 * cos of an angle in degrees
135 */
136 private static double cosDeg(double deg) {
137 return Math.cos(deg * 2.0 * Math.PI / 360.0);
138 }
139
140 /**
141 * Calculate the day of the year, where Jan 1st is day 1. Note that this
142 * method needs to know the year, because leap years have an impact here
143 */
144 private static int getDayOfYear(int year, int month, int day) {
145 int n1 = 275 * month / 9;
146 int n2 = (month + 9) / 12;
147 int n3 = (1 + ((year - 4 * (year / 4) + 2) / 3));
148 int n = n1 - (n2 * n3) + day - 30;
149 return n;
150 }
151
152 /**
153 * Get time difference between location's longitude and the Meridian, in
154 * hours. West of Meridian has a negative time difference
155 */
156 private static double getHoursFromMeridian(double longitude) {
157 return longitude / DEG_PER_HOUR;
158 }
159
160 /**
161 * Gets the approximate time of sunset or sunrise In _days_ since midnight
162 * Jan 1st, assuming 6am and 6pm events. We need this figure to derive the
163 * Sun's mean anomaly
164 */
165 private static double getApproxTimeDays(int dayOfYear,
166 double hoursFromMeridian, int type) {
167 if (type == TYPE_SUNRISE) {
168 return dayOfYear + ((6.0 - hoursFromMeridian) / 24);
169 } else /* if (type == TYPE_SUNSET) */{
170 return dayOfYear + ((18.0 - hoursFromMeridian) / 24);
171 }
172 }
173
174 /**
175 * Calculate the Sun's mean anomaly in degrees, at sunrise or sunset, given
176 * the longitude in degrees
177 */
178 private static double getMeanAnomaly(int dayOfYear, double longitude,
179 int type) {
180 return (0.9856 * getApproxTimeDays(dayOfYear,
181 getHoursFromMeridian(longitude), type)) - 3.289;
182 }
183
184 /**
185 * Calculates the Sun's true longitude in degrees. The result is an angle
186 * gte 0 and lt 360. Requires the Sun's mean anomaly, also in degrees
187 */
188 private static double getSunTrueLongitude(double sunMeanAnomaly) {
189 double l = sunMeanAnomaly + (1.916 * sinDeg(sunMeanAnomaly))
190 + (0.020 * sinDeg(2 * sunMeanAnomaly)) + 282.634;
191
192 // get longitude into 0-360 degree range
193 if (l >= 360.0) {
194 l = l - 360.0;
195 }
196 if (l < 0) {
197 l = l + 360.0;
198 }
199 return l;
200 }
201
202 /**
203 * Calculates the Sun's right ascension in hours, given the Sun's true
204 * longitude in degrees. Input and output are angles gte 0 and lt 360.
205 */
206 private static double getSunRightAscensionHours(double sunTrueLongitude) {
207 double a = 0.91764 * tanDeg(sunTrueLongitude);
208 double ra = 360.0 / (2.0 * Math.PI) * Math.atan(a);
209 // get result into 0-360 degree range
210 // if (ra >= 360.0) ra = ra - 360.0;
211 // if (ra < 0) ra = ra + 360.0;
212
213 double lQuadrant = Math.floor(sunTrueLongitude / 90.0) * 90.0;
214 double raQuadrant = Math.floor(ra / 90.0) * 90.0;
215 ra = ra + (lQuadrant - raQuadrant);
216
217 return ra / DEG_PER_HOUR; // convert to hours
218 }
219
220 /**
221 * Gets the cosine of the Sun's local hour angle
222 */
223 private static double getCosLocalHourAngle(double sunTrueLongitude,
224 double latitude, double zenith) {
225 double sinDec = 0.39782 * sinDeg(sunTrueLongitude);
226 double cosDec = cosDeg(asinDeg(sinDec));
227
228 double cosH = (cosDeg(zenith) - (sinDec * sinDeg(latitude)))
229 / (cosDec * cosDeg(latitude));
230
231 // Check bounds
232
233 return cosH;
234 }
235
236 /**
237 * Gets the cosine of the Sun's local hour angle for default zenith
238 */
239 private static double getCosLocalHourAngle(double sunTrueLongitude,
240 double latitude) {
241 return getCosLocalHourAngle(sunTrueLongitude, latitude, ZENITH);
242 }
243
244 /**
245 * Calculate local mean time of rising or setting. By `local' is meant the
246 * exact time at the location, assuming that there were no time zone. That
247 * is, the time difference between the location and the Meridian depended
248 * entirely on the longitude. We can't do anything with this time directly;
249 * we must convert it to UTC and then to a local time. The result is
250 * expressed as a fractional number of hours since midnight
251 */
252 private static double getLocalMeanTime(double localHour,
253 double sunRightAscensionHours, double approxTimeDays) {
254 return localHour + sunRightAscensionHours - (0.06571 * approxTimeDays)
255 - 6.622;
256 }
257
258 /**
259 * Get sunrise or sunset time in UTC, according to flag.
260 *
261 * @param year
262 * 4-digit year
263 * @param month
264 * month, 1-12 (not the zero based Java month
265 * @param day
266 * day of month, 1-31
267 * @param longitude
268 * in degrees, longitudes west of Meridian are negative
269 * @param latitude
270 * in degrees, latitudes south of equator are negative
271 * @param zenith
272 * Sun's zenith, in degrees
273 * @param type
274 * type of calculation to carry out {@link #TYPE_SUNRISE} or
275 * {@link #TYPE_SUNRISE}.
276 *
277 * @return the time as a double. If an error was encountered in the
278 * calculation (expected behavior for some locations such as near
279 * the poles, {@link Double.NaN} will be returned.
280 */
281 private static double getTimeUTC(int year, int month, int day,
282 double longitude, double latitude, double zenith, int type) {
283 int dayOfYear = getDayOfYear(year, month, day);
284 double sunMeanAnomaly = getMeanAnomaly(dayOfYear, longitude, type);
285 double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly);
286 double sunRightAscensionHours = getSunRightAscensionHours(sunTrueLong);
287 double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, latitude,
288 zenith);
289
290 double localHourAngle = 0;
291 if (type == TYPE_SUNRISE) {
292 if (cosLocalHourAngle > 1) { // no rise. No need for an Exception
293 // since the calculation
294 // will return Double.NaN
295 }
296 localHourAngle = 360.0 - acosDeg(cosLocalHourAngle);
297 } else /* if (type == TYPE_SUNSET) */{
298 if (cosLocalHourAngle < -1) {// no SET. No need for an Exception
299 // since the calculation
300 // will return Double.NaN
301 }
302 localHourAngle = acosDeg(cosLocalHourAngle);
303 }
304 double localHour = localHourAngle / DEG_PER_HOUR;
305
306 double localMeanTime = getLocalMeanTime(localHour,
307 sunRightAscensionHours, getApproxTimeDays(dayOfYear,
308 getHoursFromMeridian(longitude), type));
309 double pocessedTime = localMeanTime - getHoursFromMeridian(longitude);
310 while (pocessedTime < 0.0) {
311 pocessedTime += 24.0;
312 }
313 while (pocessedTime >= 24.0) {
314 pocessedTime -= 24.0;
315 }
316 return pocessedTime;
317 }
318
319 }