Вычисление элементов орбиты ИСЗ по координатам и составляющим скорости

Вычисляем постоянные площадей
с1=y·z-y∙z
с1=882016.682·-3601.67--1695.30·3929595.924=3485130195.010м2с
с2=z·x-z∙x
с2=3929595.924·-649.28--3601.67·-24172141.973=-89611469980.829м2с
с3=x·y-x∙y
с3=-24172141.973·-1695.30--649.28·882016.682=41551822749.975м2с
c=с12+с22+с32
c=3485130195.0102+-89611469980.8292+41551822749.9752=98837825039.033м2с
2. Вычисляем фокальный параметр
p=c2μ
p=98837825039.03323.9860044∙1014=24508040.328 м
3. Вычисляем угол наклона плоскости орбиты к плоскости экватора
i=arccosс3c
i=arccos41551822749.975 98837825039.033=65º08΄23,6''
4. Вычисляем долготу восходящего узла Ω
sinΩ=c1c∙sini
sinΩ=1153485130195.010 98837825039.033·sin65°08΄23,6''=0.038862184
cosΩ=c2c∙sini
cosΩ=-89611469980.829 98837825039.033·sin65°08΄23,6''=0.99924458
Ω=2º13΄37,9''
5 . Вычисляем квадрат скорости спутника и геоцентрическое расстояние до спутника
v2=x2+y2+z2
v2=-649.282+-1695.302+-3601.672=16267644.527м2с2
r=x2+y2+z2
r=-24172141.9732+390882016.6822+3929595.92492=24505348.908 м
6.Вычисляем постоянную энергии
h=v2-2∙μr
h=16267644.527 -2∙3.9860044∙101424505348.908=-16264064.489м2с2
7. Вычисляем большую полуось орбиты спутника
a=-μh
a=-3.9860044∙1014-16264064.489=24508045.961 м
8. Вычисляем составляющие вектора Лапласа
D=x∙x+y∙y+z∙z
D=-24172141.973·-649.28+882016.682·-1695.30+3929595.924·-3601.67=
D=46116651.714 м
D=μr+h
D =3.9860044∙1014124505348.908+-16264064.489=1790.019
f1=D∙x-D∙x
f1=1790.019·-24172141.973-46116651.714·-649.28=
f1=-13325936371.786
f2=D∙y-D∙y
f2=1790.019·882016.682-46116651.714·-1695.30=
f2=79760603681.649
f3=D∙z-D∙z
f3=1790.019·3929595.924-46116651.714·-3601.67=
f3=173130974517.870
f=f12+f22+f32
f=-13325936371.7862+79760603681.6492+173130974517.8702=191085501326.902 м
9