地球固定地心坐标系中地面站坐标的计算

Question

我试图实现一个计算，以确定相对于GEO卫星，将天线指向给定地面站坐标所需的Az和El角。所需的数学方程的来源是国际通信卫星组织的一份文件，可以在下面找到:

https://celestrak.org/NORAD/elements/supplemental/IESS_412_Rev_3.pdf

计算过程见2.5章

我的代码是这样的:

"

import math as m
#import numpy
import mpmath as mp

#Constants
EFc = 1 / 298.257 #flattening constant which recognizes that the polar radius is less than the equatorial radius by
#1 part in 298.257
R = 6378.140 #Radius of earth (km)
print("Radius of earth (km)", R, "km")
aGSO= 42164.57  #Circumference of earth(km)
print("Circumference of earth(km)", aGSO, "km")
#Variable Declaration
Pss = float(input("Input Location of Satellite in degrees:")) #Location of geostationary satellite(degrees)
if Pss > 0:
print("Your input:", Pss,"° East")
else:
print("Your input:", abs(Pss),"° West")

PE= float(input("Input Longitude of ES in degrees:"))   #Longitude of the earth station antenna(degrees)
print("Your input:", PE, "°")
LE= float(input("Input Latitude of ES in degrees:"))    #Latitude of the earth station antenna(degrees)
print("Your input:", LE, "°" )
HS = int(input("If North Hemi input 180:"))


#Calculations
Rs = R / m.sqrt(1-EFc*(2-EFc)*(m.sin(m.radians(LE)))**2)
ho = float(input("Input for height over seaside level in km:"))
Ra = (Rs+ho)*m.cos(m.radians(LE)) #ES radial distance from earth rotation axis
Rz = Rs*((1-EFc)**2)*m.sin(m.radians(LE)) #ES distance above earth equtorial plane
r = aGSO - Rs
rx = aGSO*m.cos(0)*m.cos(m.radians((Pss-PE)))-Ra
ry = aGSO*m.cos(0)*m.sin(m.radians((Pss-PE)))
rz = aGSO*m.sin(0)-Rz
dr_north = -rx*m.sin(m.radians(LE)) + rz*m.cos(m.radians(LE))
dr_zenith = rx*m.cos(m.radians(LE)) + rz*m.sin(m.radians(LE))
SR = m.sqrt(rx**2 + ry**2 + rz**2) #Slantrange in km
Az = m.atan(ry/dr_north) #in Rad
if HS == 180:
Azd = HS + m.degrees(Az)
else:
Azd = m.degrees(Az)
ELg = m.degrees(m.atan(dr_zenith/m.sqrt((dr_north**2+ry**2))))

x = ELg + m.degrees(0.589)
a = 0.58804392
a1 = -0.17941557
a2 = 0.29906946 * 10e-1
a3 = -0.25187400 * 10e-2
a4 = 0.82622101 * 10e-4
if ELg > 10.2:
ELob = ELg + 0.01617 * mp.cot(m.radians(ELg))
elif ELg < 10.2:
ELob = ELg + a + a1*x + a2*x**2 + a3*x**3 + a4*x**4

print("Rstation in km:", Rs, "km")
print()
print("Slantrange in km:",SR,"km")
print()
#print("Ra in km:", Ra, "km")
#print("Rz in km:", Rz, "km")
print("Azimuth angle in degrees:", Azd, "°")
print()
#print(rx)
#print(ry)
#print(rz)
#print(dr_north)
#print(dr_zenith)
print("Elevation angle in degrees:", ELob,"°")

"

到目前为止还好，但是我得到了一个奇怪的错误，让我们说一个特定输入的意外结果。

例如，如果我试图将天线指向轨道位置为50°West的卫星，则输出如下:

地球周长(km) 42164.57 km卫星输入位置:-50度你的输入:西经50.0度输入经度为ES，单位为度:11.66你的输入:11.66°输入ES的纬度，以度为单位:48.26你的输入:48.26°如果北半球输入180:180输入海面以上高度，单位为公里:0.68海拔单位:6390.059844850744 km

以km为单位的Slantrange: 40596.37254163326 km

方位角:248.10660270274292°

仰角:1471.9705470065662°

仰角完全超出合理范围

卫星轨道位置西经49°它工作得很好。地球半径(km) 6378.14 km地球周长(千米)42164.57千米卫星输入位置:-49度你的输入:西经49.0度输入经度为ES，单位为度:11.66你的输入:11.66°输入ES的纬度，以度为单位:48.26你的输入:48.26°如果北半球输入180:180输入海面以上高度，单位为公里:0.68海拔单位:6390.059844850744 km

