3
$\begingroup$

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:

but the answers are not consistent.

I would need a reference if it is available.

$\endgroup$
1
  • $\begingroup$ As mentioned 8n the answer below, you need at least one more piece of information, like the distance to then satellite. If you don't want to guess, you can use three observations to compute the orbit, and get the distance from that. $\endgroup$ Dec 5, 2023 at 23:22

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ The observer is not on the surface but on zero height. $\endgroup$
    – Uwe
    Dec 7, 2023 at 14:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.