地球固定地心坐标系中地球站坐标的计算
Calculation of the Earth Station coordinates in the Earth Fixed Geocentric Coordinate System
我尝试进行计算以确定将天线指向相对于 GEO 卫星的给定地球站坐标所需的 Az 和 El 角。所需数学方程式的来源是来自 Intelsat 的文件，可在以下位置找到：
https://celestrak.com/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° 的卫星，则输出如下：
地球周长(km) 42164.57 km
以度为单位的卫星输入位置：-50
您的输入： 50.0 ° West
在degrees:11.66中输入ES的经度
您的输入：11.66°
在 degrees:48.26 中输入 ES 的纬度
您的输入：48.26°
如果北半球输入 180:180
在 km:0.68 中输入高于海边的高度
Rstation 公里：6390.059844850744 公里
斜距公里：40596.37254163326 公里
以度为单位的方位角：248.10660270274292°
仰角度数：1471.9705470065662°
仰角完全超出合理范围
直到卫星轨道位置 49° West 它工作得很好。
地球半径（公里） 6378.14 公里
地球周长(km) 42164.57 km
以度为单位输入卫星位置：-49
您的输入：49.0°西
在degrees:11.66中输入ES的经度
您的输入：11.66°
在 degrees:48.26 中输入 ES 的纬度
您的输入：48.26°
如果北半球输入 180:180
在 km:0.68 中输入高于海边的高度
Rstation 公里：6390.059844850744 公里
斜距公里：40528.75697667383 公里
以度为单位的方位角：247.2745538278463°
仰角度数：10.5904943454838°
我只是不明白为什么会这样。可能有人给我提示。
我相信你代码中的这一行
x = ELg + m.degrees(0.589)
应改为：
x = ELg + 0.589
你的 Elg 是度数。 0.589 也是度数。函数 m.degrees(0.589) 将尝试将弧度转换为度数。但是在您提供的文档中，0.589 是度数。
就像 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,"°")
我尝试进行计算以确定将天线指向相对于 GEO 卫星的给定地球站坐标所需的 Az 和 El 角。所需数学方程式的来源是来自 Intelsat 的文件，可在以下位置找到：
https://celestrak.com/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° 的卫星，则输出如下：
地球周长(km) 42164.57 km 以度为单位的卫星输入位置：-50 您的输入： 50.0 ° West 在degrees:11.66中输入ES的经度 您的输入：11.66° 在 degrees:48.26 中输入 ES 的纬度 您的输入：48.26° 如果北半球输入 180:180 在 km:0.68 中输入高于海边的高度 Rstation 公里：6390.059844850744 公里
斜距公里：40596.37254163326 公里
以度为单位的方位角：248.10660270274292°
仰角度数：1471.9705470065662°
仰角完全超出合理范围
直到卫星轨道位置 49° West 它工作得很好。 地球半径（公里） 6378.14 公里 地球周长(km) 42164.57 km 以度为单位输入卫星位置：-49 您的输入：49.0°西 在degrees:11.66中输入ES的经度 您的输入：11.66° 在 degrees:48.26 中输入 ES 的纬度 您的输入：48.26° 如果北半球输入 180:180 在 km:0.68 中输入高于海边的高度 Rstation 公里：6390.059844850744 公里
斜距公里：40528.75697667383 公里
以度为单位的方位角：247.2745538278463°
仰角度数：10.5904943454838°
我只是不明白为什么会这样。可能有人给我提示。
我相信你代码中的这一行
x = ELg + m.degrees(0.589)
应改为：
x = ELg + 0.589
你的 Elg 是度数。 0.589 也是度数。函数 m.degrees(0.589) 将尝试将弧度转换为度数。但是在您提供的文档中，0.589 是度数。
就像 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,"°")