LoginSignup
0
1

More than 3 years have passed since last update.

LEO衛星の短周期摂動と平均軌道要素

Last updated at Posted at 2020-03-16

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
1