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

博文

计算日出日落时刻

已有 7668 次阅读 2014-1-5 14:47 |系统分类:科研笔记|关键词:学者

英文版本参见 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")

         

 

         




https://m.sciencenet.cn/blog-255662-756136.html

上一篇:AAplus天文算法C++库在R中的实现:skycalc程序包
下一篇:苍茫的大地-香港麦理浩径第八段

2 程雅畅 zdlh

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

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

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

GMT+8, 2022-12-6 18:32

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部