Elevation, Azimuth coherence

[pablotejada]

Hello @brandon-rhodes , I am using skyfield to simuate a satellite LEO constellation as it provides a very nice interface. I get some good results but I wanted to check those with you, as some of them are not clear on the azimuth of a station with respect to a satellite.

The station is Goonhilly "latitude": 50.05, "longitude": -5.18, and the satllite is at altitude 8065km at an inclination of 45 degrees.

Here is the satellite ground path of one period I get using satellite.at(t) and then subpoint():

sat_positions = satellite.at(t)
subpoint = sat_positions.subpoint()

However, when I plot the elevation and azimuth histogram for the station of Goonhilly and the first satellite of the constellation , I get the following:

The figure shows we should point the antenna at northern or completely East/West directions (Azimuths 50-90 deg and 270-310 deg) which does not make sense considering the ground path of the satellite never goes above Goonhilly.

Here is how I define the satellites of the constellation (for this study I only take the first one): Input:

{
      "walker_type": "delta",
      "orbital_planes": 1,
      "sats_per_plane": 5,
      "altitude": 8065,
      "inclination": 45
  }

Code:

MU = 3.986e14 #[m^3s^-2] 
EARTH_R = 6378*1000 #[m]

altitude = altitude*1000
T = 2*np.pi*np.sqrt((EARTH_R + altitude)**3/MU) #[seconds]

epoch = 26800 # has to be the same for all satellites
bstar = 0
ndot = 0.0
nddot = 0
ecco = 0
argpo = 0

inclination = np.deg2rad(inclination)
mo_sats = np.linspace(0,360, sats_per_plane, endpoint=False) #[deg]
no_kozai = 2*np.pi*60/T #[rad/min]
mo_sats = np.linspace(0,360, sats_per_plane, endpoint=False) #[deg]
# RAAN definition for walker type
if walker_type=="star":
    raan_limit_hour = 12
elif walker_type=="delta":
    raan_limit_hour = 24
else:
    raise Exception(f"Wrong walker type: {walker_type} ")
nodeo_planes = np.linspace(0,raan_limit_hour, orbital_planes, endpoint=False) #[hours]
satellites : Dict[str, EarthSatellite] = {}
  for i, nodeo in enumerate(nodeo_planes): # iterating orbital planes
      for j, mo in enumerate(mo_sats): # iterating staelites through the sma orbital planes
          print()
          print(f'RAAN: {np.deg2rad(nodeo)}')
          print(f'Mean Anomaly: {np.deg2rad(mo)}')
          satrec = Satrec()
          satrec.sgp4init(
              WGS84,             # gravity model
              'i',               # 'a' = old AFSPC mode, 'i' = improved mode
              i+j,               # satnum: Satellite number
              26800,             # epoch: days since 1949 December 31 00:00 UT
              bstar,             # bstar: drag coefficient (/earth radii)
              ndot,              # ndot: ballistic coefficient (revs/day)
              nddot,             # nddot: second derivative of mean motion (revs/day^3)
              ecco,              # ecco: eccentricity
              argpo,             # argpo: argument of perigee (radians)
              inclination,       # inclo: inclination (radians)
              np.deg2rad(mo),    # mo: mean anomaly (radians)
              no_kozai,          # no_kozai: mean motion (radians/minute)
              nodeo*2*np.pi/24,  # nodeo: right ascension of ascending node (radians)
          )

          sat_name = f'SAT_shell_{shell_id}_ecco_{ecco}_inclo_{inclination}_RAAN_{nodeo}_mo_{mo}'
          satellites[sat_name] = EarthSatellite.from_satrec(satrec,ts) 

And finally, this is how I compute the elevation (altitudes) and azimuths:

ts = load.timescale()
date_from = datetime(2024,5,1)
minute_range = 6*31*24*60 # x  months 1 min resolution
time_res=1# 1 minute
t = ts.utc(date_from.year, date_from.month, date_from.day, date_from.hour, range(0,minute_range, time_res))
t_str = t.utc_strftime('%Y-%m-%d %H:%M')

sat_0 = satellites[list(satellites)[0]]
station_geo_pos = wgs84.latlon(station['latitude'], station['longitude'])

difference = sat_0 - station_geo_pos

