In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import georinex as gr
import xarray as xr
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
import pathlib
import scipy.signal
import gzip
import scipy.stats
plt.rcParams.update({'figure.max_open_warning': 0})
In [2]:
mu = 3.986004418e14
c = 299792458
pi = 3.1415926535898
def largest_eccentricity(data, const, ignore_igsos = True):
ignore = ['G21', 'E14', 'E18'] # G21 has really large ecc
if ignore_igsos:
ignore.extend(['C06', 'C07', 'C08', 'C09', 'C10', 'C12', 'C13', 'C16'])
svs = np.array([sv.item() for sv in data.coords['sv'] if sv.item()[0] == const and '_' not in sv.item() and not sv.item() in ignore])
eccs = np.array([data.sel(sv = sv).mean().item() for sv in svs])
return svs[np.argmax(eccs)]
def ecc(nav):
return nav['Eccentricity'].mean().item()
def sqrtA(nav):
return nav['sqrtA'].mean().item()
def n0(nav):
return np.sqrt(mu)/sqrtA(nav)**3
def relativistic_bias(nav):
return 2*np.sqrt(mu)/c**2*sqrtA(nav)*ecc(nav)
def relativistic_drift(nav):
return relativistic_bias(nav) * n0(nav)
def data_time(data):
return (data.coords['time'] - data.coords['time'][0]).values.astype('float') * 1e-9
def detrend(data, detrend_degree):
x = data.dropna('time')
t = data_time(x)
p_fit = np.polyfit(t, x, detrend_degree)
return x - np.polyval(p_fit, t)
def fourier_analysis(data, n0, detrend_degree = None, **kwargs):
freqs = np.linspace(0, 4.5, 10000)
amplitudes = np.empty_like(freqs)
x = data.dropna('time')
t = data_time(x)
if detrend_degree is None:
y = x
else:
y = detrend(x, detrend_degree)
for j, f in enumerate(freqs):
amplitudes[j] = 2*np.abs((y * np.exp(-1j * f * n0 * t)).mean().item())
plt.plot(freqs, amplitudes, **kwargs)
def fourier_extract(data):
detrends = {'SVclockBias' : 4, 'SVclockDrift' : 0}
data = data.dropna('time')
t = data_time(data)
M0 = data['M0'][0].item()
fourier = np.empty(5, dtype = 'complex')
lo = 2 * np.exp(-1j * (n0(data) * t + M0))
for j, series in enumerate(['SVclockBias', 'SVclockDrift']):
x = detrend(data[series], detrends[series])
fourier[j] = (x * lo).mean().item()
fourier[2] = relativistic_bias(data)
fourier[3] = relativistic_drift(data)
fourier[4] = n0(data)
return fourier
def sp3_load(path, svn):
psvn = bytes(f'P{svn}', encoding = 'ascii')
with gzip.open(path) as sp3:
ls = sp3.readlines()
return np.array([float(l.strip(b'\n').split()[-1]) for l in ls if l.startswith(psvn)][:-1])
def load_sp3_clk(svn):
sp3s = pathlib.Path('sp3').glob('COD0MGXFIN_*_01D_05M_ORB.SP3.gz')
sp3_start = np.datetime64('2019-11-01T00:00:00')
clk_us = np.concatenate([sp3_load(sp3, svn) for sp3 in sp3s])
t_clk = sp3_start + np.timedelta64(1, 'ns') * np.arange(clk_us.size) * 5 * 60 * 1e9
return xr.DataArray(clk_us * 1e-6, coords = [t_clk], dims = ['time'])
Determine satellites with largest eccentricity in each constellation.
In [3]:
eccs = gr.load('brdc/BRDC00IGS_R_20193050000_01D_MN.rnx.gz')['Eccentricity']
For GPS we ignore G21, which has the largest eccentricity, for reasons explained below.
In [4]:
largest_eccentricity(eccs, 'G')
Out[4]:
In [5]:
largest_eccentricity(eccs, 'E')
Out[5]:
For Beidou, several IGSOs have much more eccentricity than any of the MEOs. We pick separately the IGSO with largest eccentricity and the MEO with largest eccentricity.
In [6]:
largest_eccentricity(eccs, 'C', ignore_igsos = False)
Out[6]:
In [7]:
largest_eccentricity(eccs, 'C')
Out[7]:
Comparison between the eccentricities of G21 (largest) and G02 (next largest).
In [8]:
eccs.sel(sv = 'G21').mean().item()
Out[8]:
In [9]:
eccs.sel(sv = 'G02').mean().item()
Out[9]:
Uncomment to regenerate .nc files from RINEX files. It takes around 6 minutes per satellite to run on my machine.
In [10]:
#%%time
#keep_variables = ['SVclockBias', 'SVclockDrift', 'SVclockDriftRate', 'sqrtA', 'Eccentricity', 'M0']
#svns = ['G02', 'G21', 'E31', 'C06', 'C11', 'R16'] + \
# ['E12', 'E19', 'E14', 'E18', 'E26', 'E24', 'E30', 'E07', 'E08', 'E09', 'E01', 'E02', 'E03',\
# 'E04', 'E05', 'E21', 'E27', 'E31', 'E36', 'E13', 'E15', 'E33'] # E11 omitted since it was out of service, E25 also out of service
#for svn in svns:
# rnxs = pathlib.Path('brdc').glob('BRDC00IGS_R_*.rnx.gz')
# data = xr.concat((gr.load(rnx)[keep_variables].sel(sv = svn).dropna('time', 'all') for rnx in rnxs), dim = 'time')
# data.to_netcdf(f'nc/{svn}.nc', group = 'NAV')
The orbital data for GLONASS is not taken from the ephemeris, but rather from the almanac on 2019-12-01.
In [11]:
data_gps = gr.load('nc/G02.nc')
data_gps_G21 = gr.load('nc/G21.nc')
data_gal = gr.load('nc/E31.nc')
data_bei_igso = gr.load('nc/C06.nc')
data_bei_meo = gr.load('nc/C11.nc')
data_glo = gr.load('nc/R16.nc')
clk_gps = load_sp3_clk('G02')
clk_gal = load_sp3_clk('E31')
clk_bei_igso = load_sp3_clk('C06')
clk_bei_meo = load_sp3_clk('C11')
clk_bei_meo[clk_bei_meo > 0.1] = np.nan # clean some spurious values
clk_glo = load_sp3_clk('R16')
n0_glo = 2*pi/40544.293 # not available from broadcast keplerian parameters
Compute relativistic amplitude terms for clock bias $$2\frac{\sqrt{\mu}}{c^2}e\sqrt{A}$$ and clock drift $$2\frac{\sqrt{\mu}}{c^2}e\sqrt{A}n_0,$$ since the relativistic clock correction is $$-2\frac{\sqrt{\mu}}{c^2}e\sqrt{A}\sin(E)$$ and $$E \approx M \approx M_0 + n(t-t_0) \approx M_0 + n_0(t-t_0).$$
In [12]:
relativistic_bias(data_gal), relativistic_drift(data_gal)
Out[12]:
In [13]:
relativistic_bias(data_gps), relativistic_drift(data_gps)
Out[13]:
In [14]:
relativistic_bias(data_bei_igso), relativistic_drift(data_bei_igso)
Out[14]:
In [15]:
relativistic_bias(data_bei_meo), relativistic_drift(data_bei_meo)
Out[15]:
In [16]:
# relativistic clock bias for R16
e_glo = 0.00272
sqrtA_glo = (np.sqrt(mu)/n0_glo)**(1/3)
2*np.sqrt(mu)/c**2*sqrtA_glo*e_glo
Out[16]:
In [17]:
plt.figure(figsize = (12,6), facecolor = 'w')
detrend(data_gal['SVclockDrift'], 0).plot(marker = '.')
detrend(data_gps['SVclockDrift'], 0).plot(marker = '.')
plt.title('Broadcast clock drift (average removed)')
plt.xlabel('UTC time')
plt.ylabel('Clock drift (s/s)')
plt.legend(['E31', 'G02']);
The cell bellow shows that the GPS clock drift is constant except for a few jumps and some tiny jumps.
In [18]:
np.sort(data_gps['SVclockDrift'].dropna('time').diff('time').values)[::-1][:50]
Out[18]:
In [19]:
plt.figure(figsize = (12,6), facecolor = 'w')
detrend(data_bei_igso['SVclockDrift'], 1).plot(marker = '.')
detrend(data_bei_meo['SVclockDrift'], 0).plot(marker = '.')
plt.title('Broadcast clock drift (linear trend/average removed)')
plt.xlabel('UTC time')
plt.ylabel('Clock drift (s/s)')
plt.legend(['C06 (IGSO)', 'C11 (MEO)']);
In [20]:
plt.figure(figsize = (12,6), facecolor = 'w')
detrend(data_gal['SVclockBias'], 4).plot()
detrend(data_gps['SVclockBias'], 4).plot()
plt.title('Broadcast clock bias (4th order polynomial removed)')
plt.xlabel('UTC time')
plt.ylabel('Clock bias (s)')
plt.legend(['E31', 'G02']);
In [21]:
plt.figure(figsize = (12,6), facecolor = 'w')
detrend(data_bei_igso['SVclockBias'], 4).plot()
detrend(data_bei_meo['SVclockBias'], 4).plot()
plt.title('Broadcast clock bias (4th order polynomial removed)')
plt.xlabel('UTC time')
plt.ylabel('Clock bias (s)')
plt.legend(['C06 (IGSO)', 'C11 (MEO)']);
In [22]:
plt.figure(figsize = (12,6), facecolor = 'w')
detrend(data_glo['SVclockBias'], 4).plot()
plt.title('Broadcast clock bias (4th order polynomial removed)')
plt.xlabel('UTC time')
plt.ylabel('Clock bias (s)')
plt.legend(['R16']);
In [23]:
plt.figure(figsize = (12,7), facecolor = 'w')
fourier_analysis(data_gal['SVclockDrift'], n0(data_gal), 0)
fourier_analysis(data_gps['SVclockDrift'], n0(data_gps), 0)
plt.title('Broadcast clock drift spectrum (average removed)')
plt.xlabel('Frequency (1/revolution)')
plt.ylabel('Amplitude (s/s)')
plt.legend(['E31', 'G02']);
In [24]:
plt.figure(figsize = (12,7), facecolor = 'w')
fourier_analysis(data_bei_igso['SVclockDrift'], n0(data_bei_igso), 1)
fourier_analysis(data_bei_meo['SVclockDrift'], n0(data_bei_meo), 0)
plt.title('Broadcast clock drift spectrum (linear trend/average removed)')
plt.xlabel('Frequency (1/revolution)')
plt.ylabel('Amplitude (s/s)')
plt.legend(['C06 (IGSO)', 'C11 (MEO)']);
In [25]:
plt.figure(figsize = (12,7), facecolor = 'w')
fourier_analysis(data_gal['SVclockBias'], n0(data_gal), 4)
fourier_analysis(data_gps['SVclockBias'], n0(data_gps), 4)
fourier_analysis(data_bei_meo['SVclockBias'], n0(data_bei_meo), 4, alpha = 0.8)
plt.title('Broadcast clock bias spectrum (4th order polynomial removed)')
plt.xlabel('Frequency (1/revolution)')
plt.ylabel('Amplitude (s)')
plt.legend(['E31', 'G02', 'C11 (MEO)']);
In [26]:
plt.figure(figsize = (12,7), facecolor = 'w')
fourier_analysis(data_gal['SVclockBias'], n0(data_gal), 4)
fourier_analysis(clk_gal, n0(data_gal), 4)
plt.title('Galileo clock bias spectrum (4th order polynomial removed)')
plt.xlabel('Frequency (1/revolution)')
plt.ylabel('Amplitude (s)')
plt.legend(['E31 (broadcast)', 'E31 (SP3)']);
In [27]:
plt.figure(figsize = (12,7), facecolor = 'w')
fourier_analysis(data_gps['SVclockBias'], n0(data_gps), 4)
fourier_analysis(clk_gps, n0(data_gps), 4)
plt.title('GPS clock bias spectrum (4th order polynomial removed)')
plt.xlabel('Frequency (1/revolution)')
plt.ylabel('Amplitude (s)')
plt.legend(['G02 (broadcast)', 'G02 (SP3)']);
In [28]:
plt.figure(figsize = (12,7), facecolor = 'w')
fourier_analysis(data_gps_G21['SVclockBias'], n0(data_gps_G21), 4)
fourier_analysis(load_sp3_clk('G21'), n0(data_gps_G21), 4)
plt.title('GPS clock bias spectrum (4th order polynomial removed)')
plt.xlabel('Frequency (1/revolution)')
plt.ylabel('Amplitude (s)')
plt.legend(['G21 (broadcast)', 'G21 (SP3)']);