sunriseset<-function(date,lat,long,timez) { radian<-pi/180 phi<-lat*radian lambda<-long/15 n<-as.numeric(date-julian(12,31,as.numeric(as.character(years(date)))-1 ) ) dst<-0 s<-0 t<- n+(18-12*s-lambda)/24 m<- (.9856*t-3.289)*radian l<- m+(1.916*sin(m)+.02*sin(2*m)+282.634)*radian delta<- .39782*sin(l) sunset<- 3600*( ((24*((1-2*s) *acos((-.01454-sin(delta)*sin(phi)) /(cos(delta)*cos(phi))) /(2*pi)+s)+12*(atan(.91746*tan(l))/pi +ifelse( (floor(l*2/pi)+4)%%4==1 | (floor(l*2/pi)+4)%%4==2,1,0)) -0.06571*t-6.622-lambda+timez+dst+48))%%24) s<-1 t<- n+(18-12*s-lambda)/24 m<- (.9856*t-3.289)*radian l<- m+(1.916*sin(m)+.02*sin(2*m)+282.634)*radian delta<- .39782*sin(l) sunrise<- 3600*( ((24*((1-2*s) *acos((-.01454-sin(delta)*sin(phi)) /(cos(delta)*cos(phi))) /(2*pi)+s)+12*(atan(.91746*tan(l))/pi +ifelse( (floor(l*2/pi)+4)%%4==1 | (floor(l*2/pi)+4)%%4==2,1,0)) -0.06571*t-6.622-lambda+timez+dst+48))%%24) return(date,lat,long,timez,sunrise,sunset) }