001/* 002 * Zmanim Java API 003 * Copyright (C) 2004-2023 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: https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html 015 */ 016package com.kosherjava.zmanim.util; 017 018import java.util.Calendar; 019 020/** 021 * Implementation of sunrise and sunset methods to calculate astronomical times. This calculator uses the Java algorithm 022 * written by <a href="htts://web.archive.org/web/20090531215353/http://www.kevinboone.com/suntimes.html">Kevin 023 * Boone</a> that is based on the <a href = "https://aa.usno.navy.mil/">US Naval Observatory's</a><a 024 * href="https://aa.usno.navy.mil/publications/asa">Astronomical Almanac</a> and used with his permission. Added to Kevin's 025 * code is adjustment of the zenith to account for elevation. This algorithm returns the same time every year and does not 026 * account for leap years. It is not as accurate as the Jean Meeus based {@link NOAACalculator} that is the default calculator 027 * use by the KosherJava <em>zmanim</em> library. 028 * 029 * @author © Eliyahu Hershfeld 2004 - 2023 030 * @author © Kevin Boone 2000 031 */ 032public class SunTimesCalculator extends AstronomicalCalculator { 033 034 /** 035 * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getCalculatorName() 036 */ 037 public String getCalculatorName() { 038 return "US Naval Almanac Algorithm"; 039 } 040 041 /** 042 * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean) 043 */ 044 public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 045 double doubleTime = Double.NaN; 046 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 047 double adjustedZenith = adjustZenith(zenith, elevation); 048 doubleTime = getTimeUTC(calendar, geoLocation, adjustedZenith, true); 049 return doubleTime; 050 } 051 052 /** 053 * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean) 054 */ 055 public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 056 double doubleTime = Double.NaN; 057 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 058 double adjustedZenith = adjustZenith(zenith, elevation); 059 doubleTime = getTimeUTC(calendar, geoLocation, adjustedZenith, false); 060 return doubleTime; 061 } 062 063 /** 064 * The number of degrees of longitude that corresponds to one hour time difference. 065 */ 066 private static final double DEG_PER_HOUR = 360.0 / 24.0; 067 068 /** 069 * @param deg the degrees 070 * @return sin of the angle in degrees 071 */ 072 private static double sinDeg(double deg) { 073 return Math.sin(deg * 2.0 * Math.PI / 360.0); 074 } 075 076 /** 077 * @param x angle 078 * @return acos of the angle in degrees 079 */ 080 private static double acosDeg(double x) { 081 return Math.acos(x) * 360.0 / (2 * Math.PI); 082 } 083 084 /** 085 * @param x angle 086 * @return asin of the angle in degrees 087 */ 088 private static double asinDeg(double x) { 089 return Math.asin(x) * 360.0 / (2 * Math.PI); 090 } 091 092 /** 093 * @param deg degrees 094 * @return tan of the angle in degrees 095 */ 096 private static double tanDeg(double deg) { 097 return Math.tan(deg * 2.0 * Math.PI / 360.0); 098 } 099 100 /** 101 * Calculate cosine of the angle in degrees 102 * 103 * @param deg degrees 104 * @return cosine of the angle in degrees 105 */ 106 private static double cosDeg(double deg) { 107 return Math.cos(deg * 2.0 * Math.PI / 360.0); 108 } 109 110 /** 111 * Get time difference between location's longitude and the Meridian, in hours. 112 * 113 * @param longitude the longitude 114 * @return time difference between the location's longitude and the Meridian, in hours. West of Meridian has a negative time difference 115 */ 116 private static double getHoursFromMeridian(double longitude) { 117 return longitude / DEG_PER_HOUR; 118 } 119 120 /** 121 * Calculate the approximate time of sunset or sunrise in days since midnight Jan 1st, assuming 6am and 6pm events. We 122 * need this figure to derive the Sun's mean anomaly. 123 * 124 * @param dayOfYear the day of year 125 * @param hoursFromMeridian hours from the meridian 126 * @param isSunrise true for sunrise and false for sunset 127 * 128 * @return the approximate time of sunset or sunrise in days since midnight Jan 1st, assuming 6am and 6pm events. We 129 * need this figure to derive the Sun's mean anomaly. 130 */ 131 private static double getApproxTimeDays(int dayOfYear, double hoursFromMeridian, boolean isSunrise) { 132 if (isSunrise) { 133 return dayOfYear + ((6.0 - hoursFromMeridian) / 24); 134 } else { // sunset 135 return dayOfYear + ((18.0 - hoursFromMeridian) / 24); 136 } 137 } 138 139 /** 140 * Calculate the Sun's mean anomaly in degrees, at sunrise or sunset, given the longitude in degrees 141 * 142 * @param dayOfYear the day of the year 143 * @param longitude longitude 144 * @param isSunrise true for sunrise and false for sunset 145 * @return the Sun's mean anomaly in degrees 146 */ 147 private static double getMeanAnomaly(int dayOfYear, double longitude, boolean isSunrise) { 148 return (0.9856 * getApproxTimeDays(dayOfYear, getHoursFromMeridian(longitude), isSunrise)) - 3.289; 149 } 150 151 /** 152 * @param sunMeanAnomaly the Sun's mean anomaly in degrees 153 * @return the Sun's true longitude in degrees. The result is an angle >= 0 and <= 360. 154 */ 155 private static double getSunTrueLongitude(double sunMeanAnomaly) { 156 double l = sunMeanAnomaly + (1.916 * sinDeg(sunMeanAnomaly)) + (0.020 * sinDeg(2 * sunMeanAnomaly)) + 282.634; 157 158 // get longitude into 0-360 degree range 159 if (l >= 360.0) { 160 l = l - 360.0; 161 } 162 if (l < 0) { 163 l = l + 360.0; 164 } 165 return l; 166 } 167 168 /** 169 * Calculates the Sun's right ascension in hours. 170 * @param sunTrueLongitude the Sun's true longitude in degrees > 0 and < 360. 171 * @return the Sun's right ascension in hours in angles > 0 and < 360. 172 */ 173 private static double getSunRightAscensionHours(double sunTrueLongitude) { 174 double a = 0.91764 * tanDeg(sunTrueLongitude); 175 double ra = 360.0 / (2.0 * Math.PI) * Math.atan(a); 176 177 double lQuadrant = Math.floor(sunTrueLongitude / 90.0) * 90.0; 178 double raQuadrant = Math.floor(ra / 90.0) * 90.0; 179 ra = ra + (lQuadrant - raQuadrant); 180 181 return ra / DEG_PER_HOUR; // convert to hours 182 } 183 184 /** 185 * Calculate the cosine of the Sun's local hour angle 186 * 187 * @param sunTrueLongitude the sun's true longitude 188 * @param latitude the latitude 189 * @param zenith the zenith 190 * @return the cosine of the Sun's local hour angle 191 */ 192 private static double getCosLocalHourAngle(double sunTrueLongitude, double latitude, double zenith) { 193 double sinDec = 0.39782 * sinDeg(sunTrueLongitude); 194 double cosDec = cosDeg(asinDeg(sinDec)); 195 return (cosDeg(zenith) - (sinDec * sinDeg(latitude))) / (cosDec * cosDeg(latitude)); 196 } 197 198 /** 199 * Calculate local mean time of rising or setting. By 'local' is meant the exact time at the location, assuming that 200 * there were no time zone. That is, the time difference between the location and the Meridian depended entirely on 201 * the longitude. We can't do anything with this time directly; we must convert it to UTC and then to a local time. 202 * 203 * @param localHour the local hour 204 * @param sunRightAscensionHours the sun's right ascention in hours 205 * @param approxTimeDays approximate time days 206 * 207 * @return the fractional number of hours since midnight as a double 208 */ 209 private static double getLocalMeanTime(double localHour, double sunRightAscensionHours, double approxTimeDays) { 210 return localHour + sunRightAscensionHours - (0.06571 * approxTimeDays) - 6.622; 211 } 212 213 /** 214 * Get sunrise or sunset time in UTC, according to flag. This time is returned as 215 * a double and is not adjusted for time-zone. 216 * 217 * @param calendar 218 * the Calendar object to extract the day of year for calculation 219 * @param geoLocation 220 * the GeoLocation object that contains the latitude and longitude 221 * @param zenith 222 * Sun's zenith, in degrees 223 * @param isSunrise 224 * True for sunrise and false for sunset. 225 * @return the time as a double. If an error was encountered in the calculation 226 * (expected behavior for some locations such as near the poles, 227 * {@link Double#NaN} will be returned. 228 */ 229 private static double getTimeUTC(Calendar calendar, GeoLocation geoLocation, double zenith, boolean isSunrise) { 230 int dayOfYear = calendar.get(Calendar.DAY_OF_YEAR); 231 double sunMeanAnomaly = getMeanAnomaly(dayOfYear, geoLocation.getLongitude(), isSunrise); 232 double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); 233 double sunRightAscensionHours = getSunRightAscensionHours(sunTrueLong); 234 double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, geoLocation.getLatitude(), zenith); 235 236 double localHourAngle = 0; 237 if (isSunrise) { 238 localHourAngle = 360.0 - acosDeg(cosLocalHourAngle); 239 } else { // sunset 240 localHourAngle = acosDeg(cosLocalHourAngle); 241 } 242 double localHour = localHourAngle / DEG_PER_HOUR; 243 244 double localMeanTime = getLocalMeanTime(localHour, sunRightAscensionHours, 245 getApproxTimeDays(dayOfYear, getHoursFromMeridian(geoLocation.getLongitude()), isSunrise)); 246 double pocessedTime = localMeanTime - getHoursFromMeridian(geoLocation.getLongitude()); 247 while (pocessedTime < 0.0) { 248 pocessedTime += 24.0; 249 } 250 while (pocessedTime >= 24.0) { 251 pocessedTime -= 24.0; 252 } 253 return pocessedTime; 254 } 255 256 /** 257 * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC) 258 * of <a href="https://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> for the given day at the given location 259 * on earth. This implementation returns solar noon as the time halfway between sunrise and sunset. 260 * {@link NOAACalculator}, the default calculator, returns true solar noon. See <a href= 261 * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of Chatzos</a> for details on solar 262 * noon calculations. 263 * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation) 264 * @see NOAACalculator 265 * 266 * @param calendar 267 * The Calendar representing the date to calculate solar noon for 268 * @param geoLocation 269 * The location information used for astronomical calculating sun times. 270 * @return the time in minutes from zero UTC. If an error was encountered in the calculation (expected behavior for 271 * some locations such as near the poles, {@link Double#NaN} will be returned. 272 */ 273 public double getUTCNoon(Calendar calendar, GeoLocation geoLocation) { 274 double sunrise = getUTCSunrise(calendar, geoLocation, 90, false); 275 double sunset = getUTCSunset(calendar, geoLocation, 90, false); 276 double noon = sunrise + ((sunset - sunrise) / 2); 277 if(noon < 0) { 278 noon += 12; 279 } 280 if(noon < sunrise) { 281 noon -= 12; 282 } 283 return noon; 284 } 285}