001    /*
002     * Zmanim Java API
003     * Copyright (C) 2004-2011 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 - 2011
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 this.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                    double adjustedZenith = zenith;
054                    if(adjustForElevation){
055                            adjustedZenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation().getElevation());
056                    } else {
057                            adjustedZenith = adjustZenith(zenith, 0);
058                    }
059    
060                    // step 1: First calculate the day of the year
061                    // NOT NEEDED in this implementation
062    
063                    // step 2: convert the longitude to hour value and calculate an
064                    // approximate time
065                    double lngHour = astronomicalCalendar.getGeoLocation().getLongitude() / 15;
066    
067                    double t = astronomicalCalendar.getCalendar().get(Calendar.DAY_OF_YEAR)
068                                    + ((6 - lngHour) / 24); // use 18 for
069                    // sunset instead
070                    // of 6
071    
072                    // step 3: calculate the sun's mean anomaly
073                    double m = (0.9856 * t) - 3.289;
074    
075                    // step 4: calculate the sun's true longitude
076                    double l = m + (1.916 * Math.sin(Math.toRadians(m)))
077                                    + (0.020 * Math.sin(Math.toRadians(2 * m))) + 282.634;
078                    while (l < 0) {
079                            double Lx = l + 360;
080                            l = Lx;
081                    }
082                    while (l >= 360) {
083                            double Lx = l - 360;
084                            l = Lx;
085                    }
086    
087                    // step 5a: calculate the sun's right ascension
088                    double RA = Math.toDegrees(Math.atan(0.91764 * Math.tan(Math
089                                    .toRadians(l))));
090    
091                    while (RA < 0) {
092                            double RAx = RA + 360;
093                            RA = RAx;
094                    }
095                    while (RA >= 360) {
096                            double RAx = RA - 360;
097                            RA = RAx;
098                    }
099    
100                    // step 5b: right ascension value needs to be in the same quadrant as L
101                    double lQuadrant = Math.floor(l / 90) * 90;
102                    double raQuadrant = Math.floor(RA / 90) * 90;
103                    RA = RA + (lQuadrant - raQuadrant);
104    
105                    // step 5c: right ascension value needs to be converted into hours
106                    RA /= 15;
107    
108                    // step 6: calculate the sun's declination
109                    double sinDec = 0.39782 * Math.sin(Math.toRadians(l));
110                    double cosDec = Math.cos(Math.asin(sinDec));
111    
112                    // step 7a: calculate the sun's local hour angle
113                    double cosH = (Math.cos(Math.toRadians(adjustedZenith)) - (sinDec * Math
114                                    .sin(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude()))))
115                                    / (cosDec * Math.cos(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude())));
116    
117                    // the following line would throw an Exception if the sun never rose.
118                    // this is not needed since the calculation will return a Double.NaN
119                    // if (cosH > 1) throw new Exception("doesnthappen");
120    
121                    // FOR SUNSET use the following instead of the above if statement.
122                    // if (cosH < -1)
123    
124                    // step 7b: finish calculating H and convert into hours
125                    double H = 360 - Math.toDegrees(Math.acos(cosH));
126    
127                    // FOR SUNSET remove "360 - " from the above
128    
129                    H = H / 15;
130    
131                    // step 8: calculate local mean time
132    
133                    double T = H + RA - (0.06571 * t) - 6.622;
134    
135                    // step 9: convert to UTC
136                    double UT = T - lngHour;
137                    while (UT < 0) {
138                            double UTx = UT + 24;
139                            UT = UTx;
140                    }
141                    while (UT >= 24) {
142                            double UTx = UT - 24;
143                            UT = UTx;
144                    }
145                    return UT;
146            }
147    
148            /**
149             * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar,
150             *      double, boolean)
151             */
152            public double getUTCSunset(AstronomicalCalendar astronomicalCalendar,
153                            /*GeoLocation geoLocation,*/ double zenith, boolean adjustForElevation) {
154                    // zenith = adjustZenithForElevation(astronomicalCalendar, zenith,
155                    // geoLocation.getElevation());
156                    // double elevationAdjustment = this.getElevationAdjustment(zenith,
157                    // geoLocation.getElevation());
158                    // double refractionAdjustment = this.getRefraction(zenith);
159                    // zenith = zenith + elevationAdjustment + refractionAdjustment;
160                    double adjustedZenith = zenith;
161                    if(adjustForElevation){
162                            adjustedZenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation().getElevation());
163                    } else {
164                            adjustedZenith = adjustZenith(zenith, 0);
165                    }
166    
167                    // step 1: First calculate the day of the year
168                    // int calendarDayOfYear = calelendar.DAY_OF_YEAR;
169    
170                    // int N=theday - date(1,1,theday.year()) + 1;
171                    int N = astronomicalCalendar.getCalendar().get(Calendar.DAY_OF_YEAR);
172    
173                    // step 2: convert the longitude to hour value and calculate an
174                    // approximate time
175                    double lngHour = astronomicalCalendar.getGeoLocation().getLongitude() / 15;
176    
177                    double t = N + ((18 - lngHour) / 24);
178    
179                    // step 3: calculate the sun's mean anomaly
180                    double M = (0.9856 * t) - 3.289;
181    
182                    // step 4: calculate the sun's true longitude
183                    double L = M + (1.916 * Math.sin(Math.toRadians(M)))
184                                    + (0.020 * Math.sin(Math.toRadians(2 * M))) + 282.634;
185                    while (L < 0) {
186                            double Lx = L + 360;
187                            L = Lx;
188                    }
189                    while (L >= 360) {
190                            double Lx = L - 360;
191                            L = Lx;
192                    }
193    
194                    // step 5a: calculate the sun's right ascension
195                    double RA = Math.toDegrees(Math.atan(0.91764 * Math.tan(Math
196                                    .toRadians(L))));
197                    while (RA < 0) {
198                            double RAx = RA + 360;
199                            RA = RAx;
200                    }
201                    while (RA >= 360) {
202                            double RAx = RA - 360;
203                            RA = RAx;
204                    }
205    
206                    // step 5b: right ascension value needs to be in the same quadrant as L
207                    double Lquadrant = Math.floor(L / 90) * 90;
208                    double RAquadrant = Math.floor(RA / 90) * 90;
209                    RA = RA + (Lquadrant - RAquadrant);
210    
211                    // step 5c: right ascension value needs to be converted into hours
212                    RA /= 15;
213    
214                    // step 6: calculate the sun's declination
215                    double sinDec = 0.39782 * Math.sin(Math.toRadians(L));
216                    double cosDec = Math.cos(Math.asin(sinDec));
217    
218                    // step 7a: calculate the sun's local hour angle
219                    double cosH = (Math.cos(Math.toRadians(adjustedZenith)) - (sinDec * Math
220                                    .sin(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude()))))
221                                    / (cosDec * Math.cos(Math.toRadians(astronomicalCalendar.getGeoLocation().getLatitude())));
222    
223                    // the following line would throw an Exception if the sun never set.
224                    // this is not needed since the calculation will return a Double.NaN
225                    // if (cosH < -1) throw new ZmanimException("doesnthappen");
226    
227                    // step 7b: finish calculating H and convert into hours
228                    double H = Math.toDegrees(Math.acos(cosH));
229                    H = H / 15;
230    
231                    // step 8: calculate local mean time
232    
233                    double T = H + RA - (0.06571 * t) - 6.622;
234    
235                    // step 9: convert to UTC
236                    double UT = T - lngHour;
237                    while (UT < 0) {
238                            double UTx = UT + 24;
239                            UT = UTx;
240                    }
241                    while (UT >= 24) {
242                            double UTx = UT - 24;
243                            UT = UTx;
244                    }
245                    return UT;
246            }
247    }