Calculating the elevation and azimuth angles of satellites with location information obtained from sp3 ephemeris files

Hi, I'm trying to calculate elevation and azimuth angles of satellites for a ground GPS receiver. I downloaded the SP3 file for the day '12-01-2015'. The ephemeris files are xs1846, ys1846, zs1846 are attached.

The ground station location is lat_receiver= 37.0;%degree, lon_receiver= 35.34;%degree.

Are SP3 ephemeris files format ECEF? How can I found the elevation and azimuth angles for that station for all 32 GPS satellite? I need to find the satellittes and the times with an elevation greater than 40 degrees?

I write a code but I am not sure if this is true. Can you help me?

clear all;close all;clc;
load xs1846; load ys1846;load zs1846;


lat_receiver= 37.0;%degree
lon_receiver= 35.34;%degree
height_receiver=0;%meter



for j=1:32
for i=1:96
xss=xs(i,j); yss=ys(i,j);zss=zs(i,j);


  [geosat_ENU_X,geosat_ENU_Y,geosat_ENU_Z] = ecef2enu(xss,yss,zss,...
        lat_receiver,lon_receiver,height_receiver,referenceEllipsoid('wgs84','m'));
    
    elev_receiver(i,j) = rad2deg(asin(geosat_ENU_Z/sqrt(geosat_ENU_X^2+geosat_ENU_Y^2+geosat_ENU_Z^2)));
    azi_NtoE_receiver(i,j) = rad2deg(asin(geosat_ENU_X/sqrt(geosat_ENU_X^2+geosat_ENU_Y^2)));
end
end

.

Answer 1

[Altaïr]

Hey @ysm ony,

The provided code indeed calculates the elevation and azimuth of satellites from their ECEF ephemeris data. The ‘asind’ function outputs a two-quadrant angle, in the range of [-90, 90] degrees. Since azimuth is a four-quadrant angle, I would suggest using the equivalent tan expression for azimuth and implement it using the ‘atan2d’ function in MATLAB. Here is the link to the documentation.

https://www.mathworks.com/help/releases/R2023a/matlab/ref/atan2d.html

However, a simpler and faster approach is available using the ‘ecef2aer’ function. The nested for loop can be replaced with the following line of code:

[azi_NtoE_receiver, elev_receiver, slantRange] = ...
    ecef2aer(xs, ys, zs, lat_receiver, lon_receiver, height_receiver, referenceEllipsoid('wgs84','m'));

You can read more about the ‘ecef2aer’ function in the following documentation page:

https://www.mathworks.com/help/releases/R2023a/map/ref/ecef2aer.html

I hope this resolves your query!