moon phase 月相

月相

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

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

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

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

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

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

地平坐標系統

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

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

月球上的臨時座標系統

月球的方位角用$\theta$表示、仰角用$\varphi$表示，則地面望向月亮的方向$e_2=(\cos \varphi \sin \theta ,\cos \varphi \cos \theta ,\sin \varphi )$。從月球位置出發，沿著同一方位角的大圓，往天頂的方向$e_3$就是$e_2$對$\varphi$的微分， 也就是$e_3=(-\sin \varphi \sin \theta ,-\sin \varphi \cos \theta ,\cos \varphi )$。令$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()