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     * &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 - 2011
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 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    }