moon phase 月相

月相

在中央氣象局月相圖，我們可以查到每一天的月相資料， 舉凡新月、滿月、上弦月、下弦月等的日子，一目了然。

新月前後三四天，只有10%左右的部分是亮的，亮的部分在右側？左側？右上角？右下角？左上角？左下角？事實上，同一天的不同時刻看，亮的位置是會改變的，前述網站沒有呈現這個現象。

基本上，月亮的形狀主要由地球、月亮、太陽的相關位置以及觀察地點所決定。月亮面對太陽的半個球是亮的，面對地球的半個球是我們(地表觀察者)可以看到的，兩者的交集就是我們看到的部分。至於哪一個角落是亮的，因為地球的自轉，觀察者的天頂方向會隨時間改變，也就是視覺上的上方會隨時間改變，造成視覺上亮的部位改變。

在地平坐標系統下，只要有太陽與月亮的位置資料，我們就可以算出月相。

下面的圖形是用附錄中的程式製作出來的。它是2016年6月3日(陰曆4月28)每一整點在台北可以看到的月相。圖中的月形，藍色代表在地平面上，紅色代表在地平面下，可以明顯看到月相隨間改變的情形。

下面的圖形是2016年6月7日(陰曆5月初三)每一整點在台北可以看到的月相。

地平坐標系統

地平坐標系統以觀察地為原點，觀測目標的位置用距離(\(\rho\))、仰角(\(\phi\))、方位角(\(\theta\))來定位，這是球體坐標。距離指的是觀測站到目標的距離，仰角指的是觀測目標的方向與地面的夾角。方位角指的觀測目標的東西南北方位，以北方為基準，向東轉的角度就是。換句話說，北方的方位角為0度，東方為90度，南方為180度，西方為270度。

將東方、北方、天頂方向分別訂為直角座標系統下的X、Y、Z軸的方向，\((x,y,z)\)為目標的直角座標。直角坐標\((x,y,z)\)與球體坐標\((\rho,\theta,\phi)\)的關係如下： \[x=\rho\cos\phi\sin\theta, y=\rho\cos\phi\cos\theta, z=\sin\phi\]

月球上的臨時座標系統

月球的方位角用\(\theta\)表示、仰角用\(\phi\)表示，則地面望向月亮的方向\(e_2=(\cos\phi\sin\theta,\cos\phi\cos\theta,\sin\phi)\)。從月球位置出發，沿著同一方位角的大圓，往天頂的方向\(e_3\)就是\(e_2\)對\(\phi\)的微分， 也就是\(e_3=(-\sin\phi\sin\theta,-\sin\phi\cos\theta,\cos\phi)\)。令\(e_1=e_2\times e_3\)， 從地表往月亮看，\(e_1\)的方向是月亮的右邊，\(e_3\)的方向是月亮的上邊，\(e_2\)的方向是月亮的後邊。以月亮中心為原點，\(\{e_1,e_2,e_3\}\)為坐標軸方向，我們建立一個臨時坐標系統。

令矩陣\(M\)的第一列為\(e_1\)，第二列為\(e_2\)，第三列為\(e_3\)。在地平坐標系統下，月球到太陽的方向向量坐標用\(S_e\)表示， 則在臨時坐標系統中，月球到太陽的方向向量坐標為\(S_m=M S_e\)。另一方面，在臨時坐標系統中，月球望向地球的方向向量為\((0,-1,0)\)。

在臨時坐標系統中，以月亮的半徑為單位，面向我們的半個月球以大圓\(C_1: x^2+z^2=1, z=0\)為邊界，\(C_1\) 的參數式為\((\cos\theta,0,\sin\theta)\)。

面向太陽的半個月亮的邊界\(C_2\)也是一個大圓，\(C_2\)就是過圓點且與\(S_m\) 垂直的平面\(E\) 跟單位球面相交的部分。令\(B_0=S_m\times (0,1,0)\)，\(B_1\)是與\(B_0\)同向的單位向量，\(B_2=S_m\times B_1\)。 如此，\(B_1,B_2\)是平面\(E\)上兩個互相垂直的單位向量，而 \(C_2\)的參數式為\(\cos\theta\cdot B_1+\sin\theta\cdot B_2\)。

\(C_1,C_2\)兩大圓相交，互相切割，得到4段圓弧，屬於\(C_1\)的兩段用\(C_{11},C_{12}\)表示，屬於\(C_2\)的兩段用\(C_{21},C_{22}\)表示。從\(C_{11},C_{12}\)中挑出靠近太陽的一段(與\(S_m\)的內積大於0)，從\(C_{21},C_{22}\)挑出靠近地球的一段(第二個坐標小於0)，所得兩段圓弧圍得月球表面的部分是月亮中亮的部分。這兩段圓弧在臨時坐標系統中\(XZ\)平面的垂直投影，所圍成的區域就是視覺上的月相。

