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

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


In [4]:
data = load_file('/home/daniel/jupyter_notebooks/dslwp/EarthMoonPosDSWLP.txt')
t = Time(mjd2unixtimestamp(data[:,0]), format='unix')

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

def angle(x):
    return np.rad2deg(np.arccos(x[:,0] / distance(x)))

def angular_radius(x, r):
    return np.rad2deg(np.arcsin(r / distance(x)))

def angle_between(x, y):
    return np.rad2deg(np.arccos(np.sum(x * y, axis = 1)/(distance(x)*distance(y))))

In [6]:
earth_angle = angle(data[:,1:4])
moon_angle = angle(data[:,4:7])
earth_radius = 6378.1
moon_radius = 1737.1
earth_angular_radius = angular_radius(data[:,1:4], earth_radius)
moon_angular_radius = angular_radius(data[:,4:7], moon_radius)

In [7]:
camera_fov_long = 18.5
camera_fov_short = 14

In [8]:
photo_times = Time(['2019-06-30 05:51:20', '2019-06-30 05:52:20', '2019-06-30 05:53:20',
                    '2019-06-30 06:29:50', '2019-06-30 06:30:50', '2019-06-30 06:31:50',
                    '2019-07-02 18:56:00', '2019-07-02 18:57:00', '2019-07-02 18:58:00',
                    '2019-07-02 19:31:45', '2019-07-02 19:32:45', '2019-07-02 19:33:45']) + TimeDelta(20, format='sec')
obs_times = Time(['2019-06-30 05:30:00', '2019-07-01 05:30:00', '2019-07-02 18:00:00',
                  '2019-07-03 06:00:00', '2019-07-04 06:30:00', '2019-07-05 07:30:00'])
obs_times_end = obs_times + TimeDelta(2*3600, format='sec')
#calibration_time = Time('2018-10-6 13:55')
#calibration_time = Time('2018-11-8 9:40')

In [9]:
plt.figure(figsize = [15,8], facecolor='w')
plt.plot(t.datetime, earth_angle, color = 'blue', linestyle = '--')
plt.plot(t.datetime, earth_angle-earth_angular_radius, color = 'blue', label='Earth')
plt.plot(t.datetime, earth_angle+earth_angular_radius, color = 'blue')
plt.plot(t.datetime, moon_angle, color = 'grey', linestyle = '--')
plt.plot(t.datetime, moon_angle-moon_angular_radius, color = 'grey', label='Moon')
plt.plot(t.datetime, moon_angle+moon_angular_radius, color = 'grey')
plt.axhline(y=camera_fov_short, linestyle=':', color = 'red')
plt.axhline(y=np.sqrt(camera_fov_short**2 + camera_fov_long**2), linestyle=':', color = 'red', label='Camera FOV (max and min)')
plt.ylim([0,45])
for j,obs in enumerate(obs_times.datetime):
    if j == 0:
        plt.axvline(x=obs, color = 'green', label='Observation times (start)', linestyle = '--')
    else:
        plt.axvline(x=obs, color = 'green', linestyle = '--')
for j,obs in enumerate(obs_times_end.datetime):
    if j == 0:
        plt.axvline(x=obs, color = 'green', label='Observation times (end)', linestyle = ':')
    else:
        plt.axvline(x=obs, color = 'green', linestyle = ':')
#plt.axvline(x=calibration_time.datetime, color = 'orange', label = 'Calibration image')
plt.grid(axis='x')
plt.legend(loc='upper right')
plt.ylabel('Angular distance from camera centre (deg)')
plt.xlabel('UTC time')
plt.title('DSLWP-B camera view');



