In [265]:
from IPython.display import Image, display
In [2]:
%matplotlib inline
%config InlineBackend.figure_format='retina'
from __future__ import absolute_import, division, print_function
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib.patches import Rectangle
from matplotlib.pyplot import GridSpec
from itertools import combinations, islice, takewhile
import seaborn as sns
import mpld3
import numpy as np
import pandas as pd
import scipy.sparse
import os, sys
import warnings
from astropy import units as u
from astropy.coordinates import SkyCoord
import statsmodels.api as sm
import warnings
warnings.filterwarnings('ignore')
sns.set_context('poster', font_scale=1.3)
from ipywidgets import interact, FloatSlider, SelectMultiple, interactive
In [6]:
# 'constants'
from scipy.constants import c # speed of light in m/s
DIP_RA = 17.3
DIP_RA_ERR = 1.0
DIP_DEC = -61.0
DIP_DEC_ERR = 10.0
DIP_AMPLITUDE = 0.97e-5
DIP_AMPLITUDE_ERR = 0.21e-5 # average of asymmetric errors
DIP_MONOPOLE = -0.178e-5
DIP_MONOPOLE_ERR = 0.084e-5
dipole = SkyCoord(DIP_RA, DIP_DEC, unit=(u.hourangle, u.deg))
# Data
vltqs = pd.read_csv('../data/vlt-transitions-new.tsv', sep='\t')
vltqs['trans'] = vltqs['transition'] + vltqs.four.astype(str)
keckqs = pd.read_csv("../data/keck-transitions-new.tsv", sep='\t')
keckqs['trans'] = keckqs.transition + keckqs.four.astype(str)
qvals = pd.read_csv("../data/qvalues.txt", sep=' ', index_col=0)
codes = pd.concat([keckqs.set_index('code')[['trans', 'wavelength']],
vltqs.set_index('code')[['trans', 'wavelength']]])
full_parse = pd.read_csv("../data/full-parse.tsv", sep='\t')
In [257]:
def observed_shifts(telescope='VLT'):
waves = []
shifts = []
for index, row in full_parse[full_parse.source.eq(telescope)].iterrows():
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
shifts.append(shifted_velocity(row.da, qval, rest_wave))
return np.array(waves), np.array(shifts)
def parse_j2000(name):
"""Takes the J2000 name stored in the results and returns it in a format astropy can understand."""
return ' '.join([name[1:3], name[3:5], name[5:7], name[7:10], name[10:12], name[12:]])
def j2000_to_theta(name):
"""Returns the angle (degrees) between the position on the sky from
a given `name` and the position of the dipole model from 2012, King."""
c = SkyCoord(parse_j2000(name), unit=(u.hourangle, u.deg))
return float(c.separation(dipole).to_string(decimal=True))
def dipole_alpha(name):
"""Returns the value of Delta alpha/alpha as given by the best fit 2012 King model for
the given name (position).
"""
theta = j2000_to_theta(name)
return (DIP_AMPLITUDE * np.cos(np.deg2rad(theta)) + DIP_MONOPOLE) * 1e6
def plot_system(system, dataframe):
plotdf = dataframe[dataframe.system == system]
fig, (ax1, ax2) = plt.subplots(figsize=(12, 10), nrows=2)
sns.regplot('wavelength', 'vshift', data=plotdf, ax=ax1)
sns.regplot('x', 'vshift', data=plotdf, ax=ax2)
fig.tight_layout()
def generate_dataset(gen_dipole_alpha=True,
wavelength_distortion=False,
seed=228,
):
df = pd.DataFrame(columns=['system',
'#J2000',
'source',
'wavelength',
'vshift',
'sigma',
'qval',
'rest_wave'
])
count = 0
abs_count = 0
for index, row in full_parse.iterrows():
waves = []
rest_waves = []
vshifts = []
qvals_list = []
for tran in row['transition'].split():
vshift = 0.0
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
if gen_dipole_alpha:
vshift += shifted_velocity(dipole_alpha(row['#J2000']),
qval,
rest_wave)
if wavelength_distortion:
vshift += distorted_velocity(row, measured_wave)
waves.append(measured_wave)
rest_waves.append(rest_wave)
vshifts.append(vshift)
qvals_list.append(qval)
vshifts = np.array(vshifts)
errors = sigmas(vshifts)
vshifts += errors * np.random.randn(len(vshifts))
vshifts = vshifts - vshifts[0]
for single in range(len(vshifts)):
abs_count += 1
df.loc[abs_count] = [int(index), #system
row['#J2000'], #j2000
row.source,
waves[single], #wavelength
vshifts[single],
errors[single],
qvals_list[single], # qvalues
rest_waves[single],
]
df['system'] = df.system.astype(int)
# For numerical stability during fitting, dividing the x-values by 1.3e13
df['x'] = -2.0 * c * df['qval'] * df['rest_wave'] / 1.0e13
return df
def shifted_velocity(del_alpha, q, lamb):
# vj =v0 + ∆α xj, xj =−2cqjλ0j,
x = -2 * c * q * lamb
return del_alpha * x * 1e-14
def sigmas(vshifts):
return np.random.rand(len(vshifts)) * 30.0 + 50.0
def VLT_distortion(measured_wave,
cutoff=4500.,
slope1=0.06,
intercept1 = -100.0,
slope2 =0.160,
intercept2=-1500.0,
):
"""Telescope dependent distortion function for the VLT sample."""
if measured_wave < cutoff:
return measured_wave * slope1 + intercept1
else:
return measured_wave * slope2 + intercept2
def Keck_distortion(measured_wave, cutoff=10000.):
"""Telescope dependent distortion function for the Keck sample."""
slope1 = .0600
intercept1 = -100
slope2 = .160
intercept2 = -1500
if measured_wave < cutoff:
return measured_wave * slope1 + intercept1
else:
return measured_wave * slope2 + intercept2
def distorted_velocity(row, measured_wave):
"""Telescope dependent distortion function for the VLT sample."""
if row.source == "VLT":
return VLT_distortion(measured_wave)
elif row.source == "Keck":
return Keck_distortion(measured_wave)
In [266]:
# Todo: make w/x process be rows instead of columns.
In [285]:
def fit_hypothesis(system=0, dataframe1=df_x, hypothesis='x'):
"""Return the chisq and the fit model object for a given dataframe and hypothesis."""
plotdf1 = dataframe1[dataframe1.system == system]
assert(hypothesis in ['x', 'w'])
if hypothesis == 'x':
X = sm.add_constant(plotdf1.x)
elif hypothesis == 'w':
X = sm.add_constant(plotdf1.wavelength)
results = sm.WLS(plotdf1.vshift, X, weights=1.0/plotdf1.sigma).fit()
chisq = np.sum((plotdf1.vshift - results.fittedvalues)**2.0 / (plotdf1.sigma) ** 2.0)
return chisq, results
def plot_hypotheses(system, dataframe1, dataframe2):
plotdf1 = dataframe1[dataframe1.system == system]
plotdf2 = dataframe2[dataframe2.system == system]
fig, ((ax2, ax4), (ax3, ax1)) = plt.subplots(figsize=(14, 10), nrows=2, ncols=2, )
chi_one, mod_one = fit_hypothesis(system=system, dataframe1=dataframe1, hypothesis='w')
ax1.errorbar(plotdf1.wavelength, plotdf1.vshift, yerr=plotdf1.sigma, ls='none',
label=r'$\chi^2$ = ' + str(np.round(chi_one, 2)), color=sns.color_palette()[2])
ax1.scatter(plotdf1.wavelength, plotdf1.vshift, color=sns.color_palette()[2], label='')
ax1.plot(plotdf1.wavelength, mod_one.fittedvalues, color='k')
chi_one, mod_one = fit_hypothesis(system=system, dataframe1=dataframe1, hypothesis='x')
ax3.errorbar(plotdf1.x, plotdf1.vshift, yerr=plotdf1.sigma, ls='none',
label=r'$\chi^2$ = ' + str(np.round(chi_one, 2)), color=sns.color_palette()[0])
ax3.scatter(plotdf1.x, plotdf1.vshift, color=sns.color_palette()[0], label='')
ax3.plot(plotdf1.x, mod_one.fittedvalues, color='k')
chi_one, mod_one = fit_hypothesis(system=system, dataframe1=dataframe2, hypothesis='w')
ax2.errorbar(plotdf2.wavelength, plotdf2.vshift, yerr=plotdf2.sigma, ls='none',
label=r'$\chi^2$ = ' + str(np.round(chi_one, 2)), color=sns.color_palette()[0])
ax2.scatter(plotdf2.wavelength, plotdf2.vshift, color=sns.color_palette()[0], label='')
ax2.plot(plotdf2.wavelength, mod_one.fittedvalues, color='k')
chi_one, mod_one = fit_hypothesis(system=system, dataframe1=dataframe2, hypothesis='x')
ax4.errorbar(plotdf2.x, plotdf2.vshift, yerr=plotdf2.sigma, ls='none',
label=r'$\chi^2$ = ' + str(np.round(chi_one, 2)), color=sns.color_palette()[2])
ax4.scatter(plotdf2.x, plotdf2.vshift, color=sns.color_palette()[2], label='')
ax4.plot(plotdf2.x, mod_one.fittedvalues, color='k')
autoAxis = ax2.axis()
rec = Rectangle((autoAxis[0],
autoAxis[2]),
(autoAxis[1]-autoAxis[0]),
(autoAxis[3]-autoAxis[2]),
fill=False,lw=2)
rec = ax2.add_patch(rec)
rec.set_clip_on(False)
autoAxis = ax3.axis()
rec = Rectangle((autoAxis[0],
autoAxis[2]),
(autoAxis[1]-autoAxis[0]),
(autoAxis[3]-autoAxis[2]),
fill=False,lw=2)
rec = ax3.add_patch(rec)
rec.set_clip_on(False)
ax3.set_title("X process X fit")
ax2.set_title("W process W fit")
ax1.set_title("X process W fit")
ax4.set_title("W process X fit")
ax2.set_ylabel("vshift [m/s]")
ax3.set_ylabel("vshift [m/s]")
ax1.set_xlabel("wavelength")
ax2.set_xlabel("wavelength")
ax3.set_xlabel("x")
ax4.set_xlabel("x")
leg = ax1.legend(handlelength=0, handletextpad=0, fancybox=True, frameon=True, facecolor='white', loc='best')
for item in leg.legendHandles:
item.set_visible(False)
leg = ax2.legend(handlelength=0, handletextpad=0, fancybox=True, frameon=True, facecolor='white', loc='best')
for item in leg.legendHandles:
item.set_visible(False)
leg = ax3.legend(handlelength=0, handletextpad=0, fancybox=True, frameon=True, facecolor='white', loc='best')
for item in leg.legendHandles:
item.set_visible(False)
leg = ax4.legend(handlelength=0, handletextpad=0, fancybox=True, frameon=True, facecolor='white', loc='best')
for item in leg.legendHandles:
item.set_visible(False)
fig.tight_layout()
In [286]:
plot_hypotheses(231, df_x, df_w)
In [188]:
df_x = generate_dataset(gen_dipole_alpha=True, wavelength_distortion=False)
df_w = generate_dataset(gen_dipole_alpha=False, wavelength_distortion=True)
In [259]:
df_x2 = generate_dataset(gen_dipole_alpha=True, wavelength_distortion=False)
df_w2 = generate_dataset(gen_dipole_alpha=False, wavelength_distortion=True)
In [291]:
plot_hypotheses(169, df_x, df_w)
In [327]:
int(False)
Out[327]:
In [322]:
False & False
Out[322]:
In [339]:
def add_binary(num1, num2):
st1 = str(num1)
st2 = str(num2)
total = []
carry = False
for index in range(1, max([len(st1), len(st2)]) + 2):
if index <= len(st1):
first = bool(int(st1[-index]))
else:
first = False
if index <= len(st2):
second = bool(int(st2[-index]))
else:
second = False
print(first, second, carry, first | second | carry)
total.append(first | second | carry)
if first & second:
carry = True
# total.insert(0, first | second | carry)
if ~(first | second):
print(index, "not first and not second")
carry = False
return ''.join([str(int(x)) for x in total[::-1]])
add_binary(1, 1)
Out[339]:
In [ ]:
In [342]:
(True | True) & False
Out[342]:
In [ ]:
In [ ]:
sns.set_context(font_scale=)
In [344]:
plot_hypotheses(264, df_x, df_w)
In [198]:
df_x[df_x.system==264][['wavelength', 'x', 'vshift', 'sigma']]
Out[198]:
In [230]:
chisqs = {# data_hypo
'x_x':[],
'x_w':[],
'w_x':[],
'w_w':[]}
data = {'x':df_x, 'w':df_w}
for system in sorted(df_w.system.unique()):
for df in ['x', 'w']:
for hypo in ['x', 'w']:
chisq, results = fit_hypothesis(system=system, dataframe1=data[df], hypothesis=hypo)
chisqs['_'.join([df, hypo])].append(chisq)
for item in chisqs:
chisqs[item] = np.array(chisqs[item])
In [245]:
np.sum(chisqs['w_x']/chisqs['w_w'] < 1.0) / len(chisqs['w_x'])
Out[245]:
In [253]:
np.sum(chisqs['x_x']/chisqs['x_w'] < 1.0) / len(chisqs['w_x'])
Out[253]:
In [247]:
full_parse.head()
Out[247]:
In [251]:
angles = []
dipole_alphas = []
measured_alphas = []
for index, row in full_parse.iterrows():
name = row['#J2000']
angles.append(j2000_to_theta(name))
dipole_alphas.append(dipole_alpha(name))
measured_alphas.append(row.da)
In [264]:
plt.scatter(angles, dipole_alphas)
plt.scatter(angles, measured_alphas, alpha=0.5)
Out[264]:
In [263]:
plt.scatter(angles, measured_alphas)
Out[263]:
In [ ]:
# Todo: ipywidgets minimization of spline + offsets for vshift
In [ ]:
In [329]:
alpha = []
distortion = []
for (index, g) in df2.groupby('system'):
X = sm.add_constant(g.x)
model = sm.WLS(g.vshift, X, weights=1.0/g.sigma)
results = model.fit()
chisq = np.sum((g.vshift - results.fittedvalues)**2.0 / (g.sigma) ** 2.0)
if chisq < 1e7:
alpha.append(chisq)
X = sm.add_constant(g.wavelength)
model = sm.WLS(g.vshift, X, weights=1.0/g.sigma)
results = model.fit()
chisq = np.sum((g.vshift - results.fittedvalues)**2.0 / (g.sigma) ** 2.0)
if chisq < 1e10:
distortion.append(chisq)
alpha = np.array(alpha)
distortion = np.array(distortion)
fig, ax = plt.subplots(figsize=(12, 8))
ax.hist(alpha, alpha=0.5,
bins=np.linspace(0, 5, 50),
label='Alpha hypothesis')
ax.hist(distortion, alpha=0.5,
bins=np.linspace(0, 5, 50),
label='Distortion hypothesis')
# ax.hist(alpha / distortion,
# bins=np.linspace(0, 5, 50),
# label='Wavelength',
# alpha=0.5)
# ax.set_xlim(0, 5.0, )
ax.legend()
fig.tight_layout()
In [ ]:
We have many measurements of $\alpha$ on both the Keck and VLT telescopes. Measurements with spectrographs both telescopes reveal statistically significant (though small in absolute terms) differences between the expected positions of absorption lines and the measured positions.
There are many possible hypotheses that could explain these differences. To list just two:
We can consider any velocity shift measurements as a hypothesis test. Previous studies started with the assumption that any shift in absorption lines would come as a result of a change in $\alpha$. We have the reported $\frac{\Delta \alpha}{\alpha}$ values for both the Keck and VLT samples. If we invert the hypothesis: use the measured alpha values as a system of velocity shifts, we can effectively discriminate which hypothesis best fits the data.
The goal will be to simulate what the velocity shifts look like if the dipole model is correct, and compare that to what we currently see.
The end goal of this analysis is ultimately a change in how $\frac{\Delta \alpha}{\alpha}$ is measured. The proposed recipe is:
vpfit already allows for this).This will likely have rather low power for any single absorption system, but for ensembles, I think that it's likely the only case
In [5]:
species = set([spec[0] for spec in qvals.index.str.split('I')])
colors = sns.color_palette(n_colors=len(species))
plot_colors = {}
for index, specie in enumerate(sorted(species)):
plot_colors[specie] = colors[index]
In [11]:
def plot_absorption(specie=['Al', 'Mg'],
z=1.3,
daa=0.05,
):
fig, ax = plt.subplots(figsize=(12, 6))
ax.hlines(1, 0, 10000)
for tran, row in qvals[qvals.index.str.startswith(tuple(specie))].iterrows():
ax.vlines((1.0 + z) * row.wave + row.qval * ((daa + 1)**2.0-1.0),
0, 1, color=plot_colors[tran[:2]])
ax.vlines((1.0 + z) * row.wave,
0, 1, color=plot_colors[tran[:2]], alpha=0.2)
ax.set_xlim(3000, 8e3)
ax.set_ylim(0, 1.1)
ax.set_ylabel("Normalized Flux")
ax.set_xlabel(r"Wavelength $\AA$")
for spec in specie:
ax.vlines(-1, -1, 0, color=plot_colors[spec], label=spec)
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
fig.tight_layout()
plot_absorption(specie=['Al', 'Mg', 'Fe', 'Ni'], )
In [51]:
def plot_shift(specie=['Al', 'Mg'],
z=1.3,
daa=0.05,
):
fig, ax = plt.subplots(figsize=(12, 6))
ax.hlines(0, 0, 10000)
for tran, row in qvals[qvals.index.str.startswith(tuple(specie))].iterrows():
ax.scatter((1.0 + z) * row.wave,
# row.qval * ((daa + 1)**2.0-1.0),
shifted_velocity(daa, row.qval, (1.0 + z) * row.wave),
color=plot_colors[tran[:2]], zorder=3)
ax.set_xlim(3000, 8e3)
ax.set_ylim(-200, 200)
ax.set_ylabel("Velocity Shift [m/s]")
ax.set_xlabel(r"Observed Wavelength $\AA$")
for spec in specie:
ax.vlines(-1, -1, 0, color=plot_colors[spec], label=spec, zorder=-3)
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
fig.tight_layout()
plot_shift(specie=['Al', 'Mg', 'Fe', 'Ni'], )
In [52]:
redshift_slider = FloatSlider(value=2.1, min=0.0, max=3.5)
alpha_slider = FloatSlider(value=0.0, min=-5.0, max=5.0, step=0.1)
species_select = SelectMultiple(options=['Al', 'Cr', 'Fe', 'Mg', 'Mn', 'Ni', 'Si'],
description="Species",
value=['Al', 'Fe']
)
# w = interactive(plot_absorption,
# specie=species_select,
# z=redshift_slider,
# daa=alpha_slider,
# )
shift_interactive = interactive(plot_shift,
specie=species_select,
z=redshift_slider,
daa=alpha_slider,
)
In [53]:
shift_interactive
In [42]:
full_parse.head()
Out[42]:
In [15]:
def plot_shifts(telescope='VLT'):
fig, ax = plt.subplots(figsize=(12, 6))
w, s = observed_shifts(telescope=telescope)
sns.regplot(w, s, lowess=True, scatter_kws={'alpha':0.2}, ax=ax)
ax.set_ylim(-200, 200)
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title(telescope)
fig.tight_layout()
In [16]:
plot_shifts(telescope='VLT')
In [17]:
plot_shifts(telescope='Keck')
In [59]:
def plot_all(title='Keck + VLT'):
fig, ax = plt.subplots(figsize=(12, 8))
telescope='VLT'
w, s = observed_shifts(telescope=telescope)
sns.regplot(w, s, lowess=True, scatter_kws={'alpha':0.5}, ax=ax, label=telescope)
telescope='Keck'
w, s = observed_shifts(telescope=telescope)
sns.regplot(w, s, lowess=True, scatter_kws={'alpha':0.5}, ax=ax, label=telescope)
ax.set_xlim(2000, 10000)
ax.set_ylim(-200, 200)
ax.legend()
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title(title)
fig.tight_layout()
In [60]:
plot_all()
In [61]:
def observed_shifts_by_sample(sample='A'):
waves = []
shifts = []
for index, row in full_parse[full_parse['sample'] == sample].iterrows():
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
shifts.append(shifted_velocity(row.da, qval, rest_wave))
return np.array(waves), np.array(shifts)
def plot_samples():
fig, ax = plt.subplots(figsize=(12, 8))
telescope='VLT'
for sample in sorted(full_parse['sample'].unique()):
w, s = observed_shifts_by_sample(sample=sample)
sns.regplot(w, s, lowess=True, scatter_kws={'alpha':0.5}, ax=ax, label=sample)
ax.set_ylim(-400, 400)
ax.set_xlim(2000, 10000)
ax.legend()
# frame = ax.get_frame()
# frame.set_alpha(1.0)
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title('Samples')
fig.tight_layout()
In [62]:
plot_samples()
In [54]:
w
In [74]:
def plot_sim(telescope='VLT', use_dipole=True):
fig, ax = plt.subplots(figsize=(12, 6))
w, s = observed_shifts(telescope=telescope)
sns.regplot(w, s, lowess=True, scatter_kws={'alpha':0.05}, ax=ax, label='Implied by Measured VLT')
w, s = simulated_shifts(telescope=telescope, use_dipole=use_dipole)
sns.regplot(w, s, lowess=True, scatter_kws={'alpha':0.05}, ax=ax, label='Implied by Dipole VLT')
ax.set_ylim(-200, 200)
ax.legend()
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title('Simulated velocity shifts for ' + telescope)
fig.tight_layout()
In [75]:
plot_sim()
In [ ]:
In [ ]:
In [69]:
def simulated_shifts(telescope='VLT', use_dipole=True):
waves = []
shifts = []
for index, row in full_parse[full_parse.source.eq(telescope)].iterrows():
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
# da = row.da
if use_dipole:
da = dipole_alpha(row['#J2000'])
else:
da = -3.0
shifts.append(shifted_velocity(da, qval, rest_wave))
return np.array(waves), np.array(shifts)
plot_sim(use_dipole=False)
In [ ]:
In [ ]:
In [80]:
Out[80]:
In [99]:
telescope = 'VLT'
row = full_parse[full_parse.source.eq(telescope)].sample(1,
random_state=2
).iloc[0]
waves = []
shifts = []
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
# da = row.da
# if use_dipole:
da = dipole_alpha(row['#J2000'])
# else:
# da = -3.0
shifts.append(shifted_velocity(da, qval, rest_wave))
wav, shift = np.array(waves), np.array(shifts)
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(wav, shift)
ax.set_ylim(-200, 200)
# ax.legend()
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title('Simulated velocity shifts for ' + telescope)
fig.tight_layout()
In [ ]:
telescope = 'VLT'
row = full_parse[full_parse.source.eq(telescope)].sample(1,
random_state=2
).iloc[0]
waves = []
shifts = []
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
# da = row.da
# if use_dipole:
da = dipole_alpha(row['#J2000'])
# else:
# da = -3.0
shifts.append(shifted_velocity(da, qval, rest_wave))
wav, shift = np.array(waves), np.array(shifts)
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(wav, shift)
ax.set_ylim(-200, 200)
# ax.legend()
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title('Simulated velocity shifts for ' + telescope + " " + str(np.round(da, 2)))
fig.tight_layout()
In [113]:
telescope = 'VLT'
row = full_parse[full_parse.source.eq(telescope)].sample(1,
random_state=3
).iloc[0]
waves = []
shifts = []
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
# da = row.da
# if use_dipole:
da = dipole_alpha(row['#J2000'])
# else:
# da = -3.0
shifts.append(shifted_velocity(da, qval, rest_wave))
wav, shift = np.array(waves), np.array(shifts)
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(wav, shift)
ax.set_ylim(-200, 200)
# ax.legend()
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title('Simulated velocity shifts for ' + telescope + " " + str(da))
fig.tight_layout()
In [112]:
telescope = 'VLT'
row = full_parse[full_parse.source.eq(telescope)].sample(1,
random_state=8
).iloc[0]
waves = []
shifts = []
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
da = dipole_alpha(row['#J2000'])
shifts.append(shifted_velocity(da, qval, rest_wave))
wav, shift = np.array(waves), np.array(shifts)
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(wav, shift)
waves = []
shifts = []
for tran in row['transition'].split():
rest_wave = qvals.loc[codes.loc[tran].trans].wave
measured_wave = rest_wave * (1 + row.zabs)
qval = qvals.loc[codes.loc[tran].trans].qval
waves.append(measured_wave)
da = row.da
shifts.append(shifted_velocity(da, qval, rest_wave))
wav, shift = np.array(waves), np.array(shifts)
ax.scatter(wav, shift)
ax.set_xlabel(r"Observed Wavelength [$\AA$]")
ax.set_ylabel("Velocity [m/s]")
ax.set_title('Simulated velocity shifts for ' + telescope + " " + str(da))
fig.tight_layout()
In [ ]:
In [ ]:
In [ ]: