Sunset and Sunrise 日出與日落

介紹

在中央氣象局，我們可以查到台灣地區每天日出、日落的時刻。 在timeanddate網站則可查到世界各地的資料。這些資料是職業天文工作者，利用觀察資料配合模型計算出的。
我們用地心赤道坐標系統、太陽軌跡、地球自轉等資料、觀察地A的經緯度及時區等資料，估算A地點的整年日出日落時間。其中，太陽的軌跡與地球自轉的資料，取自運算平台 R 的套件astrolibR。用附錄中的程式，我們算出台北與布宜諾斯艾利斯2016年每月1日的日出日落時分，與網站上的資料相比，相差全在1分鐘之內。

模型與坐標

地球公轉的平面叫做黃道面，地球赤道所在的平面叫做赤道面，黃道面與赤道面的交角\(\epsilon\)就是所謂的地軸傾角。兩面的交線與交角都會隨時間改變。 參考下面動畫，目前地軸傾角約為23.44度。

我們採用所謂的 地心赤道坐標系統 ，將地球中心固定放在原點，赤道平面就是XY平面，地心往北極方向為Z軸，地球以Z軸正向為軸自轉。右手握拳拇指伸出，拇指方向為Z軸時，四指的方向是自轉的方向。請參考下面的動畫

太陽軌道所在的平面叫做黃道面，黃道面與赤道面的交角23.44度，就是地軸傾角\(\epsilon\)。一年之中太陽有兩次運行到赤道面上，分別是春分與秋分。春分之前，太陽在赤道面的下方，之後在赤道面的上方。在春分或秋分時刻，地球南北半球各洽有一半是白天，一半是夜晚。春分之後太陽照到北半球的部分開始大於南半球，並且差距越來越大，直到夏至。此時，太陽在離赤道面最高的地方。夏至後，差距開始縮小，直到秋分，太陽照到地球南北半球的部分再度相等。秋分前後，太陽由赤道面的上方跑到赤道面的下方，秋分之後太陽照到北半球的部分開始小於南半球，並且差距越來越大，直到冬至。此時，太陽在赤道面下方離赤道面最遠的地方。
我們將春分那一瞬間，地球指向太陽的方向為X軸，以地球到太陽的距離為單位，將當時太陽的坐標定為(1,0,0)， 地球的半徑用R表示。有了Z軸與X軸，Y軸方向的訂定如一般右手系的直角坐標系統。地球的半徑用R表示。
在地心赤道坐標系統中，假設物體的直角座標為\((x,y,z)\)，該物體在XY平面上投影的幅角叫做赤經\(\alpha\),\((x,y,z)\)與XY平面的夾角叫做赤緯\(\delta\)， \((x,y,z)\)的長度用\(\rho\)表示，則有下列關係： \[x=\rho \cos\delta \cos\alpha, y=\rho \cos\delta \sin\alpha, z=\rho \sin\delta\]

Julian 日數(JD)

太陽的軌跡以及地球自轉的狀態，需要與時間建立關係。在天文學上，用Julian日數來代表時間。所謂Julian日數(JD)，是以西元前4713年1月1日中午的時刻起算的日數。西元前4713年1月1日12:00的JD=0，西元前4712年1月1日00:00的JD=-0.5。astrolibR套件中的函數jdcnv(year,month,day,hour)可用來算JD。例如，jdcnv(1,2,3,4.1)=1721458.671 代表西元1年2月3日4:06的JD為1721458.671。注意，jdcnv不能用來算西元前的JD。

地面上定點的赤經

A地的經緯度為\(E_A\) 經 \(N_A\) 緯，地球自轉一圈時，A地也繞圓圈轉一圈，圓圈的半徑是\(R \cos N_A\)。 這裡，有一點要注意，地球並不是一天自轉一圈。地球轉動，同一經線所在半平面，連續兩次掃過太陽的時間間隔才是一天。由於太陽一天繞行地球1/365.2425圈，所以地球每天要轉動比一圈多一點點的366.2425/365.2425圈，才能讓同一經線半平面再次掃過太陽。換句話說，地球一年之中自轉了366.2425圈。令\(\theta_A\)代表A地的赤經，則A地的直角坐標\(\vec{A}\)為\((R \cos N_A\cos \theta_A,R\cos N_A\sin\theta_A,R\sin N_A)\)。
\(\theta_A\)隨時間改變，我們可以用astrolibR套件中的函數ct2lst(long,tz,jd)計算之，函數值以時角為單位。這裡的時角是角度單位而非時間單位，1時角度相當於15度。1時角，細分為60分，每分又細分為60秒。例如，ct2lst(0,0,jdcnv(2016,3,20,4.5))=16.3761，代表格林威治時間2016年3月20日4:30， 倫敦(東經0度，時區GMT+0)的赤經為16.3761時角。ct2lst(121.537,8,jdcnv(2016,3,20,4.5)) =0.4786，代表格林威治時間2016年3月20日4:30， 台北(東經121.537度，時區GMT+8)的赤經為0.4786時角。