relative_pos_geocentric = difference.at(t)

elev, az, dist = relative_pos_geocentric.altaz()

And the histogram:

elev_threshold = 0
elev_round_value, az_round_value = 5,5
los_df = pd.DataFrame(columns=["Elevation", "Azimuth", "Distance", "LoS", "Elev_rounded", "Az_rounded"])

los_df.Elevation, los_df.Azimuth, los_df.Distance = elev.degrees, az.degrees, dist.km

los_df.Elev_rounded, los_df.Az_rounded = np.ceil(elev.degrees / elev_round_value) * elev_round_value, np.ceil(az.degrees / az_round_value) * az_round_value
los_df.loc[los_df.Elevation > elev_threshold, "LoS"] = 1.0
los_df.loc[los_df.Elevation <=elev_threshold, "LoS"] = 0.0

los_df = los_df.set_index(pd.to_datetime(t_str))

import matplotlib.pyplot as plt

visisbility_df = los_df.loc[los_df.Elevation > elev_threshold]
elev_bin_limits = np.arange(0,90,5)
az_bin_limits = np.arange(0,360,5)

hist, x_edges, y_edges = np.histogram2d(visisbility_df.Elevation, visisbility_df.Azimuth, bins=[elev_bin_limits, az_bin_limits])

plt.figure(figsize=(10, 8))
plt.imshow(hist, origin='lower', extent=[az_bin_limits.min(), az_bin_limits.max(), elev_bin_limits.min(), elev_bin_limits.max()], aspect='auto', cmap='viridis')
plt.colorbar(label='Frequency')
plt.ylabel('Elevation')
plt.xlabel('Azimuth')
plt.title(f'{station_study} 2D Histogram of Elevation and Azimuth')
plt.show()

I have a feeling I am donig something wrong, maybe the coordinate system is not consistent but I have not found what may be causing the issue. Any help is much appreciated!

Tanks in advance!

[brandon-rhodes]

The figure shows we should point the antenna at northern…

I'm assuming that you meant to type "southern", because of the glow around 270° on the histogram.

the satllite is at altitude 8065km…

Do you mean 8,000 km above the ground? That's more than the radius of the Earth, so I would certainly expect it to be visible from Goonhilly. If you mean instead that its orbit's radius is 8,000 km, then that would still put it at more than 1,500 km above the Earth's surface, and thus plausibly still visible from Goonhilly. I would have to draw an actual triangle with those dimensions to test whether it would really be below the horizon. Could you share the math you're doing that assures us that it doesn't come over Goonhilly's horizon?

[pablotejada]

Thanks @brandon-rhodes for the quick feedback!

Regarding the first comment, the most frequent values (in yellow) of azimuth are 60-80 degrees and 270-300 degrees. From documentation: "the azimuth Angle measures east along the horizon from geographic north (so 0 degrees means north, 90 is east, 180 is south, and 270 is west)." From this histogram then, an antenna at Goonhilly should be frequently pointed at north-east, east, north-west and west azimuth directions. That seems extrange given the satellite ground path of the first figure.

Regarding the the second comment, the satellite is at an altitude (above the ground) of 8065 km. To compute the period of the orbit I use the following formula: T = 2*np.pi*np.sqrt((EARTH_R + altitude)**3/MU) #[seconds] which allows to define no_kozai = 2*np.pi*60/T #[rad/min]. The rest of parameters are defined as 0:

bstar = 0
ndot = 0.0
nddot = 0
ecco = 0
argpo = 0

Now, I can see indeed that the satellite will rise above the horizon at Goonhilly, also the histogram shows that by the points where elevation values are above 0 degrees. However, what I can't seem to understand, is why the azimuth takes values on those ranges of 60-80 degrees and 270-300 degrees. Looking at the ground path of the satellite in the first figure, we would expect to see the satellite from Goonhilly always at southern directions, hence azimuths between [100, 250] and definetly never below 90 or above 270 degrees.

Again, any comments are much appreciated, maybe I am just interpreting wrong the results.

Thanks in advance and have a great rest of the day,

Pablo

[brandon-rhodes]

The relative bearing between any two points at the same latitude like, say, 50°, is always north-of-straight-east and north-of-straight-west, because the great circle that's the shortest distance between loops up into higher latitudes to reach the other point. Could something similar be happening here?

