# Computing the time of sunrise, sunset and transit (skycalc)

Note this is the translation of http://blog.sciencenet.cn/blog-255662-756136.html

This document describes how to compute the time of sunrise, sunset and transit for any place on earth using the skycalc package.

If the longitude and latitude for a place are known, to compute the time of sunrise, sunset and transit, these are the factors considered by the skycalc package:

1. The Right ascension and the Declination of the sun on the day before and after the day. Then the time of sunrise, sunset, transit is interpolated for a specific place (Please refer to the source code of the CAAElliptical_Calculate_Sun function).

2. In astronomy, conventionally, the longitude in the eastern hemisphere is negative (<0), and in the western hemisphere positive (>0).

3. Time zone and daylight saving time also need to be considered during computation so the time could be converted into local time.

## R code

library(skycalc)
## Loading required package: Rcpp
### Convert the julian day into year, month, day, YYYY-MM-DD

Julian2Date <- function(jd){
int2 <- function(v){ return(floor(v))}

r =list()
D = int2(jd+0.5)  F= jd+0.5-D
if(D >=2299161){
c= int2((D-1867216.25)/36524.25)
D= D +1+c- int2(c/4)
}
D  = D +1524
r$Y = int2((D -122.1)/365.25)# Year D = D - int2(365.25*r$Y)
r$M = int2(D/30.601) #Month D = D - int2(30.601*r$M)
r$D = D #Day if(r$M>13){
r$M = r$M -13
r$Y = r$Y -4715
}else{
r$M = r$M - 1
r$Y = r$Y -4716
}  F = F *24
r$h=int2(F) # Hour F = F - r$h  F = F *60
r$m=int2(F) # Minute F = F -r$m  F = F *60
r$s=F # Second return(r) } # A wrapper function for CAARiseTransitSet_Calculate() sunRTScalc <-function(year, month, day, longitude, latitude, zone, bGregorianCalendar =TRUE, type =c("rise", "transit", "set")){ type <- match.arg(type) JD = CAADate_DateToJD(year, month, day, bGregorianCalendar) SunDetails1 = CAAElliptical_Calculate_Sun(JD -1) Alpha1 = SunDetails1$ApparentGeocentricRA
Delta1 = SunDetails1$ApparentGeocentricDeclination; SunDetails2 = CAAElliptical_Calculate_Sun(JD) Alpha2 = SunDetails2$ApparentGeocentricRA
Delta2 = SunDetails2$ApparentGeocentricDeclination; SunDetails3 = CAAElliptical_Calculate_Sun(JD +1) Alpha3 = SunDetails3$ApparentGeocentricRA
Delta3 = SunDetails3$ApparentGeocentricDeclination RiseTransitSetTime =CAARiseTransitSet_Calculate(JD, Alpha1, Delta1, Alpha2, Delta2, Alpha3,Delta3, longitude, latitude, -0.8333) if(type =="rise"){ rtsJD =(JD +(RiseTransitSetTime$details.Rise/24.00))
}   if(type =="set"){
rtsJD =(JD +(RiseTransitSetTime$details.Set/24.00)) } if(type =="transit"){ rtsJD =(JD +(RiseTransitSetTime$details.Transit/24.00))
}

lclJD = rtsJD -(zone/24.00)
return(Julian2Date(lclJD))
}

## Example 1

Compute the time for sunrise, sunset for Boston (Longitude:74.73057W, Latitude: 39.275787N), USA on 2009-08-09

sunRTScalc(year =2009, month =8, day =9,
longitude =74.73057,  latitude =39.275787,
zone =4, bGregorianCalendar=TRUE, type ="rise")
## $Y ## [1] 2009 ## ##$M
## [1] 8
##
## $D ## [1] 9 ## ##$h
## [1] 6
##
## $m ## [1] 6 ## ##$s
## [1] 27.4320968985558

Explanation: As sunset at Boston happens at (local time: 2009-8-9 20h01m39m (UT:2009-8-10: 00h01m39m), of course, you do not know the precise time before calling sunRTScalc, but as soon as you get the result, you need to realise this), therefore you need to call the function as sunRTScalc(year =2009, month =8, day =10,…), because the inputs (year, month, day) are actually converted from Julian Day.

sunRTScalc(year =2009, month =8, day =10,
longitude =74.73057,  latitude =39.275787,
zone = 4, bGregorianCalendar = TRUE, type ="set")
## $Y ## [1] 2009 ## ##$M
## [1] 8
##
## $D ## [1] 9 ## ##$h
## [1] 20
##
## $m ## [1] 1 ## ##$s
## [1] 40.2297654747963

## Example 2

Calculate the time of sunrise and sunset at the Tian'anmen Square in Beijing (Longitude: 116.391325E, Latitude: 39.907356N) on 2014-01-05.

Explanation: As the sun rises at 7h36m15.2s on 2014-01-05 (local time), the UT is actually 2014-01-04 23h36m15.2s. We, therefore, will use 2014-01-04 as inputs. Note the Time zone of Beijing is 8 hours ahead of UT.

sunRTScalc(year =2014, month =1, day =4,
longitude =-116.391325,  latitude =39.907356,
zone =-8, bGregorianCalendar =TRUE, type ="rise")
## $Y ## [1] 2014 ## ##$M
## [1] 1
##
## $D ## [1] 5 ## ##$h
## [1] 7
##
## $m ## [1] 36 ## ##$s
## [1] 15.2419799566269