太陽的仰角與昇落方程式

A地的地平面向上單位法向量\(\vec{n}\)為\(\frac{1}{R}\vec{A}\)。用\(\vec{S}=(s_x,s_y,s_z)\)代表太陽的坐標， \(\vec{S}-\vec{A}\)就是從A地觀察太陽的方向。內積\(\vec{n}\cdot(\vec{S}-\vec{A})\)可用來判定太陽與地平面的相關位置。內積大於0時，代表太陽在地平面上方，是白天；小於0代表太陽在地平面下方，是夜晚； 等於0代表洽在地平面，也就是日出日落時分。因此，A地日出日落時刻滿足方程式 \(\vec{n}\cdot(\vec{S}-\vec{A})=0\)。也就是\(\vec{n}\cdot\vec{S}=R\)， R的值4.25875e-05很小，用0來近似，變成\(\vec{n}\cdot\vec{S}=0\)，其中\(\vec{n}\cdot\vec{S}\)是太陽中心位置與天頂方向的夾角餘弦，夾角餘弦等於0代表太陽中心落在地平線。因為太陽上緣露出地平就算日出，以至於太陽中心仍在地平下方就已看到太陽。又因為光線進入大氣層會發生折射，即使太陽上縁略低於地平就可看見陽光。根據 維基百科，太陽圓盤視角為32分。一般認為，折射影響的角度有34分，再加上太陽的圓盤又貢獻了16分，太陽中心在地平下方角度50分時為日出日落時刻，也就是\(\vec{n}\cdot\vec{S}=\cos(90^{\circ}50')\)。所以，日出日落時刻是解下列方程式(簡稱昇落方程式) \[s_x\cos N_A\cos \theta_A+s_y\cos N_A\sin\theta_A+s_z\sin N_A +\sin(50')=0\] 要用上述方程式求解日出日落時刻，我們已知\(\theta_A\)隨時間的變化，還需要知道太陽隨時間變化的位置。

太陽的軌跡

太陽的位置可用astrolibR套件中的函數sunpos(jd)求得。引數jd是JD，函數值的成分$ra,$dec分別是太陽位置的赤經、赤緯。例如，令s等於sunpos(jdcnv(2016,5,19,12))，則s$ra是格林威治時間2016年5月19日中午12點的時候太陽位置的赤經而s$dec是赤緯，單位是度。

台北的日出日落

將太陽的赤經緯以及台北的緯度及赤經資料，代入昇落方程式求解，解出來的JD就是台北的日出日落時間。數值計算上，我們只要計算一天之中，每一分鐘台北昇落方程式左邊的數值，前後兩分鐘數值由負變正，就找到日出時刻。前後兩分鐘數值由正變負，就找到日落時刻。不過這個時間是格林威治時間，換算成台北時區的時間還要再加上8小時。
附錄中的函數sunriseset(E,N,timezone,month,day)，其函數值就是日出日落的時刻。所需引數，E是觀察地點的經度，東經為正，西經為負，N是緯度，北緯為正，南緯為負，兩者都是度度量。timezone是觀察地的時區，以小時為單位，台北時區為GMT+8，所以台北的timezone為8，布宜諾斯艾利斯的時區為GMT-4，但因實施夏令時間，布宜諾斯艾利斯的timezone為-3。由 latlong 查出城市的經緯度，由 timeanddate 查出時區，就可以用函數sunriseset算出2016年month月份day日的日出日落時刻。
用sunriseset(E,N,timezone,month,day)函數算出台北2016年每月1日的日出日落時間如下表，日出/日落欄位是sunriseset函數計算所得， sunrise/sunset欄位取自 中央氣象局 。表中，可以看出sunriseset的計算與中央氣象局的資料最多相差1分。

