This is an implementation of the United States Naval Observatory's algorithm for calculating sunrise/sunset based on latitude, longitude, and current timezone. It is a translation of a solution published on CodeProject by infectuz.ar, which is for MS SQL Server. A big thanks to him, we needed this solution in a project that enables managing solar farms.
DELIMITER //
CREATE FUNCTION fn_get_sunrise_sunset(
_localDate DATETIME,
_latitude DECIMAL(18,3),
_longitude DECIMAL(18,3),
_gmtOffset INT,
_mode ENUM('SUNRISE', 'SUNSET')
)
RETURNS DATETIME #DECIMAL(18,3)
LANGUAGE SQL
DETERMINISTIC
CONTAINS SQL
SQL SECURITY DEFINER
BEGIN
DECLARE _zenith DECIMAL(18,3) DEFAULT 90.83;
/*
https://edwilliams.org/sunrise_sunset_algorithm.htm
Source:
Almanac for Computers, 1990
published by Nautical Almanac Office
United States Naval Observatory
Washington, DC 20392
Inputs:
day, month, year: date of sunrise/sunset
latitude, longitude: location for sunrise/sunset
zenith: Sun's zenith for sunrise/sunset
offical = 90 degrees 50'
civil = 96 degrees
nautical = 102 degrees
astronomical = 108 degrees
NOTE: longitude is positive for East and negative for West
NOTE: the algorithm assumes the use of a calculator with the
trig functions in "degree" (rather than "radian") mode. Most
programming languages assume radian arguments, requiring back
and forth convertions. The factor is 180/pi. So, for instance,
the equation RA = atan(0.91764 * tan(L)) would be coded as RA
= (180/pi)*atan(0.91764 * tan((pi/180)*L)) to give a degree
answer with a degree input for L.
*/
DECLARE _n1, _n2, _n3, _n INT;
DECLARE _dayOfYear INT;
DECLARE _longitudeAsHour DECIMAL(18,3);
DECLARE _t DECIMAL(18,3);
DECLARE _meanAnomaly DECIMAL(18,3);
DECLARE _sunLongitude DECIMAL(18,3);
DECLARE _sunRightAscention DECIMAL(18,3);
DECLARE _lQuadrant DECIMAL(18,3);
DECLARE _rQuadrant DECIMAL(18,3);
DECLARE _sinDeclination DECIMAL(18,3);
DECLARE _cosDeclination DECIMAL(18,3);
DECLARE _cosH DECIMAL(18,3);
DECLARE _h DECIMAL(18,3);
DECLARE _meanTime DECIMAL(18,3);
DECLARE _ut DECIMAL(18,3);
/*
1. first calculate the day of the year
N1 = floor(275 * month / 9)
N2 = floor((month + 9) / 12)
N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
N = N1 - (N2 * N3) + day - 30
*/
SET _n1 = FLOOR(275 * MONTH(_localDate) / 9);
SET _n2 = FLOOR((MONTH(_localDate) + 9) / 12);
SET _n3 = 1 + FLOOR((YEAR(_localDate) - 4 * FLOOR(YEAR(_localDate) / 4) + 2) / 3);
SET _dayOfYear = _n1 - (_n2 * _n3) + DAY(_localDate) - 30;
/*
2. convert the longitude to hour value and calculate an approximate time
lngHour = longitude / 15
if rising time is desired:
t = N + ((6 - lngHour) / 24)
if setting time is desired:
t = N + ((18 - lngHour) / 24)
*/
SET _longitudeAsHour = _longitude / 15;
SET _t = _dayOfYear + (CASE WHEN _mode = 'SUNRISE' THEN 6 ELSE 18 END - _longitudeAsHour) / 24;
/*
3. calculate the Sun's mean anomaly
M = (0.9856 * t) - 3.289
*/
SET _meanAnomaly = (0.9856 * _t) - 3.289;
/*
4. calculate the Sun's true longitude
L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
*/
SET _sunLongitude = _meanAnomaly + (1.916 * SIN(RADIANS(_meanAnomaly))) +
(0.020 * SIN(2 * RADIANS(_meanAnomaly))) + 282.634 - 360;
/*
5a. calculate the Sun's right ascension
RA = atan(0.91764 * tan(L))
NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
*/
SET _sunRightAscention = DEGREES(ATAN(0.91764 * TAN(RADIANS(_sunLongitude))));
/*
5b. right ascension value needs to be in the same quadrant as L
Lquadrant = (floor( L/90)) * 90
RAquadrant = (floor(RA/90)) * 90
RA = RA + (Lquadrant - RAquadrant)
*/
SET _lQuadrant = FLOOR( _sunLongitude/90) * 90;
SET _rQuadrant = FLOOR(_sunRightAscention/90) * 90;
SET _sunRightAscention = _sunRightAscention + (_lQuadrant - _rQuadrant);
/*
5c. right ascension value needs to be converted into hours
RA = RA / 15
*/
SET _sunRightAscention = _sunRightAscention / 15;
/*
6. calculate the Sun's declination
sinDec = 0.39782 * sin(L)
cosDec = cos(asin(sinDec))
*/
SET _sinDeclination = 0.39782 * SIN(RADIANS(_sunLongitude));
SET _cosDeclination = COS(ASIN(_sinDeclination));
/*
7a. calculate the Sun's local hour angle
cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))
if (cosH > 1)
the sun never rises on this location (on the specified date)
if (cosH < -1)
the sun never sets on this location (on the specified date)
*/
SET _cosH = (COS(RADIANS(_zenith)) - (_sinDeclination *
SIN(RADIANS(_latitude)))) / (_cosDeclination * COS(RADIANS(_latitude)));
/*
7b. finish calculating H and convert into hours
if if rising time is desired:
H = 360 - acos(cosH)
if setting time is desired:
H = acos(cosH)
H = H / 15
*/
SET _h = CASE WHEN _mode = 'SUNRISE' THEN 360 - DEGREES(ACOS(_cosH)) ELSE DEGREES(ACOS(_cosH)) END;
SET _h = _h / 15;
/*
8. calculate local mean time of rising/setting
T = H + RA - (0.06571 * t) - 6.622
*/
SET _meanTime = _h + _sunRightAscention - (0.06571 * _t) - 6.622;
/*
9. adjust back to UTC
UT = T - lngHour
NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24
*/
SET _ut = _meanTime - _longitudeAsHour;
/*
10. convert UT value to local time zone of latitude/longitude
localT = UT + localOffset
*/
RETURN fn_get_time_from_decimal(_ut + _gmtOffset / 60, _localDate);
END
//
DELIMITER ;
SELECT fn_get_sunrise_sunset('2012-05-11',
-34.58,
-58.3,
fn_get_current_utc_offset_in_timezone('America/Buenos_Aires'),
'SUNRISE'); # Should be: 2012-05-11 18:01:19.200
SELECT fn_get_sunrise_sunset('2012-05-11',
-34.58,
-58.3,
fn_get_current_utc_offset_in_timezone('America/Buenos_Aires'),
'SUNSET'); # Should be: 2012-05-11 07:37:19.200
The function has one dependency:
DELIMITER //
CREATE FUNCTION fn_get_time_from_decimal(
_timeAsDecimal DECIMAL(28,4),
_date DATE
)
RETURNS DATETIME
LANGUAGE SQL
DETERMINISTIC
CONTAINS SQL
SQL SECURITY DEFINER
BEGIN
DECLARE _hour, _minute, _second, _millisecond INT;
DECLARE _retVal DATETIME;
SET _hour = FLOOR(_timeAsDecimal);
# concatenate hours
SET _minute = FLOOR((_timeAsDecimal - _hour)*60);
# subtract hours and convert to mins
# select mins and secs in dec form subtract mins and convert remainder to seconds:
SET _second = (((_timeAsDecimal-_hour) * 60 - _minute)*60);
SET _millisecond = ((((_timeAsDecimal-_hour) * 60 - _minute) * 60) - _second) * 1000;
IF _hour > 0 THEN
SET _retVal = _date + INTERVAL _hour HOUR;
ELSE # will subtract into the day before
SET _retVal = _date + INTERVAL 1 DAY;
SET _retVal = _retVal + INTERVAL _hour HOUR;
END IF;
SET _retVal = _retVal + INTERVAL _minute MINUTE;
SET _retVal = _retVal + INTERVAL _second SECOND;
SET _retVal = _retVal + INTERVAL _millisecond * 1000 MICROSECOND;
# Return the result of the function
RETURN _retVal;
END
//
DELIMITER ;
SELECT fn_get_time_from_decimal(-1, UTC_TIMESTAMP);
Utility functions to get the utc offset in minutes from time-zone name:
DELIMITER //
CREATE FUNCTION fn_get_current_utc_offset_in_timezone (
_timeZoneName VARCHAR(64)
)
RETURNS INT
LANGUAGE SQL
DETERMINISTIC
CONTAINS SQL
SQL SECURITY DEFINER
BEGIN
RETURN TIMESTAMPDIFF(MINUTE, UTC_TIMESTAMP, CONVERT_TZ(UTC_TIMESTAMP(), 'UTC', COALESCE(_timeZoneName, 'UTC')));
END
//
DELIMITER ;
SELECT fn_get_current_utc_offset_in_timezone('America/Buenos_Aires');
And another to get the UTC offset in minutes for the local timezone:
DELIMITER //
CREATE FUNCTION fn_get_local_utc_offset ()
RETURNS INT
LANGUAGE SQL
DETERMINISTIC
CONTAINS SQL
SQL SECURITY DEFINER
BEGIN
RETURN TIMESTAMPDIFF(MINUTE, UTC_TIMESTAMP, NOW());
END
//
DELIMITER ;
SELECT fn_get_local_utc_offset();