In [10]:
def plot_view(start = None, end = None, plot_obstimes = False, zoom_factor = 1):
    if start is None or end is None:
        sel = np.arange(t.size)
    else:
        sel = (t >= start) & (t <= end)
    plt.figure(figsize = [15,8], facecolor='w')
    plt.plot(t[sel].datetime, earth_angle[sel], color = 'blue', linestyle = '--')
    plt.plot(t[sel].datetime, earth_angle[sel]-earth_angular_radius[sel], color = 'blue', label='Earth')
    plt.plot(t[sel].datetime, earth_angle[sel]+earth_angular_radius[sel], color = 'blue')
    plt.plot(t[sel].datetime, moon_angle[sel], color = 'grey', linestyle = '--')
    plt.plot(t[sel].datetime, moon_angle[sel]-moon_angular_radius[sel], color = 'grey', label='Moon')
    plt.plot(t[sel].datetime, moon_angle[sel]+moon_angular_radius[sel], color = 'grey')
    plt.axhline(y=camera_fov_short/zoom_factor, linestyle=':', color = 'red')
    plt.axhline(y=np.sqrt(camera_fov_short**2 + camera_fov_long**2)/zoom_factor, linestyle=':', color = 'red', label='Camera FOV (max and min)')
    plt.ylim([0,25/zoom_factor])
    plt.xlim([t[sel][0].datetime, t[sel][-1].datetime])
    if plot_obstimes:
        for j,obs in enumerate(photo_times.datetime):
            if j == 0:
                plt.axvline(x=obs, color = 'green', label='Photo times', linestyle = '--')
            else:
                plt.axvline(x=obs, color = 'green', linestyle = '--')
    plt.grid(axis='x')
    plt.legend(loc='upper right')
    plt.ylabel('Angular distance from camera centre (deg)')
    plt.xlabel('UTC time')
    plt.title('DSLWP-B camera view');

In [11]:
earth_moonrim_angle = angle_between(data[:,1:4], data[:,4:7]) - moon_angular_radius

def plot_earthrise(start = None, end = None, both_edges = False, plot_obstimes = False):
    if start is None or end is None:
        sel = np.arange(t.size)
    else:
        sel = (t >= start) & (t <= end)
    plt.figure(figsize = [15,8], facecolor='w')
    if both_edges:
        plt.plot(t[sel].datetime, earth_moonrim_angle[sel] - earth_angular_radius[sel], color = 'blue', label = 'Earth rim')
        plt.plot(t[sel].datetime, earth_moonrim_angle[sel] + earth_angular_radius[sel], color = 'blue')
    else:
        plt.plot(t[sel].datetime, earth_moonrim_angle[sel], color = 'blue', label = 'Earth centre')
    plt.axhline(y=0, linestyle=':', color = 'grey', label = 'Moon rim')
    plt.ylim([-10,10])
    plt.xlim([t[sel][0].datetime, t[sel][-1].datetime])
    plt.legend()
    if both_edges == True:
        plt.gca().xaxis.set_major_locator(mdates.MinuteLocator(range(0,60,5)))
        plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%d %H:%M'))
        plt.xticks(rotation=45)
    else:
        plt.gca().xaxis.set_major_locator(mdates.DayLocator())
        plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d'))
    if plot_obstimes:
        for j,obs in enumerate(photo_times.datetime):
            if j == 0:
                plt.axvline(x=obs, color = 'green', label='Photo times (start)', linestyle = '--')
            else:
                plt.axvline(x=obs, color = 'green', linestyle = '--')
        #for j,obs in enumerate(photo_times_end.datetime):
        #    if j == 0:
        #        plt.axvline(x=obs, color = 'green', label='Observation times (end)', linestyle = ':')
        #    else:
        #        plt.axvline(x=obs, color = 'green', linestyle = ':')
    plt.grid(axis='x')
    plt.ylabel('Angular distance from Moon rim (deg)')
    plt.xlabel('UTC time')
    plt.title('DSLWP-B Earthrise view')

In [12]:
plot_view('2019-06-30 05:45', '2019-06-30 06:45', plot_obstimes = True, zoom_factor = 2)
plot_view('2019-07-02 18:45', '2019-07-02 19:50', plot_obstimes = True, zoom_factor = 2)



In [13]:
plot_earthrise(start = '2019-07-02 18:45', end = '2019-07-02 19:50', both_edges = True, plot_obstimes = True)



In [33]:
moon_separation = np.empty(photo_times.size)
moon_earth_separation = np.empty(photo_times.size)
moon_speed = np.empty(photo_times.size)
moon_earth_speed = np.empty(photo_times.size)

moon_speed_deg_s = np.diff(moon_angle-moon_angular_radius)/np.diff((t-t[0]).sec)
moon_earth_speed_deg_s = np.diff(earth_moonrim_angle-earth_angular_radius)/np.diff((t-t[0]).sec)
t_moon_speed = t[:-1] + (t[1:] - t[:-1])/2

