備忘録としてメモ。
軌道周回内の短周期摂動は地球重力ポテンシャル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