英文版本参见 http://blog.sciencenet.cn/blog-255662-1176606.html 

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

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

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

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

R代码如下：

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))

}

### 例一

### 计算美国波士顿 (西经 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")