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第一步：计算矢径长度 r k r_k rk和速率大小 v k v_k vk r k = ∑ i = 1 3 r k i 2 r_k=sqrt{sum_{i=1}^{3} {r}_{ki}^2} rk=i=1∑3rki2 , v k = ∑ i = 1 3 v k i 2 v_k=sqrt{sum_{i=1}^{3} {v}_{ki}^2} vk=i=1∑3vki2
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
v
2
=
μ
(
2
r
−
1
a
)
v^2=mu(frac{2}{r}-frac{1}{a})
v2=μ(r2−a1)
a = 1 2 r − v 2 μ a=frac{1}{frac{2}{r}-frac{v^2}{mu}} a=r2−μv21
def geta_v(v, r, u): # 活力公式求半长轴
return 1 / (2 / r - v ** 2 / u)
a=geta_v(v1, r1, u)
第三步：计算e和t时刻的E
e
sin
E
=
r
⃗
V
⃗
μ
a
esin E=frac{vec{r}vec{V}}{sqrt{mu a}}
esinE=μa
r
V
e
cos
E
=
1
−
r
a
ecos E=1-frac{r}{a}
ecosE=1−ar
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 M = E − e sin ( e ) M=E-esin(e) M=E−esin(e) M = n ( t − τ ) M=n(t-tau) M=n(t−τ)
def getM(E, e): # 平近地点角
M1 = E - e * np.sin(E)
return M1
M=getM(E,e)
第五步：计算轨道倾角i与升交点赤经
Ω
Omega
Ω
1、计算单位质量动量矩
h
⃗
vec{h}
h
h
⃗
=
r
⃗
k
×
v
⃗
k
vec{h}=vec{r}_ktimesvec{v}_k
h
=r
k×v
k
2、计算i、
Ω
Omega
Ω
cos
i
=
h
z
h
cos{i}=frac{h_z}{h}
cosi=hhz
i
=
arccos
h
z
h
i=arccos{frac{h_z}{h}}
i=arccoshhz
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]))
第六步、计算
θ
theta
θ和
ω
omega
ω
tan
θ
2
=
1
+
e
1
−
e
tan
E
2
tanfrac{theta}{2}=sqrt{frac{1+e}{1-e}}tan{frac{E}{2}}
tan2θ=1−e1+e
tan2E
sin
u
=
z
sin
i
r
sin{u}=frac{z}{sin{i}r}
sinu=sinirz
cos
u
=
y
sin
Ω
r
+
x
sin
Ω
r
cos{u}=frac{ysin{Omega}}{r}+frac{xsin{Omega}}{r}
cosu=rysinΩ+rxsinΩ
u
=
ω
+
θ
u=omega+theta
u=ω+θ
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
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度
