In [ ]:

    

from os import path
import pickle

import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('apw-notebook')
%matplotlib inline
from scipy.optimize import minimize
from comoving_rv.plot import colors


    

In [ ]:

    

night_path = '/Volumes/ProjectData/gaia-comoving-followup/processed/mdm-spring-2017/n1/'
wavelength_gp_path = path.join(night_path, 'wavelength_GP_model.pickle')

with open(wavelength_gp_path, 'rb') as f:
    model = pickle.load(f)


    

In [ ]:

    

x = model.x
y = model.y
med_x = model.x_shift


    

In [ ]:

    

def func(p):
    return (6562.8 - model.gp.predict(model.y, p[0]-model.x_shift, return_var=False, return_cov=False)[0])**2

res = minimize(func, x0=687.)
Ha_pix = res.x[0]


    

In [ ]:

    

x_grid = np.linspace(0, 1600, 1024) - med_x


    

In [ ]:

    

mu, var = model.gp.predict(y, x_grid, return_var=True)
std = np.sqrt(var)

y_mu, var = model.gp.predict(y, x, return_var=True)

# Plot the maximum likelihood model
fig,ax = plt.subplots(1, 1, figsize=(8,6))

# data
ax.scatter(x + med_x, y - y_mu, marker='o', zorder=100)

ax.plot(x_grid+med_x, mu-mu, color=colors['gp_model'], marker='')
ax.fill_between(x_grid+med_x, std, -std, color=colors['gp_model'],
                alpha=0.3, edgecolor="none", linewidth=0)

ax.set_xlabel('pixel')
ax.set_ylabel(r'wavelength residual [${\rm \AA}$]')

ax.set_xlim(x_grid.min()+med_x, x_grid.max()+med_x)
ax.set_ylim(-1, 1)

# right axis, velocity
ax_r = ax.twinx()
ax_r.set_ylim([x/6563*3e5 for x in ax.get_ylim()])
ax_r.set_ylabel(r'velocity error at ${{\rm H}}\alpha$ [{}]'
                .format((u.km/u.s).to_string(format='latex_inline')))

# top axis, wavelength
ax_t = ax.twiny()
ax_t.set_xlim(model.gp.predict(model.y, np.array([0,1600])-model.x_shift, return_var=False, return_cov=False))
ax_t.set_xlabel(r'(mean) wavelength [${\rm \AA}$]', labelpad=12)

ax.axvline(Ha_pix, color=colors['line_marker'], alpha=0.5, zorder=-10)
ax.text(Ha_pix, 0.85, r'${\rm H}\alpha$', bbox=dict(facecolor='w'),
        fontsize=18, color=colors['line_marker'], ha='center')

# shaded region
ax.axvspan(0, 250, color='#666666', alpha=0.25, zorder=100, linewidth=0)
ax.axvspan(1210, 1600, color='#666666', alpha=0.25, zorder=100, linewidth=0)

fig.tight_layout()

# fig.savefig('wavelength_gp.pdf')


    

In [ ]:

    




    

In [ ]: