In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.dates as mdates
import xarray as xr

from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

samp_rate = 100
n_avg = 1000
T_avg = n_avg/samp_rate

In [3]:
def load_timeseries(filename, freq_tune = 65, freq_correct = 0, skip = 0):
    phase = np.fromfile(filename, dtype = 'float32')
    t_start = np.datetime64(filename.split('_')[-1].rstrip('.f32'))
    phase_correct = (freq_correct/samp_rate * 2 * np.pi * np.arange(phase.size)) % (2*np.pi)\
        if freq_correct != 0 else 0
    t = t_start + np.arange(0, phase.size, n_avg) * np.timedelta64(int(1e9/samp_rate), 'ns')
    f = np.diff(np.unwrap(phase + phase_correct)[::n_avg])/(2*np.pi*T_avg) + freq_tune
    return xr.DataArray(np.concatenate((f[skip:], [np.nan])), # note we insert a nan at the end to break up the series
                        coords = {'time' : t[skip:]}, dims = ('time'))

In [4]:
f = xr.concat((load_timeseries('phase_bpsk_2020-04-09T15:31:07.376113.f32', freq_correct = -25),
               load_timeseries('phase_bpsk_2020-04-10T08:59:32.991039.f32', skip = 7),
               load_timeseries('phase_bpsk_2020-04-11T08:05:47.727518.f32', skip = 8),
               load_timeseries('phase_bpsk_2020-04-11T19:59:53.454352.f32', skip = 7),
               load_timeseries('phase_bpsk_2020-04-12T13:39:29.030021.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-04-14T06:08:03.109027.f32', skip = 7),
               load_timeseries('phase_bpsk_2020-04-16T20:29:07.486691.f32', skip = 3),
               load_timeseries('phase_bpsk_2020-04-17T12:08:22.150961.f32', skip = 0),
               load_timeseries('phase_bpsk_2020-04-19T19:47:51.236014.f32', skip = 2),
               load_timeseries('phase_bpsk_2020-04-20T06:22:33.564485.f32', skip = 11),
               load_timeseries('phase_bpsk_2020-04-21T15:55:36.667667.f32', skip = 5),
               load_timeseries('phase_bpsk_2020-04-21T21:07:41.262343.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-04-25T13:02:24.406733.f32', skip = 2),
               load_timeseries('phase_bpsk_2020-04-26T07:15:29.588783.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-04-26T09:55:35.905448.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-04-28T18:09:23.773040.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-04-30T22:26:48.936047.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-05-01T08:22:33.557412.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-05-01T15:08:57.363091.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-05-01T16:36:26.745085.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-05-09T15:20:54.830005.f32', skip = 5),
               load_timeseries('phase_bpsk_2020-05-12T15:33:52.448629.f32', skip = 3),
               load_timeseries('phase_bpsk_2020-05-18T20:28:46.633228.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-05-23T17:58:43.584807.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-05-24T08:48:36.498129.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-05-30T17:22:27.979699.f32', skip = 1),
               load_timeseries('phase_bpsk_2020-06-01T06:16:28.284323.f32', skip = 1),
              ), 'time')

In [5]:
plt.figure(figsize = (12,6), facecolor = 'w')

f.plot()

plt.title('QO-100 BPSK beacon offset from 10489.750 MHz (RX: EA4GPZ)')
plt.ylabel('Offset (Hz)')
plt.xlabel('UTC time')
plt.xticks(rotation = 45)
#plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))
plt.grid()

# eclipse calculator https://www.satellite-calculations.com/Satellite/satellite_eclipse.htm

#plt.axvspan(np.datetime64('2020-04-09T22:04:35'), np.datetime64('2020-04-09T22:31:40'),
#            color = 'grey', alpha = 0.3, label = 'Eclipse')
#plt.axvspan(np.datetime64('2020-04-10T22:08:33'), np.datetime64('2020-04-10T22:27:10'),
#            color = 'grey', alpha = 0.3)
#plt.legend();



