sunrise_sunset.R

library(RAtmosphere)
library(maptools)
library(lubridate)


# The next function comes from the r-bloggers website
# http://www.r-bloggers.com/approximate-sunrise-and-sunset-times/
suncalc.custom <- function(d,Lat=48.1442,Long=-122.7551, earth.radius=6371, dist.sun=149598000){
	rad<-function(x)pi*x/180 # convert degrees to radians
	epsilon=rad(23.45) # Radians between the xy-plane and the ecliptic plane
	L=rad(Lat) # Convert observer's latitude to radians
	
	## Calculate offset of sunrise based on longitude (min)
	## If Long is negative, then the mod represents degrees West of
	## a standard time meridian, so timing of sunrise and sunset should
	## be made later.
	timezone = -4*(abs(Long)%%15)*sign(Long)
	theta = 2*pi/365.25*(d-80)
	z.s = dist.sun*sin(theta)*sin(epsilon)
	r.p = sqrt(dist.sun^2-z.s^2)
	t0 = 1440/(2*pi)*acos((earth.radius-z.s*sin(L))/(r.p*cos(L)))
	
	##a kludge adjustment for the radius of the sun
	that = t0+5 
	
	## Adjust "noon" for the fact that the earth's orbit is not circular:
	n = 720-10*sin(4*pi*(d-80)/365.25)+8*sin(2*pi*d/365.25)
	
	## now sunrise and sunset are:
	sunrise = (n-that+timezone)/60
	sunset = (n+that+timezone)/60
	return(list("sunrise" = sunrise,"sunset" = sunset))
}

dates <- c(force_tz(ymd("2015-12-30"), tz="CET"),
					 force_tz(ymd("2015-12-30"), tz="EST")
					 )# day in year = 151
# 2 locations: Brussels and Cancun, Mexico
latitudes <- c(50.821, 21.16)
longitudes <- c(4.36, -86.85)

# use RAtmosphere
ratm.result <- suncalc(yday(dates), latitudes, longitudes)$sunrise # location = Brussels
print("result using RAtmosphere:")
print(with_tz(dates + ratm.result * 3600, "UTC"))

# use the custom function
scalc2.result <- suncalc.custom(yday(dates), latitudes, longitudes)$sunrise # location = Brussels
print("result using custom function:")
print(with_tz(dates + scalc2.result * 3600, "UTC"))

# use maptools
mpt.result <- sunriset(matrix(c(latitudes, longitudes), nrow=2), dates, direction="sunrise")
print("result using maptools")
print(with_tz(dates + as.numeric(mpt.result) * 60 * 60 * 24, "UTC"))