I tried to implement a calculation to determine the Az and El angles needed to point an antenna at a given Earth Station coordinates with respect to a GEO Satellite. The source of the mathematical equations needed is a document from Intelsat which can be found under:
https://celestrak.org/NORAD/elements/supplemental/IESS_412_Rev_3.pdf
The calculation procedure is described in chapter 2.5
And my code looks like this:
'''
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,"°")
'''
So far so good but i get a strange error ot lets say unexpected results for a specific input.
For example if i try to point the antenna towards a Satellite with orbital position 50° West, then the output looks like:
Circumference of earth(km) 42164.57 km Input Location of Satellite in degrees:-50 Your input: 50.0 ° West Input Longitude of ES in degrees:11.66 Your input: 11.66 ° Input Latitude of ES in degrees:48.26 Your input: 48.26 ° If North Hemi input 180:180 Input for height over seaside level in km:0.68 Rstation in km: 6390.059844850744 km
Slantrange in km: 40596.37254163326 km
Azimuth angle in degrees: 248.10660270274292 °
Elevation angle in degrees: 1471.9705470065662 °
The Elevation angle is completely out of a reasonable range.
Up to an orbital position of the satellite of 49° West it works just fine. Radius of earth (km) 6378.14 km Circumference of earth(km) 42164.57 km Input Location of Satellite in degrees:-49 Your input: 49.0 ° West Input Longitude of ES in degrees:11.66 Your input: 11.66 ° Input Latitude of ES in degrees:48.26 Your input: 48.26 ° If North Hemi input 180:180 Input for height over seaside level in km:0.68 Rstation in km: 6390.059844850744 km
Slantrange in km: 40528.75697667383 km
Azimuth angle in degrees: 247.2745538278463 °
Elevation angle in degrees: 10.5904943454838 °
I just do not get it why is this happening. May be someone has a hint for me.
I believe this line in your code
x = ELg + m.degrees(0.589)
Should be changed to :
x = ELg + 0.589
your Elg is in degrees. 0.589 is also degrees. the function m.degrees(0.589) will try to convert radians to degrees. But in the documentation you provided 0.589 is degrees.