suncalc<-function(d,Lat=48.1442,Long=-122.7551){ ## d is the day of year ## Lat is latitude in decimal degrees ## Long is longitude in decimal degrees (negative == West) ##This method is copied from: ##Teets, D.A. 2003. Predicting sunrise and sunset times. ## The College Mathematics Journal 34(4):317-321. ## At the default location the estimates of sunrise and sunset are within ## seven minutes of the correct times (http://aa.usno.navy.mil/data/docs/RS_OneYear.php) ## with a mean of 2.4 minutes error. ## Function to convert degrees to radians rad<-function(x)pi*x/180 ##Radius of the earth (km) R=6378 ##Radians between the xy-plane and the ecliptic plane epsilon=rad(23.45) ##Convert observer's latitude to radians L=rad(Lat) ## 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) ## The earth's mean distance from the sun (km) r = 149598000 theta = 2*pi/365.25*(d-80) z.s = r*sin(theta)*sin(epsilon) r.p = sqrt(r^2-z.s^2) t0 = 1440/(2*pi)*acos((R-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)) }
Thursday, October 18, 2007
Approximate sunrise and sunset times
This function is not perfect, but it does a reasonable job estimating sunrise and sunset times for my field site. If more accurate data are required, try here. Note, the command to calculate day of year is: strptime(x, "%m/%d/%Y")$yday+1
Subscribe to:
Post Comments (Atom)
3 comments:
thanks for providing this function, and all these tips!
Anne Ghisla
I'm using this algorithm for a variety of cities and I get NaN errors for cities which are very close to the North Pole:
city=Eureka, Nunavut, country=Canada, latitude=79.9833, longtude=-85.95
city=Alert, Nunavut, country=Canada, latitude=82.4667, longtude=-62.5
city=Nord, country=Greenland, latitude=81.6, longtude=-16.6667
For instance, if I try to calculate for Alert for today (the 52nd day of the year), I have these values in R:
> z.s
[1] -27578819
> r.p
[1] 147033909
> L
[1] 1.395972
> acos(R-z.s*sin(L))/(r.p*cos(L))
[1] NaN
Warning message:
In acos(R - z.s * sin(L)) : NaNs produced
> R-z.s*sin(L)
[1] 27164815
> acos(R-z.s*sin(L))
[1] NaN
Warning message:
In acos(R - z.s * sin(L)) : NaNs produced
>
@John I believe that is because the sun is down all day for those latitudes this time of the year. There is twilight, but no sunrise/set.
Post a Comment