张金龙的博客分享 http://blog.sciencenet.cn/u/zjlcas 物种适应性、分布与进化

博文

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

已有 3672 次阅读 2019-5-1 22:06 |个人分类:天文地理|系统分类:生活其它

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

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://wap.sciencenet.cn/blog-255662-1176606.html

上一篇:QGIS3 创建包含多边形的shape文件
下一篇:What are we doing in a herbarium? A talk at ISF
收藏 IP: 14.0.236.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-12-14 17:53

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部