### Sunset

sunRTScalc(year =2014, month =1, day =5,
longitude =-116.391325, latitude =39.907356,
zone =-8, bGregorianCalendar =TRUE, type ="set")
## $Y ## [1] 2014 ## ##$M
## [1] 1
##
## $D ## [1] 5 ## ##$h
## [1] 17
##
## $m ## [1] 3 ## ##$s
## [1] 16.8268677592278

### Transit

sunRTScalc(year =2014, month =1, day =5,
longitude =-116.391325,  latitude =39.907356,
zone =-8, bGregorianCalendar =TRUE, type ="transit")
## $Y ## [1] 2014 ## ##$M
## [1] 1
##
## $D ## [1] 5 ## ##$h
## [1] 12
##
## $m ## [1] 19 ## ##$s
## [1] 40.6033730506897

## Example 3

Calculate the time of sunset, sunrise for Tian’anmen Square, Beijing (+8hrs) (Longitude: 116.391325E, Latitude: 39.907356N) on 2014-06-20

### Sunrise

Explanation: As the time of sunrise at the site is (Local time 2014-06-20, 4h45m45.14s, UT 2014-06-19 22h45m45.14s), we should use 2014-06-19.

sunRTScalc(year =2014, month =6, day =19,
longitude =-116.391325, latitude =39.907356,
zone =-8, bGregorianCalendar =TRUE, type ="rise")
## $Y ## [1] 2014 ## ##$M
## [1] 6
##
## $D ## [1] 20 ## ##$h
## [1] 4
##
## $m ## [1] 45 ## ##$s
## [1] 45.1418057084084

### Sunset

sunRTScalc(year =2014, month =6, day =20,
longitude =-116.391325, latitude =39.907356,
zone =-8, bGregorianCalendar =TRUE, type ="set")
## $Y ## [1] 2014 ## ##$M
## [1] 6
##
## $D ## [1] 20 ## ##$h
## [1] 19
##
## $m ## [1] 46 ## ##$s
## [1] 5.37074446678162

### Transit

sunRTScalc(year =2014, month =6, day =20,
longitude =-116.391325, latitude =39.907356,
zone =-8, bGregorianCalendar =TRUE, type ="transit")
## $Y ## [1] 2014 ## ##$M
## [1] 6
##
## $D ## [1] 20 ## ##$h
## [1] 12
##
## $m ## [1] 15 ## ##$s
## [1] 54.5454159379005

## Example 4

Compute the time of sunrise, sunset for Melbourn, Australia (Longitude: -144.963171E, Latitude: -37.814247S) on 2020-1-15.

Explanation: The Sunlight saving date starts from the first Sunday of October the year before (herein, 2019), and ends at the first Sunday in April (2020). On January 15th, Australia is in Sunlight saving time, therefore during computation, the time zone is 11 hours ahead of UT. When the sun rises on 2020-14-15 in Melbourn (local time), UT is still on 2020-1-14, we, therefore, need to use sunRTScalc(year =2020, month =1, day =14,...)

### Sunrise

sunRTScalc(year =2020, month =1, day =14,
longitude =-144.963171, latitude =-37.814247,
zone =-11, bGregorianCalendar =TRUE, type ="rise")
## $Y ## [1] 2020 ## ##$M
## [1] 1
##
## $D ## [1] 15 ## ##$h
## [1] 6
##
## $m ## [1] 13 ## ##$s
## [1] 55.4257643222809

### Sunset

sunRTScalc(year =2020, month =1, day =15,
longitude =-144.963171, latitude =-37.814247,
zone =-11, bGregorianCalendar =TRUE, type ="set")
## $Y ## [1] 2020 ## ##$M
## [1] 1
##
## $D ## [1] 15 ## ##$h
## [1] 20
##
## $m ## [1] 44 ## ##$s
## [1] 6.54912650585175

### Transit

sunRTScalc(year =2020, month =1, day =15,
longitude =-144.963171, latitude =-37.814247,
zone =-11, bGregorianCalendar =TRUE, type ="transit")
## $Y ## [1] 2020 ## ##$M
## [1] 1
##
## $D ## [1] 15 ## ##$h
## [1] 13
##
## $m ## [1] 29 ## ##$s
## [1] 13.399246931076

## Session Information

sessionInfo()
## R version 3.5.3 (2019-03-11)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
## [1] skycalc_1.62 Rcpp_1.0.1
##
## loaded via a namespace (and not attached):
##  [1] compiler_3.5.3  magrittr_1.5    tools_3.5.3     htmltools_0.3.6
##  [5] yaml_2.2.0      stringi_1.4.3   rmarkdown_1.12  knitr_1.22
##  [9] stringr_1.4.0   xfun_0.6        digest_0.6.18   evaluate_0.13

Please note that skycalc is just a wrapper package for the AA+ library. The following examples are mainly derived from the documents of the AA+ CPP library written by PJ Naughter. For more information, please refer to http://www.naughter.com/aa.html

https://wap.sciencenet.cn/blog-255662-1176606.html

## 相关博文

GMT+8, 2022-8-18 07:40