Sunset and Sunrise Calculation Based on Latitude, Longitude and UTC Offset for SQL Server 2005

infectuz.ar

Rate me:

5.00/5 (2 votes)

11 May 2012CPOL

24.1K

Determining sunrise and sunset values for a certain latitude and longitude with SQL Server 2005

Introduction

The main goal of this tip is to provide a set of functionality to SQL Server 2005 databases from which sunset and sunrise values will be determined from geopositional information provided as parameters and return a DateTime value.

The main goal for this set of functions is to be used as an astronomic calculator for daylight related calculations that need to be done from DB directly.

Background

This is an implementation of the United States Naval Observatory's sunset/sunrise algorithm for SQL Server and probably (but not tested) compatible with most of the SQL based engines available in the market.

Using the Code

The usage is pretty straight forward. Two main stored procedures are included: DetermineSunrise and DetermineSunset. Each of them will expect date and geo location parameters.

SQL

DECLARE @OffSet int
SET @OffSet = dbo.udf_getUTCOffSet()

-- Buenos Aires, Argentina, 34, 35 S, 58, 22 W, 11th May 2012
--
-- Parameters:
-- -------------
-- Today's date
-- latitude
-- longitude
-- UTC offset

Exec DetermineSunrise
	'05/11/2012', 
	-34.58, 
	-58.3,
	@OffSet
	
Exec DetermineSunset 
	'05/11/2012', 
	-34.58, 
	-58.3,
	@OffSet

OUTPUT

2012-05-11 07:37:19.200
2012-05-11 18:01:19.200

These are the related functions:

SQL

CREATE FUNCTION [dbo].[TimeFromDecimal]
(
	@dTime decimal(28,4),
	@dateTime datetime
)
RETURNS
datetime
AS
BEGIN
-- Add the T-SQL statements to compute the return value here
DECLARE @iHour int, @iMin int, @iSec int, @imSec int
set @iHour = @dTime 
--concatenate hours
set @iMin = ((@dTime - @iHour)*60.) 
--subtract hours and convert to mins
--select mins and secs in dec form subtract mins and convert remainder to seconds:
set @iSec = (((@dTime-@iHour)*60-@iMin)*60)
 
set @imSec = ((((@dTime-@iHour)*60-@iMin)*60) - @iSec)* 1000
 
if @iHour > 0
    begin
        set @dateTime = dateadd(hh, @ihour, @dateTime)
    end
    else --will subtract into the day before
    begin
        set @dateTime = dateadd(d, 1, @datetime)
        set @dateTime = dateadd(hh, @ihour, @dateTime)
    end
 
set @dateTime = dateadd(mi, @imin, @dateTime)
set @dateTime = dateadd(s, @isec, @dateTime)
set @dateTime = dateadd(ms, @imsec, @dateTime)
 
-- Return the result of the function
 
RETURN
@dateTime
 
END

SQL

 CREATE FUNCTION [dbo].[udf_calcDayofYear]
(
    @dDate datetime
)
RETURNS int
AS
BEGIN
 
    DECLARE @N1 int, @N2 int, @N3 int, @N int
 
    set @N1 = floor(275 * (select month(@dDate)) / 9)
    set @N2 = floor(((select month(@dDate)) + 9) / 12)
    set @N3 = (1 + floor(((select year(@dDate)) - 
    4 * floor((select year(@dDate)) / 4) + 2) / 3))
    set @N = @N1 - (@N2 * @N3) + (select day(@dDate)) - 30
 
    RETURN @N
 
END

SQL

CREATE FUNCTION [dbo].[udf_getUTCOffSet]()
RETURNS int
AS
BEGIN
    RETURN DATEDIFF(hh, GETUTCDATE(), GETDATE()) 
END

Stored procedures described in implementation:

SQL

CREATE PROCEDURE [dbo].[DetermineSunrise]
	@localDate datetime,
	@latitude decimal(18,3),
	@longitude decimal(18,3),
	@GMT int
AS
DECLARE @ZENITH decimal(18,3)
SET @ZENITH = 90.83