https://www.oag.com/blog/great-circle-routes-flight-paths

I myself have been confused by that when trying to work out which way the illuminated edge of the Moon should point! Here's a paper explaining that situation:

https://www.seas.upenn.edu/~amyers/MoonPaperOnline.pdf

[pablotejada]

Hello @brandon-rhodes. Thanks again for your input. I had a look at the great circles and the paper and it makes sense up to a certain point. I implemented the elevation and azimuth equations and I get interesting results:

For Elevation I get 0.0* errors which are negligible. For Azimuth I get big errors. I checked the equations and I think Skyfield altaz() does not correct the Azimuth considering the relative position between the sub-satellite point and the ground station.

I took the equations from Section 2.1.6.3 of the book Satellite Communications Systems Systems, Techniques and Technology (Gerard Maral, Michel Bousquet, Zhili Sun) where it states that the actual Azimuth has to be corrected according to the following table:

For reference:

Where L is the relative longitude between the satellite and the station.

Coding the equations from the book I get these results (Code attached below):

Error on Elevation [deg]:
50        0.015515
51        0.016256
52        0.016919
53        0.017508
54        0.018028
            ...   
266906    0.017314
266907    0.016879
266908    0.016413
266909    0.015916
266910    0.015387
Name: Elevation, Length: 10236, dtype: float64

Error on Skyfield Azimuth [deg]:
50         2.259384
51         5.523434
52         8.688221
53        11.753754
54        14.720581
            ...    
266906    64.173184
266907    65.175783
266908    66.145115
266909    67.082637
266910    67.989737
Name: Azimuth, Length: 10236, dtype: float64

Error on Auxilirar Azimuth [deg]:
50        0.024748
51        0.023157
52        0.021638
53        0.020197
54        0.018837
            ...   
266906    0.009150
266907    0.009470
266908    0.009816
266909    0.010186
266910    0.010581
Name: Azimuth, Length: 10236, dtype: float64

From the results, looks like Skyfield's azimuth relates to the book's auxiliar azimuth a and we might be missing the correction considering the cases on the Table 2.3.

I recomputed the histogram and results are the following:

Code to compute Elevation and Azimuth correcting by relative position considering Table 2.3:

def compute_elev_az(sat_lat:float, sat_lon:float, h_sat: float, ground_point_lat:float, ground_point_lon:float):
    """
    Computes Elevation and Azimuth from a ground point to a satellite 
    given the lat, lon nad altitude of the satelite and lat lon of the ground point.

    Input:  
        - sat_lat, sat_lon, h_sat (float): lat, lon [deg] and altitude [m] of the satellite 
        - ground_lat, ground_lon (float): lat, lon of the station on erath surface [deg]
    Returns:
        - E (float): Elevation angle [deg]
        - A (float): Azimuth angle [deg]
        - a (floar): Auxiliar Azimuth [deg]
    """

    sat_lat = np.deg2rad(sat_lat) 
    sat_lon = np.deg2rad(sat_lon) 
    ground_point_lat = np.deg2rad(ground_point_lat) 
    ground_point_lon = np.deg2rad(ground_point_lon)    

    # ELEVATION
    L = abs(ground_point_lon - sat_lon)
    phi = np.arccos(
        np.cos(L)*np.cos(sat_lat)*np.cos(ground_point_lat) + np.sin(sat_lat)*np.sin(ground_point_lat)
    )
    r = EARTH_R + h_sat
    R = np.sqrt(EARTH_R**2 + r**2 - 2*EARTH_R*r*np.cos(phi))
    E = np.arccos(r/R * np.sin(phi)) # Elevation [rad]
    E = np.rad2deg(E)

    # AZIMUTH
    #assert any(phi>0)
    phi = abs(phi)
    a = np.arcsin(
        np.sin(L)*np.cos(sat_lat)/np.sin(phi)
    )
    a = np.rad2deg(a)

    # Position of the sub-satellite point T with respect to point P
    A = np.full_like(sat_lat, np.nan)  # Initialize with NaNs
    SE_condition = np.logical_and(sat_lat < ground_point_lat, sat_lon > ground_point_lon)
    NE_condition = np.logical_and(sat_lat > ground_point_lat, sat_lon > ground_point_lon)
    SW_condition = np.logical_and(sat_lat < ground_point_lat, sat_lon < ground_point_lon)
    NW_condition = np.logical_and(sat_lat > ground_point_lat, sat_lon < ground_point_lon)

    # Check conditions are disjoint and complete
    sum_condition = SE_condition.astype(int) + NE_condition.astype(int) + SW_condition.astype(int) + NW_condition.astype(int)  
    is_only_one_true = np.all(sum_condition == 1)
    assert is_only_one_true == True

    A = np.where(SE_condition, 180 - a, A)
    A = np.where(NE_condition, a, A)
    A = np.where(SW_condition, 180 + a, A)
    A = np.where(NW_condition, 360 - a, A)

    return E, A, a

