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 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     * &amp; 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 &copy; Eliyahu Hershfeld 2004 - 2008
036     * @author &copy; 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 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 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    }