博文
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:
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).
In astronomy, conventionally, the longitude in the eastern hemisphere is negative (<0), and in the western hemisphere positive (>0).
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
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
https://m.sciencenet.cn/blog-255662-1176606.html
上一篇：QGIS3 创建包含多边形的shape文件
下一篇：What are we doing in a herbarium? A talk at ISF