The fact that we get such close results for E and a is recomforting, but in the case of the actual Azimuth we might be missing the geometric corrections.

If you have any comments regarding this I'd love to know your opinion!

Thanks again and have a great week!

Pablo

[brandon-rhodes]

I cannot run your example code, alas, because it doesn't include its import statements, nor am I set up to call the routine with exactly the same arguments you did. Could you provide a full working Python script for me to run? Maybe you could just focus on one moment in time, and show Skyfield's number, then how you compute the number you wish Skyfield were returning. Thanks!

[pablotejada]

Hello again @brandon-rhodes ,

Here's the code to be run. It generates three figures, the last one is a plot comparing the results of Elevation and Azimuth. I have to say though, that the Elevation with these custom equations is like an absolute value of Skyfield values. In the azimuth the differences are bigger.

If you have the source from where you got the equations I would be happy to check them as well.

Thanks again and have a great weekend!

# %%
from skyfield.api import load, wgs84, EarthSatellite
from sgp4.api import Satrec, WGS84
import numpy as np
from datetime import datetime

# %%
ts = load.timescale()
date_from = datetime(2024,5,1)
minute_range = 6*31*24*60 # x  months 1 min resolution
time_res=1# 1 minute
t = ts.utc(date_from.year, date_from.month, date_from.day, date_from.hour, range(0,minute_range, time_res))
t_str = t.utc_strftime('%Y-%m-%d %H:%M')
MU = 3.986e14 #[m^3 s^-2] 
EARTH_R = 6378*1000 #[m]

# %%
def compute_elev_az(sat_lat:float, sat_lon:float, h_sat: float, ground_point_lat:float, ground_point_lon:float):
    """
    Computes Elevation and Azimuth from a ground point to a satellite 
    given the lat, lon and altitude of the satelite and lat lon of the ground point.

    Input:  
        - sat_lat, sat_lon, h_sat (float): lat, lon [deg] and altitude [m] of the satellite 
        - ground_lat, ground_lon (float): lat, lon of the station on erath surface [deg]
    Returns:
        - E (float): Elevation angle [deg]
        - A (float): Azimuth angle [deg]
        - a (floar): Auxiliar Azimuth [deg]
    """

    sat_lat = np.deg2rad(sat_lat) 
    sat_lon = np.deg2rad(sat_lon) 
    ground_point_lat = np.deg2rad(ground_point_lat) 
    ground_point_lon = np.deg2rad(ground_point_lon)    

    # ELEVATION
    L = abs(ground_point_lon - sat_lon)
    phi = np.arccos(
        np.cos(L)*np.cos(sat_lat)*np.cos(ground_point_lat) + np.sin(sat_lat)*np.sin(ground_point_lat)
    )

    r = EARTH_R + h_sat
    R = np.sqrt(EARTH_R**2 + r**2 - 2*EARTH_R*r*np.cos(phi))
    E = np.arccos(r/R * np.sin(phi)) # Elevation [rad]
    E = np.rad2deg(E)

    # AZIMUTH
    #assert any(phi>0)

    a = np.arcsin(
        np.sin(L)*np.cos(sat_lat)/np.sin(phi)
    )
    a = np.rad2deg(a)

    # Position of the sub-satellite point T with respect to point P
    A = np.full_like(sat_lat, np.nan)  # Initialize with NaNs
    SE_condition = np.logical_and(sat_lat < ground_point_lat, sat_lon > ground_point_lon)
    NE_condition = np.logical_and(sat_lat > ground_point_lat, sat_lon > ground_point_lon)
    SW_condition = np.logical_and(sat_lat < ground_point_lat, sat_lon < ground_point_lon)
    NW_condition = np.logical_and(sat_lat > ground_point_lat, sat_lon < ground_point_lon)

    # Check conditions are disjoint and complete
    sum_condition = SE_condition.astype(int) + NE_condition.astype(int) + SW_condition.astype(int) + NW_condition.astype(int)  
    is_only_one_true = np.all(sum_condition == 1)
    assert is_only_one_true == True

    A = np.where(SE_condition, 180 - a, A)
    A = np.where(NE_condition, a, A)
    A = np.where(SW_condition, 180 + a, A)
    A = np.where(NW_condition, 360 - a, A)

    return E, A, a

