How can I calculate GPS satellite's Latitude and Longitude if I have the elevation, the azimuth, the XYZ coordinates of the ground station

Question

How can I calculate a GPS satellite's latitude and longitude if the elevation, the azimuth, and the XYZ coordinates of the ground GPS station are provided?

I have tried these questions:

• Computing GEO satellite's longitude from elevation/azimuth from a given latitude/longitude?
• How can I calculate GPS satellite's Latitude and Longitude if I have the elevation and the azimuth

but the answers are not consistent.

I would need a reference if it is available.

Answer 1

I have also been looking for how to perform this calculation for a few days. I think I found a solution that works, courtesy of EASA but copying here in case of link rot.

Inputs

• Observer latitude, $\phi_o$, (instead of XYZ coordinates you mentioned, which I guess was ECEF?)
• Observer longitude, $\lambda_o$
• Elevation angle, $\theta$
• Azimuth angle, $\alpha$
• Satellite altitude, $h$, you didn't mention needing this, but could assume about 20,000 km since GPS satellites are in MEO.

Some assumptions

• Spherical Earth with radius of 6371 km
• Negligible atmospheric effects
• Observer is on the surface

Outputs

• Latitude, $\phi_s$, and longitude, $\lambda_s$, below satellite

Input angles, including latitude and longitude should be in radians depending on software you use. Then perform these calculations:

1. $\Psi = \frac{\pi}{2} - \theta - asin\left(\left(\frac{6371}{6371 + h}\right)cos\,\theta\right)$
2. $\phi_s = asin(sin\,\phi_o*cos\,\Psi + cos\,\phi_o*sin\,\Psi*cos\,\alpha)$

 3a) If either latitude>70$^\circ$ and $tan\,\Psi*cos(\alpha)>tan(\frac{\pi}{2}-\phi_o)$ or latitude<-70$^\circ$ and $tan\,\Psi*cos(\alpha+\pi)>tan(\frac{\pi}{2}+\phi_o)$:

 $\lambda_s = \lambda_o + \pi + asin(\frac{sin\,\Psi*sin\,\alpha}{cos\,\phi_s})$

 3b) Otherwise:

 $\lambda_s = \lambda_o + asin(\frac{sin\,\Psi*sin\,\alpha}{cos\,\phi_s})$

Python Implementation

import numpy as np

def calculate_nadir(observer_lat, observer_lon, satellite_alt, elevation, azimuth):
    # Constants
    Earth_radius = 6371  # Earth's radius in kilometers
    
    # Convert angles from degrees to radians
    observer_lat_rad = np.radians(observer_lat)
    observer_lon_rad = np.radians(observer_lon)
    elevation_rad = np.radians(elevation)
    azimuth_rad = np.radians(azimuth)

    # Calculate intermediate angle
    psi_alt = np.pi/2 - elevation_rad - np.arcsin(Earth_radius/(Earth_radius+satellite_alt)*np.cos(elevation_rad))
    SSP_lat = np.arcsin(np.sin(observer_lat_rad) * np.cos(psi_alt) + 
                        np.cos(observer_lat_rad) * np.sin(psi_alt) * np.cos(azimuth_rad))
    if (observer_lat>70   and  np.tan(psi_alt)*np.cos(azimuth_rad) > np.tan(np.pi/2-observer_lat_rad)) or \
       (observer_lat<-70  and  np.tan(psi_alt)*np.cos(azimuth_rad + np.pi) > np.tan(np.pi/2+observer_lat_rad)):
      SSP_lon = observer_lon_rad +np.pi - np.arcsin(np.sin(psi_alt)*np.sin(azimuth_rad)/np.cos(SSP_lat))
      
    else:
      SSP_lon = observer_lon_rad + np.arcsin(np.sin(psi_alt)*np.sin(azimuth_rad)/np.cos(SSP_lat))
    
    return np.degrees(SSP_lat), np.degrees(SSP_lon)

# Example inputs
observer_latitude = 85  # Replace with your observer latitude
observer_longitude = -10  # Replace with your observer longitude
satellite_altitude = 20000  # Replace with the satellite orbital altitude in kilometers
elevation_angle = 80  # Replace with the elevation angle in degrees
azimuth_angle = -120  # Replace with the azimuth angle in degrees

# Calculate nadir latitude and longitude
nadir_latitude, nadir_longitude = calculate_nadir(observer_latitude, observer_longitude,
                                                  satellite_altitude,
                                                  elevation_angle, azimuth_angle)

# (Optional) Check results using Skyfield API
import skyfield.api as sapi
planets = sapi.load('de430_1850-2150.bsp')
earth_obj = planets['earth']
surface_loc = earth_obj + sapi.wgs84.latlon(observer_latitude* sapi.N, observer_longitude * sapi.E)
obs_obj = earth_obj + sapi.wgs84.latlon(nadir_latitude * sapi.N, nadir_longitude * sapi.E, elevation_m=satellite_altitude*1000)
surf2obs = surface_loc.at(t_UTC).observe(obs_obj).apparent()
verification_LEA_LAA = [surf2obs.altaz()[0].degrees, surf2obs.altaz()[1].degrees]
print(verification_LEA_LAA)

Try it online! (Without Skyfield verification)

For the example above Skyfield calculates angles of 79.989604$^\circ$ and 240.001123$^\circ$ for the elevation and azimuth respectively. If greater accuracy is needed, it could be obtained by using the above as a first guess and then using a bisection type approach with Skyfield to get the matching latitude and longitude.