用skycalc计算日出日落时刻

给定经纬度， 要计算某年某月某日该地点的日出日落时刻， 考虑的因素包括：

  1. 太阳在天球上的位置 （赤经,赤纬）： 一般要获得计算当日， 以及前一日， 后一日的太阳赤经赤纬， 通过插值的方法， 获得太阳在日出日落时刻的准确位置

  2. 城市所在的时区 （天文学上的习惯是， 在格林尼治子午圈以东为负， 以西为正） ，以便将计算获得的世界时转换为当地时区的时间。

这里介绍用skycalc程序包如何计算一个地点， 给定日期的日出日落时刻。

R代码如下： 

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
64
65
66
67
68
library(skycalc)
### 儒略日转换为年月日时分秒的函数
Julian2Date <-function(jd){
int2 <-
function(v){
  return(floor(v));
}
   r =list();
  D = int2(jd+0.5);
  F= jd+0.5-D #取得日数的整数部份A及小数部分F
  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);#年数
  D  = D - int2(365.25*r$Y);  
  r$M = int2(D/30.601); #月数
  D  = D - int2(30.601*r$M);
  r$D = D; #日数
  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);
  F = F - r$h;
  F = F *60;
  r$m=int2(F);
  F = F -r$m;
  F = F *60;
  r$s=F;
  return(r);
}

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))
}
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
### 例一
### 计算美国波士顿 (西经 74.73057,  北纬 39.275787) 2009年8月9日的日出日落时刻
sunRTScalc(year =2009, month =8, day =9, longitude =74.73057,  
          latitude =39.275787, zone =4, bGregorianCalendar=TRUE, type ="rise")
### 因为波士顿的日落已经发生在格林尼治子午圈的下一天 (波士顿当地时间(西四区)2009年8月9日20h01m39m，对应世界时2009年8月10日00h01m39m)， 所里这里在日期上需要加1.

sunRTScalc(year =2009, month =8, day =10, longitude =74.73057,  
          latitude =39.275787, zone =4, bGregorianCalendar=TRUE, type ="set")


### 例二
### 计算北京天安门 (东经 -116.391325,  北纬 39.907356) 2014年1月5日的日出日落时刻
### 日出 因为北京的日出时刻发生在北京时间4h45m43.7s，对应的世界时为1月4日23h36m13.8s， 因此， 在日期上需要减去一天
sunRTScalc(year =2014, month =1, day =4, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="rise")
### 日落
sunRTScalc(year =2014, month =1, day =5, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="set")
### 中天
sunRTScalc(year =2014, month =1, day =5, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="transit")

### 例三
### 计算2014年6月20日北京天安门 (东经 -116.391325,  北纬 39.907356) 的日出日落时刻
### 日出 因为北京的日出时刻发生在北京时间4h45m43.7s，对应的世界时为6月19日22h45m43.7s， 因此， 在日期上需要减去一天
sunRTScalc(year =2014, month =6, day =19, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="rise")
### 日落
sunRTScalc(year =2014, month =6, day =20, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="set")
### 中天
sunRTScalc(year =2014, month =6, day =20, longitude =-116.391325,  
          latitude =39.907356, zone =-8, bGregorianCalendar =TRUE, type ="transit")

####例四
### 计算2020年1月15日澳大利亚墨尔本(东经 -144.963171, 南纬-37.814247) 的日出日落时刻
### 墨尔本为位于东10区，澳大利亚夏时制开始于上一年10月的第一个星期天， 下一年4月第一个星期天止。1月15日，正处于夏令时， 因此时区取-11
### 日出, 由于澳大利亚日出时， 世界时仍然处在2020年1月14日， 因此这里day取14
sunRTScalc(year =2020, month =1, day =14, longitude =-144.963171,  
          latitude =-37.814247, zone =-11, bGregorianCalendar =TRUE, type ="rise")
### 日落
sunRTScalc(year =2020, month =1, day =15, longitude =-144.963171,  
          latitude =-37.814247, zone =-11, bGregorianCalendar =TRUE, type ="set")
### 中天
sunRTScalc(year =2020, month =1, day =15, longitude =-144.963171,  
          latitude =-37.814247, zone =-11, bGregorianCalendar =TRUE, type ="transit")