直角坐标求轨道六要素

导入关机点参数 ,，以及地球引力常数

import numpy as np
rk=np.array([2222.016, -6132.926, 241.962])
vk=np.array([6.099, 0.948, 3.498])
u=3.986012*10**5

第一步：计算矢径长度和速率大小

,

def getmod_3(x):#获取模长
    return np.sqrt(np.sum(x**2))
def get_sit_2d(a,b):#获取两向量夹角
    return np.arccos(a.dot(b)/(getmod_3(a)*getmod_3(b)))
    
r1=getmod_3(rk)
v1=getmod_3(vk)

第二步：利用活力公式计算a

def geta_v(v, r, u):  # 活力公式求半长轴
    return 1 / (2 / r - v ** 2 / u)
a=geta_v(v1, r1, u)

第三步：计算e和t时刻的E

def getE(r,v,u,a):
    r1=getmod_3(r)
    E1=np.arctan((np.sqrt(1/(u*a))*(r.dot(v)))/(1-r1/a))
    #判断象限
    s=get_sit_2d(r,v)
    xx=13
    if s>0:
        xx=24
    if xx==13:
        if r1>a:
            xx=3
        else:
            xx=1
    elif xx==24:
        if r1>a:
            xx=2
        else:
            xx=4
    if xx==2:E1+=np.pi
    elif xx==3:E1+=np.pi*1.5
    elif xx==4:E1+=np.pi*2
    print(f'在第{xx}象限')
    return E1
def gete(E,r,a):
    return (1-r/a)/np.cos(E)
E=getE(rk,vk,u,a)
e=gete(E,r1,a)

在第2象限

第四步：计算平近地点角M

def getM(E, e):  # 平近地点角
    M1 = E - e * np.sin(E)
    return M1
M=getM(E,e)

第五步：计算轨道倾角i与升交点赤经

1、计算单位质量动量矩

2、计算i、

def geth(r,v):
    return np.cross(r,v)
h=geth(rk,vk)
h1=getmod_3(h)
def getio(hx, hy, hz):  # 轨道倾角与升交点赤经
    i = np.arccos(hz / np.sqrt(hx ** 2 + hy ** 2 + hz ** 2))
    o = np.arctan(-hx / hy)
    if hy<0:
        if o<0:o+=np.pi*2
    else:
        o+=np.pi
    return i, o
i,o=getio(h[0],h[1],h[2])
i=get_sit_2d(h,np.array([0,0,1]))

第六步、计算和

def getsit(e, E):  # 真近地点角
    return 2 * np.arctan(np.sqrt((1 + e) / (1 - e)) * np.tan(E / 2))
def get_omig(o, r, s):
    rabs = getmod_3(r)
    sinu=r[2]/(np.sin(i)*rabs)
    cosu=r[1]*np.sin(o)/rabs + r[0]*np.cos(o)/rabs
    u1=np.arcsin(sinu)
    u2=np.arccos(cosu)
    if abs(u1-u2)<0.001:
        print(u1/np.pi*180)
        if u1<s:
            u1=u1+np.pi-s 
            print(u1/np.pi*180)
        return u1
sit=getsit(e,E)
omega=get_omig(o, rk, sit)

4.2846636983210615
40.627804117683006

最后输出结果

print("轨道六根数如下：")
#print("矢径长度rk=%.5fkm"%r1)
#print("速率大小vk=%.5fkm/s2"%v1)
print("半长轴长度a=%.5fkm"%a)
print("偏心率e=%.5f"%e)
print("偏近点角E=%.5f度"%(E/np.pi*180))
print("平近点角M=%.5f度"%(M/np.pi*180))
print("轨道倾角i=%.5f度"%(i/np.pi*180))
print("升交点赤经Ω=%.5f度"%(o/np.pi*180))
print("真近点角θ=%.5f度"%(sit/np.pi*180))
print("近心点幅角ω=%.5f度"%(omega/np.pi*180))

轨道六根数如下：
半长轴长度a=5551.80385km
偏心率e=0.25335
偏近点角E=133.92300度
平近点角M=123.46741度
轨道倾角i=29.74530度
升交点赤经Ω=286.19410度
真近点角θ=143.65686度
近心点幅角ω=40.62780度