R - civil Sunrise and sunset time

civilSunriseSunset.R

# Copyright (c) 2016, Iñigo Gonzalez Ponce <igponce (at) gmail>
# All rights reserved.

# Redistribution and use in source and binary forms, with or without 
# modification, are permitted provided that the following conditions are met:

# 1. Redistributions of source code must retain the above copyright notice,
#    this list of conditions and the following disclaimer.

# 2. Redistributions in binary form must reproduce the above copyright notice,
#    this list of conditions and the following disclaimer in the documentation
#    and/or other materials provided with the distribution.

# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

# civilSunriseSunset.R
# Method from: http://williams.best.vwh.net/sunrise_sunset_algorithm.htm

civilSunriseSunset <- function (year, month, day, latitude, longitude, utcOffset = 0) {

	# Sunrise/Sunset Algorithm

	# 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
  
  zenith = (96 / 90) * pi
		
	# 	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.


	# 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

	# 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)
	# t.Rise <- N + ((6 - lngHour) / 24)
	#if setting time is desired:
	#  t = N + ((6 - lngHour) / 24)
	
	t <- c( N + ((6 - lngHour) / 24),  #rising time
	        N + ((18 - lngHour) / 24)   #settint time
	        )

	# 3. calculate the Sun's mean anomaly
	
	M <-(0.9856 * t) - 3.289

	# 4. calculate the Sun's true longitude
	
	L <- M + (1.916 * sin((pi/180)*M)) + (0.020 * sin((pi/180) * 2 * M)) + 282.634
	# NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
	L <- ifelse(L > 360, L-360, L)
	L <- ifelse(L < 0, L+360, L)

	# 5a. calculate the Sun's right ascension
	
	RA <- atan(0.91764 * tan((pi/180)* L))
	# NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
	RA <- ifelse (RA <0, RA+360, RA )
	RA <- ifelse (RA>360, RA-360, RA)

	# 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)

	# 5c. right ascension value needs to be converted into hours
	RA <- RA / 15

	# 6. calculate the Sun's declination
	sinDec <- 0.39782 * sin((pi/180) * L)
	cosDec <- cos((pi/180) * asin(sinDec))

	# 7a. calculate the Sun's local hour angle
	cosH <- (cos((pi/180) * zenith) - (sinDec * sin((pi/180) * latitude))) / (cosDec * cos((pi/180) * latitude))
	
	# if (cosH < -1)
	#  # the sun never sets on this location (on the specified date)
    # if (cosH > 1 )
	#  # the sun never rises on this location (on the specified date)
	
	# 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
	
	H <- c( 360 - (180/pi)* acos((pi/180)*cosH[1]), (180/pi)*acos((pi/180)*cosH[2]) ) / 15

	# 8. calculate local mean time of rising/setting
	
	T <- H + RA - (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
	UT <- ifelse(UT < 0, UT+24, UT)
	UT <- ifelse(UT > 24, UT-24, UT)

	# 10. convert UT value to local time zone of latitude/longitude
	
	localT <- UT + utcOffset

	return( localT )

}