月	日	日出	sunrise	日落	sunset
1	1	06:39	06:39	17:15	17:15
2	1	06:37	06:37	17:38	17:38
3	1	06:17	06:17	17:56	17:56
4	1	05:45	05:46	18:10	18:10
5	1	05:18	05:18	18:24	18:24
6	1	05:04	05:04	18:40	18:40
7	1	05:08	05:08	18:48	18:48
8	1	05:21	05:21	18:39	18:39
9	1	05:35	05:35	18:13	18:13
10	1	05:46	05:46	17:41	17:41
11	1	06:01	06:01	17:13	17:13
12	1	06:22	06:22	17:04	17:04

布宜諾斯艾利斯的日出日落

接著用sunriseset(E,N,timezone,month,day)函數算出布宜諾斯艾利斯2016年每月1日的日出日落時間如下表，日出/日落欄位是sunriseset函數計算所得， sunrise/sunset欄位取自 timeanddate 。表中，可以看出sunriseset的計算與timeanddate網站的資料，完全相同。

月	日	日出	sunrise	日落	sunset
1	1	05:44	05:44	20:10	20:10
2	1	06:13	06:13	20:00	20:00
3	1	06:41	06:41	19:30	19:30
4	1	07:06	07:06	18:47	18:47
5	1	07:30	07:30	18:11	18:11
6	1	07:52	07:52	17:51	17:51
7	1	08:01	08:01	17:54	17:54
8	1	07:47	07:47	18:13	18:13
9	1	07:12	07:12	18:35	18:35
10	1	06:30	06:30	18:57	18:57
11	1	05:51	05:51	19:23	19:23
12	1	05:34	05:34	19:52	19:52

附錄

options(digits=14)
library(astrolibR)
Rad<-pi/180

N.taipei<-25.017    
E.taipei<-121.537   
timezone.taipei<- 8
correctionangle<- -50/60

sunheight<-function(t,E,N,timezone){
  jd<-t-timezone/24
  sun<-sunpos(jd)
  sx<-cos(sun$dec*Rad)*cos(sun$ra*Rad)
  sy<-cos(sun$dec*Rad)*sin(sun$ra*Rad)
  sz<-sin(sun$dec*Rad)

  theta<-ct2lst(E,timezone,jd)*pi/12
  sx*cos(N*Rad)*cos(theta)+
  sy*cos(N*Rad)*sin(theta)+
  sz*sin(N*Rad)-sin(correctionangle*Rad)
}

timestring<-function(i){
  min<-i%%60
  hour<-(i-min)/60
  sprintf("%02d:%02d",hour,min)
}

sunriseset<-function(E,N,timezone,month,day){
  minutes<-seq(1,1440)/1440
  day<-jdcnv(2016,month,day,0)
  v1<-sunheight(day+minutes[1],E,N,timezone)
  sunset<-"-"
  sunrise<-"-"
  for (i in 1:1439){
    v2<-sunheight(day+minutes[i+1],E,N,timezone)
    product<-v1*v2
    if (product<0){
      if (abs(v1)<abs(v2))
        j<-i
      else
        j<-i+1
      if (v1<0)
        sunrise<-timestring(j)
      else
        sunset<-timestring(j)
    } else if (product==0) {
      if (v1==0) {
        if (v2>0)
          sunrise<-timestring(i)
        else
          sunset<-timestring(i)
      } else {
        if (v1>0)
          sunset<-timestring(i+1)
        else
          sunrise<-timestring(i+1)
      }
    }
    v1<-v2
  }
  paste0(sunrise,"\t",sunset)
}

msg<-""
for (month in 1:12){
  day<-1
    msg<-paste0(msg,month,'\t',day,'\t',sunriseset(E.taipei,N.taipei,timezone.taipei,month,day),'\n')
}
cat(msg)

E.buenosaires<--58.381559
N.buenosaires<--34.603684
timezone.buenosaires<--3

msg<-""
for (month in 1:12){
  day<-1
    msg<-paste0(msg,month,'\t',day,'\t',sunriseset(E.buenosaires,N.buenosaires,timezone.buenosaires,month,day),'\n')
}
cat(msg)