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