附錄中，利用python平台的套件ephem與matplotlib，我們製作函數 moonshape_t(lon,lat,tz,year,month,day,hour)來顯示月相。參數lon為經度，東經為正；lat為緯度，北緯為正，tz為時區(以小時為單位，台北為+8)，year為西元，hour採24小時制，可以有小數。 執行moonshape_t('121.537','25.017',8,2016,6,3,16.5)，我們就可得到，2016年6月3日下午4點半在台北(東經121.537度， 北緯25.017度)所看到的月相。

附錄

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt

def spher2cart(theta,phi):
  return np.array([math.cos(phi)*math.sin(theta),math.cos(phi)*math.cos(theta),math.sin(phi)])

def spher2cartDphi(theta,phi):
  return np.array([-math.sin(phi)*math.sin(theta),-math.sin(phi)*math.cos(theta),math.cos(phi)])


def greatCircleWithPar(N):
  # [0,1.0,0] 為觀察者望向月球的方向
  # N 為月球望向太陽的方向
  # 通過月球球心與N垂直的平面切割月球的大圓C，C可以用cos(theta)*b1+sin(theta)*b2表示，
  # 其中的b1,b2是C平面上兩個互相垂直的單位向量
  b1=np.cross(N,np.array([0,1.0,0]))
  if np.linalg.norm(b1)>0.000001:
    b1=1/np.linalg.norm(b1)*b1
    b2=np.cross(N,b1)
  else:
    b1=np.array([1.,0,0])
    b2=np.array([0,0,1.])
  #print([b1,b2])
  def fun(theta):
    return math.cos(theta)*b1+math.sin(theta)*b2
  return fun

def cutarc(points):
  arcs=[]
  arc=[]
  for i in range(len(points)):
    if arc==[]:
      arc.append(points[i])
    else:
      if np.linalg.norm(points[i]-arc[-1])<0.1:
        arc.append(points[i])
      else:
        arcs.append(arc)
        arc=[points[i]]
  arcs.append(arc)  
  return arcs    


def moonshape_t(lon,lat,tz,year,month,day,hour,show=True):
  date=ephem.Date((year,month,day,hour-tz)) # GMT+8
  site=ephem.Observer()
  site.date=date
  site.lon=lon
  site.lat=lat
  sun=ephem.Sun(site)
  moon=ephem.Moon(site)
  # e2 望向月亮的方向
  e2=spher2cart(moon.az,moon.alt)
  # e3 月亮在天球上的往上方向
  e3=spher2cartDphi(moon.az,moon.alt)
  # e1 月亮在天球上的往右方向
  e1=np.cross(e2,e3) 
  M=np.stack([e1,e2,e3])
  moon2sun=sun.earth_distance*spher2cart(sun.az,sun.alt)-moon.earth_distance*spher2cart(moon.az,moon.alt)
  s_e=1/np.linalg.norm(moon2sun)*moon2sun
  # s_m: 以 {e1,e2,e3} 為坐標軸方向，月亮為原點時，太陽的方向
  s_m=M.dot(s_e)
  if moon.alt>0:
    color='b'
  else:
    color='r'
  print(['moon.alt',moon.alt,'moon.az',moon.az])
  f=greatCircleWithPar(s_m)
  thetas=np.linspace(-math.pi,math.pi,num=200)
  plt.plot([math.cos(theta) for theta in thetas],[math.sin(theta) for theta in thetas],'w-')
  points=[f(theta) for theta in thetas if f(theta)[1]<0]
  arcs=cutarc(points)
  for arc in arcs:
    plt.plot([x[0] for x in arc],[x[2] for x in arc],color+'-')
  points=[np.array([math.cos(theta),0.0,math.sin(theta)]) for theta in thetas 
    if np.dot(np.array([math.cos(theta),0.0,math.sin(theta)]),s_m)>0]
  arcs=cutarc(points)
  for arc in arcs:
    plt.plot([x[0] for x in arc],[x[2] for x in arc],color+'-')
  plt.axis('off')
  plt.axis('equal')
  plt.title(str(year)+'/'+str(month)+'/'+str(day)+'/'+str(hour))
  if show:
    plt.show()

lon='121.537'
lat='25.017'
tz=8
moonshape_t(lon,lat,tz,2016,6,3,16)
for i in range(1,25):
  plt.subplot(6,4,i)
  moonshape_t(lon,lat,tz,2016,6,3,i,False)

plt.show()