/*
	  ZENITH
	  offical      = 90 degrees 50'
	  civil        = 96 degrees
	  nautical     = 102 degrees
	  astronomical = 108 degrees
*/
--SUNSET
	--Calculate the day of the year
	DECLARE @DayOfYear int
	SET @DayOfYear = dbo.udf_calcDayofYear(@localDate)

	-- convert the longitude to hour value and calculate an approximate time

	DECLARE @lngHour decimal(18,3)
	DECLARE @t decimal(18,3)

	SET @lngHour = @longitude / 15
	SET @t = @DayOfYear + ((6 - @lngHour) / 24)

	-- calculate the Sun's mean anomaly

	DECLARE @Mean decimal(18,3)	
	SET @Mean = (0.9856 * @t) - 3.289

	-- calculate the Sun's true longitude

	DECLARE @SunLon decimal(18,3)

	SET @SunLon = @Mean + (1.916 * sin(RADIANS(@Mean))) + 
	(0.020 * sin(2 * RADIANS(@Mean))) + 282.634 - 360


	---- calculate the Sun's right ascension

	DECLARE @SunRightAsention decimal(18,3)
	SET @SunRightAsention = DEGREES(atan(0.91764 * tan(RADIANS(@SunLon))))

	---- right ascension value needs to be in the same quadrant as L

	DECLARE @Lquadrant decimal(18,3)
	DECLARE @Rquadrant decimal(18,3)

	SET @Lquadrant  = (floor( @SunLon/90)) * 90
	SET @Rquadrant = (floor(@SunRightAsention/90)) * 90
	SET @SunRightAsention = @SunRightAsention + (@Lquadrant - @Rquadrant)

	----right ascension value needs to be converted into hours

	SET @SunRightAsention = @SunRightAsention / 15

	---- calculate the Sun's declination

	DECLARE @sinDec decimal(18,3)
	DECLARE @cosDec decimal(18,3)

	SET @sinDec = 0.39782 * sin(RADIANS(@SunLon))
	SET @cosDec = cos(asin(@sinDec))

	---- calculate the Sun's local hour angle
	DECLARE @cosH decimal(18,3)	
	SET @cosH = (cos(radians(@ZENITH)) - (@sinDec * 
	sin(radians(@latitude)))) / (@cosDec * cos(radians(@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)

	DECLARE @H decimal(18,3)
	SET @H = 360 - degrees(ACOS(@cosH))
	SET @H = @H / 15

	---- calculate local mean time of rising/setting
	DECLARE @MeanTime decimal(18,3)	
	SET @MeanTime = @H + @SunRightAsention - (0.06571 * @t) - 6.622
		
	DECLARE @UT decimal(18,3)
	SET @UT = @MeanTime - @lngHour

	SELECT dbo.TimeFromDecimal(@UT + @GMT, @localDate)

SQL

CRATE PROCEDURE [dbo].[DetermineSunset]
	@localDate datetime,
	@latitude decimal(18,3),
	@longitude decimal(18,3),
	@GMT int
AS
DECLARE @ZENITH decimal(18,3)
SET @ZENITH = 90.83

/*
	  ZENITH
	  offical      = 90 degrees 50'
	  civil        = 96 degrees
	  nautical     = 102 degrees
	  astronomical = 108 degrees
*/
--SUNSET
	--Calculate the day of the year
	DECLARE @DayOfYear int
	SET @DayOfYear = dbo.udf_calcDayofYear(@localDate)

	-- convert the longitude to hour value and calculate an approximate time

	DECLARE @lngHour decimal(18,3)
	DECLARE @t decimal(18,3)

	SET @lngHour = @longitude / 15
	SET @t = @DayOfYear + ((18 - @lngHour) / 24)

	-- calculate the Sun's mean anomaly

	DECLARE @Mean decimal(18,3)	
	SET @Mean = (0.9856 * @t) - 3.289

	-- calculate the Sun's true longitude

	DECLARE @SunLon decimal(18,3)

	SET @SunLon = @Mean + (1.916 * sin(RADIANS(@Mean))) + 
	(0.020 * sin(2 * RADIANS(@Mean))) + 282.634 - 360


	---- calculate the Sun's right ascension

	DECLARE @SunRightAsention decimal(18,3)
	SET @SunRightAsention = DEGREES(atan(0.91764 * tan(RADIANS(@SunLon))))

	---- right ascension value needs to be in the same quadrant as L

	DECLARE @Lquadrant decimal(18,3)
	DECLARE @Rquadrant decimal(18,3)

	SET @Lquadrant  = (floor( @SunLon/90)) * 90
	SET @Rquadrant = (floor(@SunRightAsention/90)) * 90
	SET @SunRightAsention = @SunRightAsention + (@Lquadrant - @Rquadrant)

	----right ascension value needs to be converted into hours

	SET @SunRightAsention = @SunRightAsention / 15

	---- calculate the Sun's declination

	DECLARE @sinDec decimal(18,3)
	DECLARE @cosDec decimal(18,3)

	SET @sinDec = 0.39782 * sin(RADIANS(@SunLon))
	SET @cosDec = cos(asin(@sinDec))

	---- calculate the Sun's local hour angle
	DECLARE @cosH decimal(18,3)	
	SET @cosH = (cos(radians(@ZENITH)) - (@sinDec * 
	sin(radians(@latitude)))) / (@cosDec * cos(radians(@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)

	DECLARE @H decimal(18,3)
	SET @H = degrees(ACOS(@cosH))
	SET @H = @H / 15

	---- calculate local mean time of rising/setting
	DECLARE @MeanTime decimal(18,3)	
	SET @MeanTime = @H + @SunRightAsention - (0.06571 * @t) - 6.622
		
	DECLARE @UT decimal(18,3)
	SET @UT = @MeanTime - @lngHour


	SELECT dbo.TimeFromDecimal(@UT + @GMT, @localDate)

Hope you find this information useful. Happy coding!

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Written By

infectuz.ar

Argentina

This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

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 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); 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'); 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();

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.