In [5]:
plt.figure(figsize = (12,6), facecolor = 'w')

offset_by_day = 0

for j,d in enumerate(list(f.resample(time = '1d'))[-4:]):
    time_of_day = d[1].time - d[1].time[0].astype('datetime64[D]') + np.datetime64('2000-01-01')
    plt.plot(time_of_day, d[1] - offset_by_day*j, label = d[0].astype('datetime64[D]'))

plt.title('QO-100 BPSK beacon offset from 10489.750 MHz (RX: EA4GPZ)')
plt.ylabel('Offset (Hz)')
plt.xlabel('UTC time of day')
plt.xticks(rotation = 45)
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
plt.gca().xaxis.set_major_locator(mdates.HourLocator())
plt.grid()
plt.legend()
# Put a legend to the right of the current axis
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5));



In [6]:
plt.figure(figsize = (12,6), facecolor = 'w')

offset_by_day = 0

f_res = f.resample(time = '1d')
with sns.color_palette('GnBu_d', len(f_res)):
    for j,d in enumerate(f_res):
        time_of_day = d[1].time - d[1].time[0].astype('datetime64[D]') + np.datetime64('2000-01-01')
        plt.plot(time_of_day, d[1] - offset_by_day*j, label = d[0].astype('datetime64[D]'))

plt.title('QO-100 BPSK beacon offset from 10489.750 MHz (RX: EA4GPZ)')
plt.ylabel('Offset (Hz)')
plt.xlabel('UTC time of day')
plt.xticks(rotation = 45)
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
plt.gca().xaxis.set_major_locator(mdates.HourLocator())
plt.grid()
#plt.legend()
# Put a legend to the right of the current axis
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5));



In [7]:
plt.figure(figsize = (12,20), facecolor = 'w')

offset_by_day = 15

with sns.color_palette('GnBu_d', len(f_res)):
    for j,d in enumerate(f_res):
        time_of_day = d[1].time - d[1].time[0].astype('datetime64[D]') + np.datetime64('2000-01-01')
        plt.plot(time_of_day, d[1] - offset_by_day*j, label = d[0].astype('datetime64[D]'))

plt.title('QO-100 BPSK beacon offset from 10489.750 MHz (RX: EA4GPZ)')
plt.ylabel('Offset (Hz)')
plt.xlabel('UTC time of day')
plt.xticks(rotation = 45)
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
plt.gca().xaxis.set_major_locator(mdates.HourLocator())
plt.grid()
plt.legend()
# Put a legend to the right of the current axis
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5));



In [8]:
bpsk = load_timeseries('phase_bpsk_2020-04-12T13:39:29.030021.f32', skip = 1).\
    sel(time = slice('2020-04-12', '2020-04-12T16:45:00'))
ea4gpz_cw = load_timeseries('phase_2020-04-12T13:40:13.670776.f32', freq_tune = 50, skip = 3)

In [9]:
plt.figure(figsize = (12,6), facecolor = 'w')
bpsk.plot()
ea4gpz_cw.plot()

plt.title('NB transponder marker beacons offsets (RX: EA4GPZ)')
plt.legend(['BPSK beacon (10489.750 MHz)', 'EA4GPZ CW (10489.995 MHz)'])
plt.ylabel('Offset (Hz)')
plt.xlabel('UTC time')
plt.xticks(rotation = 45)
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M'))
plt.grid()



In [10]:
plt.figure(figsize = (12,6), facecolor = 'w')
(bpsk-bpsk.mean()).plot()
(ea4gpz_cw-ea4gpz_cw.mean()).plot()

plt.title('NB transponder marker beacons offsets (average removed) (RX: EA4GPZ)')
plt.legend(['BPSK beacon (10489.750 MHz)', 'EA4GPZ CW (10489.995 MHz)'])
plt.ylabel('Offset (average removed) (Hz)')
plt.xlabel('UTC time')
plt.xticks(rotation = 45)
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M'))
plt.grid()