スネルの法則を使って屈折方向の計算 - 流れに従ひて己を失はず．

スネルの法則で3次元空間上の屈折方向ベクトルを求めるのをPythonで実装．

ベクトルで解く

各要素の関係は以下．ただし各ベクトルの大きさを1とする．
入射ベクトル： $v_{inter}$
法線ベクトル： $v_{norm}$
出射ベクトル(これを求める)： $v_{exit}$
入射角： $\theta_{inter}$
出射角(これもわかってない)： $\theta_{exit}$
入射側屈折率： $n_{inter}$
入射側屈折率： $n_{exit}$

f:id:Yasutchi:20190914234142j:plain

スネルの法則
$n_{exit} sin\theta_{exit} = n_{inter} sin\theta_{inter}. \tag{0}$
つまり， $\theta_{exit}$ はぐりぐりやればわかる．

ここで，入射ベクトル，法線ベクトル，出射ベクトルは同一平面上にあるから， $\alpha, \beta$ (ともにスカラ)を使って以下のように表現できる．

$v_{exit} = \alpha\ v_{inter} + \beta\ v_{norm}. \tag{1}$

まず $\alpha$ から求める．
式(1)の両辺について，法線ベクトルとの外積をとって，更にその大きさをとると，各ベクトルの大きさが1であること，最後は式(0)を用いて $\alpha$ が決まる．
$\begin{eqnarray} |v_{norm} \times v_{exit}| &=& |v_{norm} \times (\alpha\ v_{inter} + \beta\ v_{norm})|,\\|v_{norm} \times v_{exit}| &=& |v_{norm} \times \alpha\ v_{inter} + v_{norm} \times \beta\ v_{norm}|,\\|v_{norm}||v_{exit}|sin(\pi - \theta_{exit}) &=& \alpha|v_{norm}||v_{inter}|sin(\pi - \theta_{inter}) + 0,\\sin\theta_{exit} &=& \alpha sin\theta_{inter},\\ \alpha &=& \dfrac{sin\theta_{exit}}{sin\theta_{inter}} = \dfrac{n_{inter}}{n_{exit}}. \end{eqnarray}$

次に $\beta$ を求める．
式(1)の両辺について，法線ベクトルとの内積をとると，先程と同様に各ベクトルの大きさが1であることを用いて $\beta$ が決まる． $\theta_{exit}$ が残っているが，式(0)から求めることができる．
$\begin{eqnarray} v_{norm} \cdot v_{exit} &=& v_{norm} \cdot (\alpha\ v_{inter} + \beta\ v_{norm}),\\v_{norm} \cdot v_{exit} &=& \alpha\ v_{norm} \cdot v_{inter} + \beta\ v_{norm} \cdot v_{norm},\\ |v_{norm}||v_{exit}|cos(\pi - \theta_{exit}) &=& \alpha|v_{norm}||v_{inter}|cos(\pi - \theta_{inter}) + \beta,\\-cos\theta_{exit} &=& -\alpha cos\theta_{inter} + \beta,\\ \beta &=& \alpha cos\theta_{inter} -cos\theta_{exit}. \end{eqnarray}$

Python実装

import numpy as np 

def snell(v_inter,v_norm,n_inter,n_exit):
    #3次元スネルの公式の適用結果を返す関数
    #入力：入射光線の方向ベクトル，入射点での法線ベクトル，入射元の屈折率，入射先の屈折率
    
    #入力されたベクトルの正規化
    v_inter = v_inter / np.linalg.norm(v_inter)
    v_norm = v_norm / np.linalg.norm(v_norm)
    
    #係数a
    a =  n_inter/n_exit
    
    cos_t_inter = -1.0*np.dot(v_inter,v_norm)
    #量子化誤差対策
    if cos_t_inter<-1.:
        cos_t_inter = -1.
    elif cos_t_inter>1.:
        cos_t_inter = 1.
        
    sin_t_inter = (1.0 - cos_t_inter**2)**(1/2)
    sin_t_exit = sin_t_inter*a
    
    if sin_t_exit>1.0:
        #全反射する場合
        return np.zeros(3), False
    
    cos_t_exit = (1 - sin_t_exit**2)**(1/2)
    
    #係数b
    b = -1.0*cos_t_exit + a*cos_t_inter
    
    #出射光線の方向ベクトル
    v_exit = a*v_inter + b*v_norm
    #正規化
    v_exit = v_exit / np.linalg.norm(v_exit)
    
    return v_exit, True

vi = np.array([0.08584,0.17301,0.9811726])
vn = np.array([0.050, 0.060, -0.9969453])
ni = 1.0
ne = 1.5

snell(vi,vn,ni,ne)#(array([0.0401461 , 0.0948433 , 0.99468238]), True)

#方向ベクトルなどの実値参考：
#https://www.cybernet.co.jp/optical/course/optics/opt02/opt06.html

参考：http://yees.g.dgdg.jp/RayTrace/RayTraceVersion001b.pdf