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