# %%

altitude = 8065 # [km]
altitude = altitude*1000  #[m]
T = 2*np.pi*np.sqrt((EARTH_R + altitude)**3/MU) #[seconds]

epoch = 26800 # has to be the same for all satellites
bstar = 0
ndot = 0.0
nddot = 0
ecco = 0
argpo = 0
nodeo = 0
mo = 0
inclination = np.deg2rad(45)
no_kozai = 2*np.pi*60/T #[rad/min]

ts = load.timescale()
satrec = Satrec()
satrec.sgp4init(
                WGS84,             # gravity model
                'i',               # 'a' = old AFSPC mode, 'i' = improved mode
                0,                 # satnum: Satellite number
                26800,             # epoch: days since 1949 December 31 00:00 UT
                bstar,             # bstar: drag coefficient (/earth radii)
                ndot,              # ndot: ballistic coefficient (revs/day)
                nddot,             # nddot: second derivative of mean motion (revs/day^3)
                ecco,              # ecco: eccentricity
                argpo,             # argpo: argument of perigee (radians)
                inclination,       # inclo: inclination (radians)
                np.deg2rad(mo),    # mo: mean anomaly (radians)
                no_kozai,          # no_kozai: mean motion (radians/minute)
                nodeo*2*np.pi/24,  # nodeo: right ascension of ascending node (radians)
            )

sat_test = EarthSatellite.from_satrec(satrec,ts) 

# %%
T/60/60

# %%
stationn_lat = 50.05
station_lon =  -5.18
station_geo_pos = wgs84.latlon(stationn_lat, station_lon)
station_geo_pos

# %%
sat_0_geocentric = sat_test.at(t)
station_geocentric = station_geo_pos.at(t)
rel_pos = sat_0_geocentric - station_geocentric
sat_geodetic_pos = wgs84.geographic_position_of(sat_0_geocentric)

# %%
np.sqrt(sum(sat_0_geocentric.position.km**2)) - EARTH_R/1000

# %%
elev, az, dist = rel_pos.altaz()

# %%

E, A, a = compute_elev_az(
    sat_lat=sat_geodetic_pos.latitude.degrees, 
    sat_lon=sat_geodetic_pos.longitude.degrees, 
    h_sat=altitude, 
    ground_point_lat= stationn_lat,
    ground_point_lon= station_lon
)

# %%
import matplotlib.pyplot as plt
az_bin_limits = np.arange(0,360,5)

elev_bin_limits = np.arange(0,90,5)
hist, x_edges, y_edges = np.histogram2d(
    elev.degrees,
    az.degrees,
    bins=[elev_bin_limits, az_bin_limits]
)

plt.figure(figsize=(10, 8))
plt.imshow(hist, origin='lower', extent=[az_bin_limits.min(), az_bin_limits.max(), elev_bin_limits.min(), elev_bin_limits.max()], aspect='auto', cmap='viridis')
plt.colorbar(label='Frequency')
plt.ylabel('Elevation')
plt.xlabel('Azimuth')
plt.title('2D Histogram of Elevation and Azimuth Skyfield')
plt.show()

# %%
import matplotlib.pyplot as plt
az_bin_limits = np.arange(0,360,5)

elev_bin_limits = np.arange(0,90,5)
hist, x_edges, y_edges = np.histogram2d(
    E,
    A,
    bins=[elev_bin_limits, az_bin_limits]
)

