In [1]:
%matplotlib inline

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

c = astropy.constants.c.value

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 = 25
    data = np.fromfile(path, sep=' ')
    return data.reshape((data.size // ncols, ncols))

The file below is obtained by running the dslwp_sme.script GMAT script.


In [4]:
def load_data(path, start, end):
    data = load_file(path)
    t = Time(mjd2unixtimestamp(data[:,0]), format='unix')
    sel = (Time(start) <= t) & (t <= Time(end))
    t = t[sel]
    data = data[sel,:]
    return t, data

In [5]:
t, data = load_data('/home/daniel/jupyter_notebooks/dslwp/SME.txt', '2018-10-07T09:50:00', '2018-10-07T12:10:00')
t2, data2 = load_data('/home/daniel/jupyter_notebooks/dslwp/SME2.txt', '2018-10-19T16:50:00', '2018-10-19T18:50:00')

In [6]:
def rangerate(x):
    return np.sum(x[:,:3] * x[:,3:], axis=1) / np.sqrt(np.sum(x[:,:3]**2, axis=1)) * 1e3

def rangerate2doppler(x, freq = 436.4e6):
    return -x * freq / c

In [7]:
direct_doppler = rangerate2doppler(rangerate(data[:,1:7]))
direct_doppler_435 = rangerate2doppler(rangerate(data[:,1:7]), freq = 435.4e6)
direct_t = t

To compute the Moonbounce Doppler, we first compute the length of the reflection path by laying a fine grid of points on the Moon surface and choosing as reflection point the grid point that minimizes the sum of distances between the point and the satellite and Groundstation. The Doppler is computed as the time derivative of the length of the reflection path.


In [8]:
def compute_moonbounce(data):
    long, lat = np.meshgrid(np.linspace(-np.pi, np.pi, 1000), np.linspace(-np.pi/2, np.pi/2, 1000))
    moon_radius = 1737
    moon_x = moon_radius * np.cos(long) * np.cos(lat)
    moon_y = moon_radius * np.sin(long) * np.cos(lat)
    moon_z = moon_radius * np.sin(lat)
    min_dist = np.empty(data.shape[0])
    reflection_point_long = np.empty(data.shape[0])
    reflection_point_lat = np.empty(data.shape[0])
    for j in range(data.shape[0]):
        dist = np.sqrt((data[j,19]-moon_x)**2 + (data[j,20]-moon_y)**2 + (data[j,21]-moon_z)**2) + np.sqrt((data[j,13]-moon_x)**2 + (data[j,14]-moon_y)**2 + (data[j,15]-moon_z)**2) 
        min_dist[j] = np.min(dist)
        minpoint = np.unravel_index(np.argmin(dist), dist.shape)
        reflection_point_long[j] = long[minpoint]
        reflection_point_lat[j] = lat[minpoint]
    return min_dist, reflection_point_long, reflection_point_lat

In [9]:
min_dist, reflection_point_long, reflection_point_lat = compute_moonbounce(data)
min_dist_2, reflection_point_long_2, reflection_point_lat_2 = compute_moonbounce(data2)

In [10]:
plt.figure(figsize = [15,8], facecolor='w')
plt.plot(np.rad2deg(reflection_point_long), np.rad2deg(reflection_point_lat))
plt.plot(np.rad2deg(reflection_point_long_2), np.rad2deg(reflection_point_lat_2))
plt.imshow(mpimg.imread('/home/daniel/GMAT/R2018a/data/graphics/texture/Moon_HermesCelestiaMotherlode.png'), extent=[-180,180,-90,90])
plt.title('Specular reflection point')
plt.xlabel('Longitude (deg)')
plt.ylabel('Latitude (deg)')
plt.legend(['2018-10-07', '2018-10-19']);



In [11]:
moonbounce_rangerate = np.diff(min_dist)/np.diff(mjd2unixtimestamp(data[:,0])) * 1e3
moonbounce_doppler = rangerate2doppler(moonbounce_rangerate)
moonbounce_doppler_435 = rangerate2doppler(moonbounce_rangerate, 435.4e6)
moonbounce_t = Time(mjd2unixtimestamp((data[:-1,0] + data[1:,0])/2), format='unix')

In [12]:
plt.figure(figsize = [15,8], facecolor='w')
plt.plot(direct_t.datetime, direct_doppler)
plt.plot(moonbounce_t.datetime, moonbounce_doppler)
plt.title('DSLWP-B Doppler at 436.4MHz as seen in Dwingeloo')
plt.xlabel('UTC Time')
plt.ylabel('Doppler (Hz)')
plt.legend(['Direct', 'Moonbounce']);



In [13]:
def waterfall_plot(image_path, start_time, timespan, doppler_offset = 750, band435 = False):
    img = mpimg.imread(image_path)
    img_start = Time(start_time)
    img_end = img_start + TimeDelta(timespan, format='sec')
    plt.figure(figsize = [15,8], facecolor='w')
    plt.plot(direct_t.datetime, direct_doppler if not band435 else direct_doppler_435)
    plt.plot(moonbounce_t.datetime, moonbounce_doppler if not band435 else moonbounce_doppler_435)
    if band435:
        for offset in 1000 + 312.5*np.arange(4):
            plt.plot(direct_t.datetime, direct_doppler + offset, color='C0', linewidth=1, linestyle=':')
            plt.plot(moonbounce_t.datetime, moonbounce_doppler + offset, color='C1', linewidth=1, linestyle=':')
    freq = '436.4MHz' if not band435 else '435.4MHz'
    plt.title(f'DSLWP-B Doppler at {freq} and signal recorded in Dwingeloo (offset {-doppler_offset}Hz)')
    plt.xlabel('UTC Time')
    plt.ylabel('Doppler (Hz)')
    plt.legend(['Direct', 'Moonbounce'])

    min_freq = -1000 if not band435 else 900-1800
    max_freq = 1000 if not band435 else 900+1800
    plt.imshow(img, extent = (mdates.date2num(img_start.datetime), mdates.date2num(img_end.datetime), min_freq-doppler_offset, max_freq-doppler_offset),\
               aspect = 'auto');

In [14]:
waterfall_plot('dslwp_sme0.png', '2018-10-07T09:59:13', 2280.8313)