001    /*
002     * Zmanim Java API
003     * Copyright (C) 2004-2010 Eliyahu Hershfeld
004     *
005     * This library is free software; you can redistribute it and/or modify it under the terms of the
006     * GNU Lesser General Public License as published by the Free Software Foundation; either
007     * version 2.1 of the 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 Lesser 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    
021    import net.sourceforge.zmanim.AstronomicalCalendar;
022    
023    /**
024     * Implementation of sunrise and sunset methods to calculate astronomical times.
025     * This implementation is a port of the C++ algorithm written by Ken Bloom for
026     * the sourceforge.net <a href="http://sourceforge.net/projects/zmanim/">Zmanim</a>
027     * project. Ken's algorithm is based on the US Naval Almanac algorithm. Added to
028     * Ken's code is adjustment of the zenith to account for elevation. Originally released
029     * under the GPL, it has been released under the LGPL as of April 8, 2010.
030     *
031     * @author &copy; Chanoch (Ken) Bloom 2003 - 2004
032     * @author &copy; Eliyahu Hershfeld 2004 - 2010
033     * @version 1.1
034     */
035    public class ZmanimCalculator extends AstronomicalCalculator {
036            private String calculatorName = "US Naval Almanac Algorithm";
037            public String getCalculatorName(){
038                    return calculatorName; //"US Naval Almanac Algorithm";
039            }
040    
041            /**
042             * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(AstronomicalCalendar,
043             *      double, boolean)
044             */
045            public double getUTCSunrise(AstronomicalCalendar astronomicalCalendar,
046                            /*GeoLocation geoLocation,*/ double zenith, boolean adjustForElevation) {
047                    // zenith = adjustZenithForElevation(astronomicalCalendar, zenith,
048                    // geoLocation.getElevation());
049                    // double elevationAdjustment = this.getElevationAdjustment(zenith,
050                    // geoLocation.getElevation());
051                    // double refractionAdjustment = this.getRefraction(zenith);
052                    // zenith = zenith + elevationAdjustment + refractionAdjustment;
053                    if(adjustForElevation){
054                            zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation().getElevation());
055                    } else {
056                            zenith = adjustZenith(zenith, 0);
057                    }
058    
059                    // step 1: First calculate the day of the year
060                    // NOT NEEDED in this implementation
061    
062                    // step 2: convert the longitude to hour value and calculate an
063                    // approximate time
064                    double lngHour = astronomicalCalendar.getGeoLocation().getLongitude() / 15;
065    
066                    double t = astronomicalCalendar.getCalendar().get(Calendar.DAY_OF_YEAR)
067                                    + ((6 - lngHour) / 24); // use 18 for
068                    // sunset instead
069                    // of 6
070    
071                    // step 3: calculate the sun's mean anomaly
072                    double m = (0.9856 * t) - 3.289;
073    
074                    // step 4: calculate the sun's true longitude
075                    double l = m + (1.916 * Math.sin(Math.toRadians(m)))
076                                    + (0.020 * Math.sin(Math.toRadians(2 * m))) + 282.634;
077                    while (l < 0) {
078                            double Lx = l + 360;
079                            l = Lx;
080                    }
081                    while (l >= 360) {
082                            double Lx = l - 360;
083                            l = Lx;
084                    }
085    
086                    // step 5a: calculate the sun's right ascension
087                    double RA = Math.toDegrees(Math.atan(0.91764 * Math.tan(Math
088                                    .toRadians(l))));
089    
090                    while (RA < 0) {
091                            double RAx = RA + 360;
092                            RA = RAx;
093                    }
094                    while (RA >= 360) {
095                            double RAx = RA - 360;
096                            RA = RAx;
097                    }
098    
099                    // step 5b: right ascension value needs to be in the same quadrant as L
100                    double lQuadrant = Math.floor(l / 90) * 90;
101                    double raQuadrant = Math.floor(RA / 90) * 90;
102                    RA = RA + (lQuadrant - raQuadrant);
103    
104                    // step 5c: right ascension value needs to be converted into hours
105                    RA /= 15;
106    
107                    // step 6: calculate the sun's declination
108                    double sinDec = 0.39782 * Math.sin(Math.toRadians(l));
109                    double cosDec = Math.cos(Math.asin(sinDec));
110    
111                    // step 7a: calculate the sun's local hour angle
112                    double cosH = (Math.cos(Math.toRadians(zenith)) - (sinDec * Math
113                                    .sin(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude()))))
114                                    / (cosDec * Math.cos(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude())));
115    
116                    // the following line would throw an Exception if the sun never rose.
117                    // this is not needed since the calculation will return a Double.NaN
118                    // if (cosH > 1) throw new Exception("doesnthappen");
119    
120                    // FOR SUNSET use the following instead of the above if statement.
121                    // if (cosH < -1)
122    
123                    // step 7b: finish calculating H and convert into hours
124                    double H = 360 - Math.toDegrees(Math.acos(cosH));
125    
126                    // FOR SUNSET remove "360 - " from the above
127    
128                    H = H / 15;
129    
130                    // step 8: calculate local mean time
131    
132                    double T = H + RA - (0.06571 * t) - 6.622;
133    
134                    // step 9: convert to UTC
135                    double UT = T - lngHour;
136                    while (UT < 0) {
137                            double UTx = UT + 24;
138                            UT = UTx;
139                    }
140                    while (UT >= 24) {
141                            double UTx = UT - 24;
142                            UT = UTx;
143                    }
144                    return UT;
145            }
146    
147            /**
148             * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar,
149             *      double, boolean)
150             */
151            public double getUTCSunset(AstronomicalCalendar astronomicalCalendar,
152                            /*GeoLocation geoLocation,*/ double zenith, boolean adjustForElevation) {
153                    // zenith = adjustZenithForElevation(astronomicalCalendar, zenith,
154                    // geoLocation.getElevation());
155                    // double elevationAdjustment = this.getElevationAdjustment(zenith,
156                    // geoLocation.getElevation());
157                    // double refractionAdjustment = this.getRefraction(zenith);
158                    // zenith = zenith + elevationAdjustment + refractionAdjustment;
159    
160                    if(adjustForElevation){
161                            zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation().getElevation());
162                    } else {
163                            zenith = adjustZenith(zenith, 0);
164                    }
165    
166                    // step 1: First calculate the day of the year
167                    // int calendarDayOfYear = calelendar.DAY_OF_YEAR;
168    
169                    // int N=theday - date(1,1,theday.year()) + 1;
170                    int N = astronomicalCalendar.getCalendar().get(Calendar.DAY_OF_YEAR);
171    
172                    // step 2: convert the longitude to hour value and calculate an
173                    // approximate time
174                    double lngHour = astronomicalCalendar.getGeoLocation().getLongitude() / 15;
175    
176                    double t = N + ((18 - lngHour) / 24);
177    
178                    // step 3: calculate the sun's mean anomaly
179                    double M = (0.9856 * t) - 3.289;
180    
181                    // step 4: calculate the sun's true longitude
182                    double L = M + (1.916 * Math.sin(Math.toRadians(M)))
183                                    + (0.020 * Math.sin(Math.toRadians(2 * M))) + 282.634;
184                    while (L < 0) {
185                            double Lx = L + 360;
186                            L = Lx;
187                    }
188                    while (L >= 360) {
189                            double Lx = L - 360;
190                            L = Lx;
191                    }
192    
193                    // step 5a: calculate the sun's right ascension
194                    double RA = Math.toDegrees(Math.atan(0.91764 * Math.tan(Math
195                                    .toRadians(L))));
196                    while (RA < 0) {
197                            double RAx = RA + 360;
198                            RA = RAx;
199                    }
200                    while (RA >= 360) {
201                            double RAx = RA - 360;
202                            RA = RAx;
203                    }
204    
205                    // step 5b: right ascension value needs to be in the same quadrant as L
206                    double Lquadrant = Math.floor(L / 90) * 90;
207                    double RAquadrant = Math.floor(RA / 90) * 90;
208                    RA = RA + (Lquadrant - RAquadrant);
209    
210                    // step 5c: right ascension value needs to be converted into hours
211                    RA /= 15;
212    
213                    // step 6: calculate the sun's declination
214                    double sinDec = 0.39782 * Math.sin(Math.toRadians(L));
215                    double cosDec = Math.cos(Math.asin(sinDec));
216    
217                    // step 7a: calculate the sun's local hour angle
218                    double cosH = (Math.cos(Math.toRadians(zenith)) - (sinDec * Math
219                                    .sin(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude()))))
220                                    / (cosDec * Math.cos(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude())));
221    
222                    // the following line would throw an Exception if the sun never set.
223                    // this is not needed since the calculation will return a Double.NaN
224                    // if (cosH < -1) throw new ZmanimException("doesnthappen");
225    
226                    // step 7b: finish calculating H and convert into hours
227                    double H = Math.toDegrees(Math.acos(cosH));
228                    H = H / 15;
229    
230                    // step 8: calculate local mean time
231    
232                    double T = H + RA - (0.06571 * t) - 6.622;
233    
234                    // step 9: convert to UTC
235                    double UT = T - lngHour;
236                    while (UT < 0) {
237                            double UTx = UT + 24;
238                            UT = UTx;
239                    }
240                    while (UT >= 24) {
241                            double UTx = UT - 24;
242                            UT = UTx;
243                    }
244                    return UT;
245            }
246    }