plt.figure(figsize=(10, 8))
plt.imshow(hist, origin='lower', extent=[az_bin_limits.min(), az_bin_limits.max(), elev_bin_limits.min(), elev_bin_limits.max()], aspect='auto', cmap='viridis')
plt.colorbar(label='Frequency')
plt.ylabel('Elevation')
plt.xlabel('Azimuth')
plt.title('2D Histogram of Elevation and Azimuth Skyfield')
plt.show()

# %%
import pandas as pd
elev_az_df = pd.DataFrame(
    data= {
        "Elevation_Skyfield": elev.degrees,
        "Elevation_custom": E,
        "Azimuth_Skyfield": az.degrees,
        "Azimuth_aux": a,
        "Azimuth_corrected": A,
        "Error_Elevation": abs(elev.degrees - E),
        "Error_Azimuth": abs(az.degrees - A),
        "Sat_Latitude": sat_geodetic_pos.latitude.degrees,
        "Sat_Longitude": sat_geodetic_pos.longitude.degrees
    }
)
elev_az_df

# %%
elev_az_df.loc[(elev_az_df.Elevation_Skyfield < 5) & (elev_az_df.Elevation_Skyfield > -5), ["Elevation_Skyfield", "Elevation_custom", "Sat_Latitude", "Sat_Longitude"]].iloc[10:20]

# %%
import matplotlib.pyplot as plt

indexes = slice(1,500)

df = elev_az_df.loc[elev_az_df.Elevation_Skyfield > 0].iloc[indexes]
# df = elev_az_df[indexes]

# Create a subplot with 1 row and 2 columns
fig, ax = plt.subplots(2, 1)

# Plot the first line

ax[0].plot(df.Elevation_Skyfield, label="Skyfield")
ax[0].plot(df.Elevation_custom, label="Custom")
ax[0].set_title('Elevation')
ax[0].set_ylabel('[deg]')
ax[0].legend()

# Plot the second line
ax[1].plot(df.Azimuth_Skyfield, label="Skyfield")
ax[1].plot(df.Azimuth_corrected, label="Custom_Corrected")
ax[1].plot(df.Azimuth_aux, label="Custom_Aux")
ax[1].set_title('Azimuth')
ax[1].set_ylabel('[deg]')
ax[1].legend()

# Adjust layout to prevent overlap
plt.tight_layout()

# Display the plot
plt.show()

# %%
error_elev = max(df.Error_Elevation)
error_az = max(df.Error_Azimuth)

print("Max erorr elevation: ", error_elev)
print("Max erorr azimuth: ", error_az)

[EndlessDex]

So let me start by saying that Skyfield is VERY accurate in its satellite propagation and does NOT need any post-processing/modifications to yield exact satellite look angles.

Your initial question was regarding why you saw the hotspots you did. I believe this can be most clearly answered by looking at a polar plot of the look angles.

Code