for j, photo in enumerate(photo_times):
    moon_earth_separation[j] = np.interp((photo-t[0]).jd, (t-t[0]).jd, earth_moonrim_angle-earth_angular_radius)
    moon_separation[j] = np.interp((photo-t[0]).jd, (t-t[0]).jd, moon_angle-moon_angular_radius)
    moon_speed[j] = np.interp((photo-t_moon_speed[0]).jd, (t_moon_speed-t_moon_speed[0]).jd, moon_speed_deg_s)
    moon_earth_speed[j] = np.interp((photo-t_moon_speed[0]).jd, (t_moon_speed-t_moon_speed[0]).jd, moon_earth_speed_deg_s)
    print(photo, 'Moon separation:', moon_separation[j])
    print(photo, 'Moon-Earth separation:', moon_earth_separation[j])


2019-06-30 05:51:40.000 Moon separation: 4.678689588126566
2019-06-30 05:51:40.000 Moon-Earth separation: -11.682121647254913
2019-06-30 05:52:40.000 Moon separation: 2.938698709556035
2019-06-30 05:52:40.000 Moon-Earth separation: -12.478516390848121
2019-06-30 05:53:40.000 Moon separation: 1.1603829345735788
2019-06-30 05:53:40.000 Moon-Earth separation: -13.204947111115477
2019-06-30 06:30:10.000 Moon separation: -0.4663544189253721
2019-06-30 06:30:10.000 Moon-Earth separation: 27.01692824324825
2019-06-30 06:31:10.000 Moon separation: 2.007805767735777
2019-06-30 06:31:10.000 Moon-Earth separation: 29.347256182031156
2019-06-30 06:32:10.000 Moon separation: 4.498122327999905
2019-06-30 06:32:10.000 Moon-Earth separation: 31.698237612270844
2019-07-02 18:56:20.000 Moon separation: 4.814369877387829
2019-07-02 18:56:20.000 Moon-Earth separation: 2.9534104241287276
2019-07-02 18:57:20.000 Moon separation: 2.6903515085191056
2019-07-02 18:57:20.000 Moon-Earth separation: 0.8530852354495602
2019-07-02 18:58:20.000 Moon separation: 0.506036755688968
2019-07-02 18:58:20.000 Moon-Earth separation: -1.3071980169644475
2019-07-02 19:32:05.000 Moon separation: -0.5256301384398874
2019-07-02 19:32:05.000 Moon-Earth separation: -1.4328252183367844
2019-07-02 19:33:05.000 Moon separation: 2.675694390495612
2019-07-02 19:33:05.000 Moon-Earth separation: 1.74577976975264
2019-07-02 19:34:05.000 Moon separation: 5.887683618325296
2019-07-02 19:34:05.000 Moon-Earth separation: 4.935548221765514

In [29]:
photo_separation = np.array([np.nan, 6.469961843009586, 4.700435772416988, np.nan, -1.3292181866248751, 1.3121612143594001] + 6*[np.nan])
photo_moon_earth_separation = np.array(6*[np.nan] + [np.nan, 1.32, np.nan, np.nan, 1.06, 4.36])

In [16]:
angle_difference = photo_separation - moon_separation
angle_difference


Out[16]:
array([        nan,  3.53126313,  3.54005284,         nan, -3.33702395,
       -3.18596111,         nan,         nan,         nan,         nan,
               nan,         nan])

In [17]:
time_difference = angle_difference/moon_speed
time_difference


Out[17]:
array([          nan, -120.44259247, -118.1510677 ,           nan,
        -80.62889492,  -76.54770661,           nan,           nan,
                 nan,           nan,           nan,           nan])

In [18]:
separation_error = 3.25
time_error = np.abs(3.25/moon_speed)
time_error


Out[18]:
array([113.29931411, 110.84940735, 108.47040641,  79.10985236,
        78.52622938,  78.08634117,  93.10260834,  90.52601615,
        88.03569411,  61.07355722,  60.77454005,  60.66756668])

In [32]:
angle_difference_moon_earth = photo_moon_earth_separation - moon_earth_separation
angle_difference_moon_earth


Out[32]:
array([        nan,         nan,         nan,         nan,         nan,
               nan,         nan,  0.46691476,         nan,         nan,
       -0.68577977, -0.57554822])

In [35]:
angle_difference_moon_earth/moon_earth_speed


Out[35]:
array([         nan,          nan,          nan,          nan,
                nan,          nan,          nan, -13.15119332,
                nan,          nan, -12.9144165 , -10.81768429])