\(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
ロッシュポテンシャルとラグランジュ点
連星で、2つの成分星を\(\,M_1,\,M_2\,\)、星の重心間の距離を\(\,a\,\)とし、\(M_1\,\)星の重心を原点、2つの星の重心を結ぶ線を\(\,x\,\)軸、軌道面内に\(\,y\,\)軸とし、軌道は円軌道を仮定する。ここで任意の点\(\,(x,y,z)\,\)に働く力は2つの星の万有引力と遠心力の全ポテンシャル\(\,\Psi\,\)で表される。
\(\quad\Psi=-G\dfrac{M_1}{r_1}-G\dfrac{M_2}{r_2}-\dfrac{\Omega^2}{2}\left\{\left(x-\dfrac{M_2a}{M_1+M_2}\right)^2+y^2\right\}\qquad\cdots\,\)(1)
ここで
\(\quad r_1^2=x^2+y^2+z^2\quad,\,\,r_2^2=(a-x)^2+y^2+z^2\qquad\cdots\,\)(2)
ここで、\(\psi(x,y,z)\,\)を次のように定義する。
\(\quad\dfrac{GM_1}{a}\psi(x,y,z)=\Psi\qquad\cdots\,\)(3)
質量比\(\,q=M_2/M_1\,\)を使い、軌道面に限定して考えると
\(\quad\psi=-\dfrac{1}{r_1}-\dfrac{q}{r_2}-\dfrac{(1+q)(x^2+y^2)}{2}-\dfrac{q^2}{2(1+q)}\qquad\cdots\,\)(4)
これを 3D サーフェス表示と等高線で表示する。
また、ラグランジュ点は式(4)を微分して、\(\partial\psi/\partial x=0\,\)とした以下の式から求められる
\(\quad -\dfrac{x}{r_1^3}+\left(\dfrac{1-x}{r_2^3}-1\right)q+(1+q)x=0\qquad\cdots\,\)(5)
\(L_1\,\)を求めるには、\(r_1=x\,,\quad r_2=1-x\,\,\)を代入して通分すると
\(\quad -(1+q)x^5+(2+3q)x^4-(1+3q)x^3+x^2-2x+1=0\qquad\cdots\,\)(6)
\(L_2\,\)を求めるには、\(r_1=x\,,\quad r_2=x-1\,\,\)を代入して通分すると
\(\quad (1+q)x^5-(2+3q)x^4+(1+3q)x^3-(1+2q)x^2+2x-1=0\qquad\cdots\,\)(7)
(6)、(7)式をニュートン法で解く。
Python の実行例は 地球－月系の計算結果である。
◆ 参考文献
- W Ser型食連星はくちょう座V367星の 測光ならびに分光観測：小木美奈子(2012)、岡山理科大学
- LagrangePoint詳解-2: c Hiromu Nakagawa
import numpy as np
from matplotlib import pyplot as plt
import japanize_matplotlib
import warnings
import sys
# 警告非表示
warnings.simplefilter('ignore')
mx = 40
my = 30
pot = np.zeros((my,mx)) # ポテンシャルの配列
aa = np.zeros(6) # 5次方程式の係数
xx = np.linspace(-4.,4.,40)
yy = np.linspace(-2.,4.,30)
X, Y = np.meshgrid(xx,yy)
def newton(x,a): # 5次方程式をニュートン法で解く
x1 = x # 初期値
for i in range(20):
y1 = a[5]*x1**5 + a[4]*x1**4 + a[3]*x1**3 + a[2]*x1**2 + a[1]*x1 + a[0]
y2 = 5.*a[5]*x1**4 + 4.*a[4]*x1**3 + 3.*a[3]*x1**2 + 2.*a[2]*x1 + a[1]
if abs(y2) < 1.e-4:
return 0
dx = y1/y2
if abs(dx) < 1.e-4: # 差分が一定値以下なら
return x1
x1 = x1 -dx
return x1
def poten(x,y,q): # ロッシュポテンシャルを計算する
r1 = x**2 + y**2
r2 = (1.-x)**2 + y**2
zz = 1./np.sqrt(r1) + q/np.sqrt(r2) + .5*(1.+q)*r1 -q*x + .5*q**2/(1+q)
return -zz
q = float(input( '質量比( M2/M1) : '))
if q < 1.e-2 or q > 0.95:
sys.exit() # 入力エラーで終了
# ラグランジュ点の計算
aa[0] = 1
aa[1] = -2
aa[2] = 1
aa[3] = -1 - 3 * q
aa[4] = 2 + 3 * q
aa[5] = -1 - q
x1 = 0.5 # L1 の初期値
x2 = newton(x1,aa)
print('\n ---> L1 の座標：','{:.3f}'.format(x2),' (伴星のX座標を1)')
aa[0] = -1
aa[1] = 2
aa[2] = -1 - 2 * q
aa[3] = 1 + 3 * q
aa[4] = -2 - 3 * q
aa[5] = 1 + q
x1 = 1.2 # L2 の初期値
x2 = newton(x1,aa)
print('\n ---> L2 の座標：','{:.3f}'.format(x2),' (伴星のX座標を1)')
# ロッシュポテンシャルの計算
pot = poten(X, Y, q)
fig=plt.figure(figsize=(15, 10))
ax = plt.subplot(111, projection='3d')
ax.plot_surface(X, Y, pot, cmap='viridis')
ax.set_zlim(-30,0)
cs = ax.contourf(X,Y,pot,zdir='z',offset=-30, cmap='autumn')
plt.colorbar(cs)
plt.show()
質量比( M2/M1) : 0.0123 ---> L1 の座標： 0.849 (伴星のX座標を1) ---> L2 の座標： 1.168 (伴星のX座標を1)
