||
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.
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)) }
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
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
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
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
Calculate the time of sunset, sunrise for Tian’anmen Square, Beijing (+8hrs) (Longitude: 116.391325E, Latitude: 39.907356N) on 2014-06-20
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
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
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
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,...)
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
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
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
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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-12-14 17:53
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社