```python import skyfield.api as sf from skyfield.constants import ERAD import sgp4.api as sgp4 import numpy as np import matplotlib.pyplot as plt def build_satellite(epoch_time, altitude_km, inclination_deg): altitude_m = altitude_km * 1000.0 MU = 3.986004418e14 # m^3 * s^−2 period = 2.0 * np.pi * np.sqrt((ERAD + altitude_m) ** 3 / MU) # [seconds] TLE_EPOCH_OFFSET = 2433281.5 tle_epoch = epoch_time.to_astropy().to_value("jd") - TLE_EPOCH_OFFSET no_kozai = (2.0 * np.pi) / (period / 60.0) # [rad/min] ts = sf.load.timescale() satrec = sgp4.Satrec() satrec.sgp4init( sgp4.WGS84, # gravity model "i", # 'a' = old AFSPC mode, 'i' = improved mode 0, # satnum: Satellite number tle_epoch, # epoch: days since 1949 December 31 00:00 UT 0.0, # bstar: drag coefficient (/earth radii) 0.0, # ndot: ballistic coefficient (revs/day) 0.0, # nddot: second derivative of mean motion (revs/day^3) 0.0, # ecco: eccentricity 0.0, # argpo: argument of perigee (radians) np.deg2rad(inclination_deg), # inclo: inclination (radians) 0.0, # mo: mean anomaly (radians) no_kozai, # no_kozai: mean motion (radians/minute) 0.0, # nodeo: right ascension of ascending node (radians) ) return sf.EarthSatellite.from_satrec(satrec, ts) def propagate_satellite(time_range, satellite, station): look_vector = satellite - station alt, az, _ = look_vector.at(time_range).altaz() # Only care about when above horizon az = az.degrees alt = alt.degrees clipped_az = az[alt > 0] clipped_alt = alt[alt > 0] sat_vector = satellite.at(time_range) lat, lon = sf.wgs84.latlon_of(sat_vector) lat = lat.degrees lon = lon.degrees clipped_lat = lat[alt > 0] clipped_lon = lon[alt > 0] return (clipped_alt, clipped_az, clipped_lat, clipped_lon) ts = sf.load.timescale() date_from = ts.utc(2024, 5, 1) range_days = 1 time_res = 1 # 1 minute date_from_dt = date_from.utc_datetime() t_range = ts.utc( date_from_dt.year, date_from_dt.month, date_from_dt.day, date_from_dt.hour, range(0, range_days * 24 * 60, time_res), ) satellite = build_satellite(date_from, 8065, 45) station = sf.wgs84.latlon( latitude_degrees=50.05, longitude_degrees=-5.18 ) alt, az, lat, lon = propagate_satellite(t_range, satellite, station) # Make a plot to show local apparent trajectories fig, ax = plt.subplots(subplot_kw={"projection": "polar"}) ax.set_rlim(bottom=90, top=0) ax.scatter(np.deg2rad(az), alt, alpha=0.5) ax.set_theta_zero_location("N") ax.set_theta_direction(-1) ax.grid(True) plt.title(f"Look angles over {range_days} day(s)") plt.show() # Make a plot to show local apparent trajectories fig, ax = plt.subplots(subplot_kw={"projection": "polar"}) ax.set_rlim(bottom=90, top=0) ax.scatter(np.deg2rad(lon), lat, alpha=0.5) ax.scatter(np.deg2rad(station.longitude.degrees), station.latitude.degrees) ax.set_theta_zero_location("N") ax.set_theta_direction(-1) ax.grid(True) plt.title(f"Lat/Lon over {range_days} day(s)") plt.show() ```

As you can see, in one day there are three periods of visibility (passes). Each pass starts at the horizon, peaks, and then sets at the horizon. With up being north and the center of the plot being zenith.

Here is a plot of the same passes but in the geodetic (lat/lon) frame. The orange dot is your station. This is from the perspective of looking down on the earth from the north pole with The equator being the outer border and up being the prime meridian.

I believe that given these plots it is evident why you find yourself looking due east and west regularly. This is the effect of "great circle" that Brandon mentioned to you in his second reply to you. The ground plot you initially displayed is incredibly misleading for the understanding of satellite look angles. 2D polar plots much more closely match the reality of the 3D earth and orbit paths. Especially because you are orbiting at almost 2x the earths radius which means you can see the satellite when its ground track is over 90deg longitude away from your station.

Here's some more plots for larger date ranges where you can really see the nodes forming.

 

And even longer!!

 

[EndlessDex]

I am going to be honest and say that I did not really read your replies after you started trying to correct the azimuth. You were WAY in the weeds and applying some stuff incorrectly.

Please let me know if my post doesn't clear up everything for you.

[pablotejada]

Hi @EndlessDex,

Thank you for the detailed explanation. I totally agree that skyfield is very accurate on the calculations and it has been very helpful so far. On my previous comments, I just tried to understand mathematically the great circle concept between a station and a satellite that Brandon suggested by comparing skyfield results to the analytical equations found in text books. The Elevation seemed to match perfectly but not the azimuth, but seeing the polar plots you appended makes total sense now.

Thanks again for the responses and for the great work!

[brandon-rhodes]

@EndlessDex — Thanks for jumping in and providing some diagrams! After trying to describe the great-circle situation in my first couple of comments, I didn't have the free cycles this summer to read through the formulae in the screenshot above and try to pursue this conversation further. Your idea of following up with explicit diagrams was a great idea, and it looks like the original question is now resolved to @pablotejada’s satisfaction. Thanks!