0

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.

Community
  • 1
  • 1
blast
  • 13
  • 4

2 Answers2

1

Like Atanas Atanasov already responded there is the first mistake:

x = ELg + m.degrees(0.589)

In addition i found another two mistakes in the code:

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

it should be:

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

and the definition of the coefficients should be like this:

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

After this adaptations have been done the code is working now as expected.

So the final version looks like:

# -*- 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,"°")
blast
  • 13
  • 4
0

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.

Atanas Atanasov
  • 359
  • 1
  • 10