備忘録としてメモ。
軌道周回内の短周期摂動は地球重力ポテンシャルJ2項が支配的
平均軌道要素 = 接触軌道要素 - 短周期摂動項
数式は、古在大先生の論文から。
def osc2mean(osc_elm):
a,e,i,omega,ra,M = osc_elm[1:]
A = J2*RE**2
E = Ecalc(e, M)
r = a*(1-e*np.cos(E))
nu = 2*np.arctan(np.sqrt((1+e)/(1-e))*np.tan(E/2))
p = a*(1-e**2)
da = (A/a)*((a/r)**3-(1/((1-e**2)**(3/2)))+(-(a/r)**3 + (1/((1-e**2)**(3/2))) + (a/r)**3*np.cos(2*omega+2*nu))*(3/2)*np.sin(i)**2)
de = (A/4)*(-2/(a**2*e*np.sqrt(1-e**2)) + 2*a*(1-e**2)/(e*r**3)
+(3/(a**2*e*np.sqrt(1-e**2))-3*a*(1-e**2)/(e*r**3)
-3*(1-e**2)*np.cos(nu+2*omega)/(p**2)-3*np.cos(2*nu+2*omega)/(a**2*e*(1-e**2))
+3*a*(1-e**2)*np.cos(2*nu+2*omega)/(e*r**3)-(1-e**2)*np.cos(3*nu+2*omega)/p**2)*np.sin(i)**2)
di = (A*np.sin(2*i)/(8*p**2))*(3*np.cos(2*nu+2*omega)+3*e*np.cos(2*omega+nu)+e*np.cos(2*omega+3*nu))
domega = (3*A/(2*p**2))*((2-(5/2)*np.sin(i)**2)*(nu-M+e*np.sin(nu))
+ (1-(3/2)*np.sin(i)**2)*((1/e)*(1-(1/4)*e**2)*np.sin(nu)+(1/2)*np.sin(2*nu)+(e/12)*np.sin(3*nu))
- (1/e)*((1/4)*np.sin(i)**2+(1/2-(15/16)*np.sin(i)**2)*e**2)*np.sin(nu+2*omega)
+ (e/16)*np.sin(i)**2*np.sin(nu-2*omega)-(1/2)*(1-(5/2)*np.sin(i)**2)*np.sin(2*nu+2*omega)
+ (1/e)*((7/12)*np.sin(i)**2-(1/6)*(1-(19/8)*np.sin(i)**2)*e**2)*np.sin(3*nu+2*omega)
+ (3/8)*np.sin(i)**2*np.sin(4*nu+2*omega)+(e/16)*np.sin(i)**2*np.sin(5*nu+2*omega))
draan = (-A*np.cos(i)/(4*p**2))*(6*(nu-M+e*np.sin(nu))-3*np.sin(2*omega+2*nu)-3*e*np.sin(2*omega+nu)-e*np.sin(2*omega+3*nu))
dM = (3*A*np.sqrt(1-e**2)/(2*e*p**2))*(-(1-(3/2)*np.sin(i)**2)*((1-e**2/4)*np.sin(nu)+(e/2)*np.sin(2*nu)+(e**2/12)*np.sin(3*nu))
+ np.sin(i)**2*((1/4)*(1+5*e**2/4)*np.sin(nu+2*omega)-(e**2/16)*np.sin(nu-2*omega)
-(7/12)*(1-e**2/28)*np.sin(3*nu+2*omega)
-(3*e/8)*np.sin(4*nu+2*omega)-(e**2/16)*np.sin(5*nu+2*omega)))
am = a - da
em = e - de
im = i - di
omegam = omega - domega
ram = ra - draan
Mm = M - dM
mean_elm = [0 for i in range(7)]
mean_elm[0] = osc_elm[0]
mean_elm[1:] = am, em, im, omegam, ram, Mm
return mean_elm
なお、平均近点角$M$、動径$r$、離心近点角$E$及び真近点角$\nu$は以下のような関係。離心近点角は解析的に計算できないのでニュートンラフソン法
E - e\sin E = M \\
r = a(1-e\cos E) \\
\nu = 2 \arctan (\tan (E/2) \sqrt (1+e)/(1-e) )
def Ecalc(e, M):
E_init = M - e * math.sin(M)
iteration_num = 10
for i in range(iteration_num):
if i == 0:
E = E_init
E_before = E
else:
E = E - (E - e * math.sin(E) - M) / (1 - e * math.cos(E))
if E - E_before < 0.001:
break
E_before = E
return E
ケプラリアンからカルテシアンへの変換は以下のようにする。
def Kep2Cal(elm):
a, e, i, omega, Omega, M = elm[1:]
mu = 3.986004418 * (10 ** 14) # m^3 s^-2 # 重力定数
T = 2 * math.pi * math.sqrt(a ** 3 / mu) # 軌道周期
n = 2 * math.pi / T
# 離心近点角の計算
E = Ecalc(e, M)
r = np.zeros((3))
f = 2 * math.atan(math.sqrt((1 + e) / (1 - e)) * math.tan(E / 2))
theta = omega + f
r_scalar = a * (1 - e * math.cos(E))
r = r_scalar * np.array([math.cos(Omega) * math.cos(theta) - math.sin(Omega) * math.sin(theta) * math.cos(i)
, math.sin(Omega) * math.cos(theta) + math.cos(Omega) * math.sin(theta) * math.cos(i)
, math.sin(theta) * math.sin(i)])
Xi = r[0]
Yi = r[1]
Zi = r[2]
v = np.zeros((3))
# 角運動量
dEdt = (1 / r_scalar) * math.sqrt(mu / a)
h = a ** 2 * math.sqrt(1 - e ** 2) * (1 - e * math.cos(E)) * dEdt
v = -(mu / h) * np.array([math.cos(Omega) * (math.sin(theta) + e * math.sin(omega)) + math.sin(Omega) * (
math.cos(theta) + e * math.cos(omega)) * math.cos(i)
, math.sin(Omega) * (math.sin(theta) + e * math.sin(omega)) - math.cos(Omega) * (
math.cos(theta) + e * math.cos(omega)) * math.cos(i)
, -(math.cos(theta) + e * math.cos(omega)) * math.sin(i)])
VXi = v[0]
VYi = v[1]
VZi = v[2]
cal = [0 for i in range(7)]
cal[0] = elm[0]
cal[1:] = Xi, Yi, Zi, VXi, VYi, VZi
return cal