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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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

Further reading

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