Slantrange单位:40528.75697667383 km

方位角:247.2745538278463°

仰角(以度为单位):10.5904943454838°

我就是不明白为什么会发生这种事。也许有人能给我一个提示。

Answer 1

就像Atanas一样，Atanasov已经回应了，这是第一个错误:

x = ELg + m.degrees(0.589)

此外，我还发现了代码中的另外两个错误:

Rz = Rs*((1-EFc)**2)*m.sin(m.radians(LE)) #ES distance above earth equtorial plane

应该是:

Rz = (Rs*((1-EFc)**2)+ho)*m.sin(m.radians(LE)) #ES distance above earth equtorial plane

，系数的定义应该是这样的:

a2 = 0.29906946 * 1e-1
a3 = -0.25187400 * 1e-2
a4 = 0.82622101 * 1e-4

完成这些调整后，代码现在按预期工作了。

所以最终版本看起来像:

# -*- coding: utf-8 -*-
"""
Created on Mon May 31 09:57:50 2021
@author: o27503515
"""
import math as m
import numpy as np
import mpmath as mp

#Constants
EFc = 1 / 298.257 #flattening constant which recognizes that the polar radius is less than the equatorial radius by
#1 part in 298.257
R = 6378.140 #Radius of earth (km)
print("Radius of earth (km)", R, "km")
aGSO= 42164.57  #Circumference of earth(km)
print("Circumference of earth(km)", aGSO, "km")
#Variable Declaration
Pss = float(input("Input Location of Satellite in degrees:")) #Location of geostationary satellite(degrees)
if Pss > 0:
print("Your input:", Pss,"° East")
else:
print("Your input:", abs(Pss),"° West")

PE= float(input("Input Longitude of ES in degrees:"))   #Longitude of the earth station antenna(degrees)
print("Your input:", PE, "°")
LE= float(input("Input Latitude of ES in degrees:"))    #Latitude of the earth station antenna(degrees)
print("Your input:", LE, "°" )
HS = int(input("If North Hemi input 180:"))


#Calculations
Rs = R / m.sqrt(1-EFc*(2-EFc)*(m.sin(m.radians(LE)))**2)
ho = float(input("Input for height over seaside level in km:"))
Ra = (Rs+ho)*m.cos(m.radians(LE)) #ES radial distance from earth rotation axis
Rz = (Rs*((1-EFc)**2)+ho)*m.sin(m.radians(LE)) #ES distance above earth equtorial plane
r = aGSO - Rs
rx = aGSO*m.cos(0)*m.cos(m.radians((Pss-PE)))-Ra
ry = aGSO*m.cos(0)*m.sin(m.radians((Pss-PE)))
rz = aGSO*m.sin(0)-Rz
dr_north = -rx*m.sin(m.radians(LE)) + rz*m.cos(m.radians(LE))
dr_zenith = rx*m.cos(m.radians(LE)) + rz*m.sin(m.radians(LE))
SR = m.sqrt(rx**2 + ry**2 + rz**2) #Slantrange in km
Az = m.atan(ry/dr_north) #in Rad
if HS == 180: #In case ES is in North Hemisphere
Azd = HS + m.degrees(Az)
else:
Azd = m.degrees(Az)
ELg = m.degrees(m.atan(dr_zenith/m.sqrt((dr_north**2+ry**2))))

x = ELg + 0.589 
a0 = 0.58804392
a1 = -0.17941557
a2 = 0.29906946 * 1e-1
a3 = -0.25187400 * 1e-2
a4 = 0.82622101 * 1e-4
if ELg > 10.2:
ELob = ELg + 0.01617 * mp.cot(m.radians(ELg))
else:
ELob = ELg + a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4

print("Rstation in km: %.2f" % Rs, "km")
print()
print("Slantrange in km: %.2f" % SR,"km")
print()
#print("Ra in km:", Ra, "km")
#print("Rz in km:", Rz, "km")
print("Azimuth angle in degrees: %.2f" % Azd, "°")
print()
#print(rx)
#print(ry)
#print(rz)
#print(dr_north)
#print(dr_zenith)
print("Elevation angle in degrees: %.2f" % ELob,"°")

Answer 2

我相信这是你代码中的一行x = ELg + m.degrees(0.589)

应该改成:

x = ELg + 0.589

你的生命是以度为单位的。0.589也是度。函数m.degrees(0.589)将尝试将弧度转换为角度。但是在文档中你提供的0.589是度。