In [1]:

    

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from astropy.time import Time, TimeDelta


    

In [2]:

    

mjd_unixtimestamp_offset = 10587.5
seconds_in_day = 3600 * 24

def mjd2unixtimestamp(m):
    return (m - mjd_unixtimestamp_offset) * seconds_in_day


    

In [3]:

    

def load_file(path):
    ncols = 1+4*3+1
    data = np.fromfile(path, sep=' ')
    return data.reshape((data.size // ncols, ncols))


    

In [4]:

    

data = load_file('/home/daniel/jupyter_notebooks/dslwp/vlbi_report.txt')
t = Time(mjd2unixtimestamp(data[:,0]), format='unix')
ta = data[:,-1]


    

In [5]:

    

def distance(x):
    return np.sqrt(np.sum(x**2, axis=1))

def elevation(x):
    return np.rad2deg(np.arcsin(x[:,2]/distance(x)))


    

In [6]:

    

groundstations = ['Dwingeloo', 'Shahe', 'Harbin', 'Wakayama']
elev_mask = 5
obs_times_start = Time(['2019-06-03 03:05', '2019-06-04 07:00', '2019-06-05 07:00', '2019-06-06 07:00', '2019-06-07 08:00'])
obs_times_end = obs_times_start + TimeDelta(2*3600, format='sec')


    

In [7]:

    

def plot_vlbi_window(start, end):
    sel = (t >= start) & (t <= end)
    fig, ax1 = plt.subplots(figsize = [15,5], facecolor='w')
    ax2 = ax1.twinx()
    min_elev = np.empty(data[sel].shape[0])
    min_elev[:] = np.inf
    for j, gs in enumerate(groundstations):
        x = data[sel, 1 + j*3 : 1 + (j+1)*3]
        elev = elevation(x)
        min_elev[elev < min_elev] = elev[elev < min_elev]
        ax1.plot(t[sel].datetime, elev, label = gs)
    ax2.plot(t[sel].datetime, ta[sel], color = 'purple', linestyle = '--', label = 'True anomaly')
    ax1.set_ylim([0,90])
    ax2.set_ylim([0,360])
    ax1.set_xlim((t[sel][0].datetime, t[sel][-1].datetime))
    for j,obs in enumerate(obs_times_start.datetime):
        if j == 0:
            ax1.axvline(x=obs, color = 'red', label='Observation times (start)', linestyle = '--')
        else:
            ax1.axvline(x=obs, color = 'red', linestyle = '--')
    for j,obs in enumerate(obs_times_end.datetime):
        if j == 0:
            ax1.axvline(x=obs, color = 'red', label='Observation times (end)', linestyle = ':')
        else:
            ax1.axvline(x=obs, color = 'red', linestyle = ':')
    ax1.axhline(y = elev_mask, color = 'grey', linestyle = ':', label = 'Elevation mask')
    ax1.fill_between(t[sel].datetime, 0, 90, where = min_elev >= elev_mask, facecolor = 'green', alpha = 0.1)
    ax1.fill_between(t[sel].datetime, 0, 90, where = min_elev < elev_mask, facecolor = 'red', alpha = 0.1)
    ax1.grid(axis='x')
    ax2.grid(axis='y')
    #plt.legend(loc='upper right')
    ax1.legend()
    ax1.set_ylabel('Elevation (deg)')
    ax2.set_ylabel('True anomaly (deg)', color = 'purple')
    ax1.set_xlabel('UTC time')
    ax1.set_title(f'DSLWP-B elevation and true anomaly {t[sel][0].datetime.strftime("%Y-%m-%d")}')
    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
    ax1.xaxis.set_major_locator(mdates.HourLocator())


    

In [8]:

    

for j in range(int((t[-1] - t[0]).jd)):
    plot_vlbi_window(t[0]+j, t[0]+j+1)