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 GNU Lesser General
006 * Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option)
007 * any later version.
008 *
009 * This library is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; without even the implied
010 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
011 * details.
012 * You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to
013 * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA,
014 * or connect to: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html
015 */
016package net.sourceforge.zmanim.util;
017
018import java.util.Calendar;
019
020/**
021 * Implementation of sunrise and sunset methods to calculate astronomical times. This implementation is a port of the
022 * C++ algorithm written by Ken Bloom for the sourceforge.net <a
023 * href="http://sourceforge.net/projects/zmanim/">Zmanim</a> project. Ken's algorithm is based on the US Naval Almanac
024 * algorithm. Added to Ken's code is adjustment of the zenith to account for elevation. Originally released under the
025 * GPL, it has been released under the LGPL as of April 8, 2010.
026 * 
027 * @author &copy; Chanoch (Ken) Bloom 2003 - 2004
028 * @author &copy; Eliyahu Hershfeld 2004 - 2011
029 * @version 1.1
030 */
031public class ZmanimCalculator extends AstronomicalCalculator {
032        private String calculatorName = "US Naval Almanac Algorithm";
033
034        /**
035         * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getCalculatorName()
036         */
037        public String getCalculatorName() {
038                return this.calculatorName;
039        }
040
041        /**
042         * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean)
043         */
044        public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
045                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
046                double adjustedZenith = adjustZenith(zenith, elevation);
047
048                // step 1: First calculate the day of the year
049                // NOT NEEDED in this implementation
050
051                // step 2: convert the longitude to hour value and calculate an approximate time
052                double lngHour = geoLocation.getLongitude() / 15;
053
054                double t = calendar.get(Calendar.DAY_OF_YEAR) + ((6 - lngHour) / 24); // use 18 for sunset instead of 6
055
056                // step 3: calculate the sun's mean anomaly
057                double m = (0.9856 * t) - 3.289;
058
059                // step 4: calculate the sun's true longitude
060                double l = m + (1.916 * Math.sin(Math.toRadians(m))) + (0.020 * Math.sin(Math.toRadians(2 * m))) + 282.634;
061                while (l < 0) {
062                        double Lx = l + 360;
063                        l = Lx;
064                }
065                while (l >= 360) {
066                        double Lx = l - 360;
067                        l = Lx;
068                }
069
070                // step 5a: calculate the sun's right ascension
071                double RA = Math.toDegrees(Math.atan(0.91764 * Math.tan(Math.toRadians(l))));
072
073                while (RA < 0) {
074                        double RAx = RA + 360;
075                        RA = RAx;
076                }
077                while (RA >= 360) {
078                        double RAx = RA - 360;
079                        RA = RAx;
080                }
081
082                // step 5b: right ascension value needs to be in the same quadrant as L
083                double lQuadrant = Math.floor(l / 90) * 90;
084                double raQuadrant = Math.floor(RA / 90) * 90;
085                RA = RA + (lQuadrant - raQuadrant);
086
087                // step 5c: right ascension value needs to be converted into hours
088                RA /= 15;
089
090                // step 6: calculate the sun's declination
091                double sinDec = 0.39782 * Math.sin(Math.toRadians(l));
092                double cosDec = Math.cos(Math.asin(sinDec));
093
094                // step 7a: calculate the sun's local hour angle
095                double cosH = (Math.cos(Math.toRadians(adjustedZenith)) - (sinDec * Math.sin(Math.toRadians(geoLocation
096                                .getLatitude())))) / (cosDec * Math.cos(Math.toRadians(geoLocation.getLatitude())));
097
098                // step 7b: finish calculating H and convert into hours
099                double H = 360 - Math.toDegrees(Math.acos(cosH));
100
101                // FOR SUNSET remove "360 - " from the above
102
103                H = H / 15;
104
105                // step 8: calculate local mean time
106
107                double T = H + RA - (0.06571 * t) - 6.622;
108
109                // step 9: convert to UTC
110                double UT = T - lngHour;
111                while (UT < 0) {
112                        double UTx = UT + 24;
113                        UT = UTx;
114                }
115                while (UT >= 24) {
116                        double UTx = UT - 24;
117                        UT = UTx;
118                }
119                return UT;
120        }
121
122        /**
123         * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean)
124         */
125        public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
126                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
127                double adjustedZenith = adjustZenith(zenith, elevation);
128
129                // step 1: First calculate the day of the year
130                int N = calendar.get(Calendar.DAY_OF_YEAR);
131
132                // step 2: convert the longitude to hour value and calculate an approximate time
133                double lngHour = geoLocation.getLongitude() / 15;
134
135                double t = N + ((18 - lngHour) / 24);
136
137                // step 3: calculate the sun's mean anomaly
138                double M = (0.9856 * t) - 3.289;
139
140                // step 4: calculate the sun's true longitude
141                double L = M + (1.916 * Math.sin(Math.toRadians(M))) + (0.020 * Math.sin(Math.toRadians(2 * M))) + 282.634;
142                while (L < 0) {
143                        double Lx = L + 360;
144                        L = Lx;
145                }
146                while (L >= 360) {
147                        double Lx = L - 360;
148                        L = Lx;
149                }
150
151                // step 5a: calculate the sun's right ascension
152                double RA = Math.toDegrees(Math.atan(0.91764 * Math.tan(Math.toRadians(L))));
153                while (RA < 0) {
154                        double RAx = RA + 360;
155                        RA = RAx;
156                }
157                while (RA >= 360) {
158                        double RAx = RA - 360;
159                        RA = RAx;
160                }
161
162                // step 5b: right ascension value needs to be in the same quadrant as L
163                double Lquadrant = Math.floor(L / 90) * 90;
164                double RAquadrant = Math.floor(RA / 90) * 90;
165                RA = RA + (Lquadrant - RAquadrant);
166
167                // step 5c: right ascension value needs to be converted into hours
168                RA /= 15;
169
170                // step 6: calculate the sun's declination
171                double sinDec = 0.39782 * Math.sin(Math.toRadians(L));
172                double cosDec = Math.cos(Math.asin(sinDec));
173
174                // step 7a: calculate the sun's local hour angle
175                double cosH = (Math.cos(Math.toRadians(adjustedZenith)) - (sinDec * Math.sin(Math.toRadians(geoLocation
176                                .getLatitude())))) / (cosDec * Math.cos(Math.toRadians(geoLocation.getLatitude())));
177
178                // step 7b: finish calculating H and convert into hours
179                double H = Math.toDegrees(Math.acos(cosH));
180                H = H / 15;
181
182                // step 8: calculate local mean time
183
184                double T = H + RA - (0.06571 * t) - 6.622;
185
186                // step 9: convert to UTC
187                double UT = T - lngHour;
188                while (UT < 0) {
189                        double UTx = UT + 24;
190                        UT = UTx;
191                }
192                while (UT >= 24) {
193                        double UTx = UT - 24;
194                        UT = UTx;
195                }
196                return UT;
197        }
198}