In [1]:
import glob
import sys
from copy import deepcopy
import numpy as np
import scipy.constants
import scipy.stats
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import astropy
from astropy.table import Table, Column
import astropy.units as u
from astropy.io import fits
from astropy.modeling import models, fitting
import astropy.time
import astropy.io.ascii
import spectrum
from spectrum.spectrum import Spectrum
sys.path.append('/data/hguenther/Dropbox/code/python/')
import TWHya
%matplotlib inline
In [15]:
datapath = '/data/hguenther/obs/HST/TW_Hya/data/'
posterdir = '/data/hguenther/Dropbox/my_poster/14_CS18/'
In [16]:
colorlist = [(1,0,0), (1,.2,0), (1., .5,0), 'b', 'g','c','m', 'grey', 'chartreuse', 'k',
'sienna', 'lime', 'lightsalmon', '0.2', '0.4', '0.6', '0.8', (.2,1.,.2), 'r','g']
def plot_spectra(*args):
#plt.clf()
for arg in args:
for tab, c in zip(arg, colorlist):
line = plt.plot(tab['WAVELENGTH'].T, tab['FLUX'].T*1e12, color = c)
line[0].set_label('{0:s}'.format(tab.meta['DATE-OBS'])) # label only the first line same color
plt.xlabel(r'$\lambda [\AA]$')
plt.ylabel(r'flux $\left[10^{-12} \frac{\rm{erg}}{\rm{s\;cm}^2 \AA}\right]$')
def plot_spectra_diff(*args):
plt.clf()
for arg in args:
for tab, c in zip(arg[1:], colorlist):
f = interp1d(fuv[0]['WAVELENGTH'], fuv[i].FLUX[1], 'nearest', bounds_error = False)
w = tab['WAVELENGTH'].T
flux = tab['FLUX'].T
line = plt.plot(w, flux-f(w), color = c)
line[0].set_label(tab.meta['DATE-OBS']) # label only the first line same color
plt.xlabel(r'$\lambda [\AA]$')
plt.ylabel(r'flux [erg/s/cm$^2/\AA$]')
In [17]:
filelist = glob.glob(datapath + '*sum.fits')
fuv = []
nuv = []
for f in filelist:
tab = Spectrum.read(f)
if max(tab.disp) > 2000 * u.Angstrom:
nuv.append(tab)
else:
fuv.append(tab)
In [19]:
#Merge FUV and NUV of each observation into one array for easier plotting
nf = []
for f, n in zip(fuv, nuv):
disp = np.hstack([f.disp.value, n.disp.value]) * f.disp.unit
nf.append(spectrum.spectrum.coadd_simple([f,n], dispersion=disp, bounds_error=False))
nf[-1].meta = n.meta
In [20]:
plt.figure(figsize=(20,4))
plot_spectra(nf)
plt.ylim([0,1.5])
plt.xlim([1380, 2830])
plt.legend(loc='center')
fig = plt.gcf()
ax = plt.gca()
ax.set_axis_bgcolor('0.95')
ax.text(2500,1.0, 'COS NUV channel', fontsize='xx-large')
ax.text(1400,1.0, 'COS FUV channel', fontsize='xx-large')
plt.setp(ax.get_xticklabels(), fontsize='xx-large')
ax.set_xlabel(ax.get_xlabel(), fontsize='xx-large')
plt.setp(ax.get_yticklabels(), fontsize='xx-large')
ax.set_ylabel(ax.get_ylabel(), fontsize='xx-large')
#ann = ax.annotate('variable\ncontinuum', xy=(1700,.1), xycoords='data',
# xytext=(1800, 0.6), textcoords='data',
# size=20, va="center",
# bbox=dict(boxstyle="round", fc=(1.0, 0.7, 0.7), ec="none"),
# arrowprops=dict(arrowstyle="wedge,tail_width=1.",
# fc=(1.0, 0.7, 0.7), ec="none",
# patchA=None,
# patchB=None,
# relpos=(0.2, 0.5),
# )
# )
plt.title('TW Hya observed with HST/COS for 10 orbits', fontsize='xx-large')
fig.subplots_adjust(left=0.05, right=0.99, bottom=0.15)
plt.savefig(posterdir + 'longspec.png', dpi=600, max_dpi=None, facecolor='0.8', edgecolor='none')
In [184]:
fig = plt.figure(figsize=(12,5))
fig = plot_spectra(fuv)
#temp = plt.xlim([1458,1461])
#temp = plt.xlim([1526,1530])
#temp = plt.xlim([1500,1530])
temp = plt.xlim([1500,1630])
plt.ylim([0,1])
Out[184]:
The spectra from 2011-04-11 and 2011-04-18 have the lowest flux and one of them (the blue one in the plot above) is also the spectrum with the largest wavelength offset compared to he mean (the other is only marginally shifted compared with the mean).
Interestingly, while the continuum in those two exposures has only about a third of the flux in the other exposures, the flux in the $H_2$ lines is virtually unchanged. Lines that are presumably stellar in origin (C IV, Mg II) are also the weakeast when the continuum is weak, but then again we would expect the lines to fluctuate with the continuum if both are powered by accretion.
In [22]:
s = 'LOS_LIMB' # change to other fits header keyword for Sun angle, Moon angle etc.
for f in glob.glob(datapath+'*20_jif.fits'):
hdulist = fits.open(f)
print '{0}: limb angles in the 3 exp: {1} {2} {3}'.format(f, hdulist[1].header[s],
hdulist[2].header[s], hdulist[3].header[s])
hdulist.close()
I checked the jit and jif files using the code above and by manually reading the entries. These files contain the jitter information from the FGS. There is nothing special about those two exposures. They have the same range of jitter, the same shape of jitter, the same range of angles between LOS and Sun, Moon, Earth's limb and brightness of the closest Earth's limb.
In [23]:
# Look at the corrtag file in more detail and see if there is any hint of an instrumental
# cause for the different flux levels.
n = []
meandark = []
dqflags = []
n0=[]
n512=[]
n4096=[]
nwave=[]
flist = glob.glob(datapath+'*corrtag_a.fits')
for f in flist:
events = Table.read(f, hdu='EVENTS')
n.append(len(events))
dqflags.append(set(events['DQ']))
n0.append((events['DQ']==0).sum())
n512.append((events['DQ']==512).sum())
n4096.append((events['DQ']==4096).sum())
nwave.append((events['WAVELENGTH'] > 1500).sum())
timeline = Table.read(f, hdu='TIMELINE')
meandark.append(np.mean(timeline['DARKRATE']))
In [24]:
# compare spectrum with low and high cont rate
evhigh = Table.read(flist[7], hdu='EVENTS')
evlow = Table.read(flist[10], hdu='EVENTS')
hist = plt.hist(evhigh['RAWY'] - evhigh['YCORR'])
hist = plt.hist(evlow['RAWY'] - evlow['YCORR'])
I looked a features in the TW Hya spectrum. Here, I find a difference of about 3 pixel. This means that the wavelengths will be wrong by 6-3=3 pixel = 6 km/s. This is well within the COS specifications: http://www.stsci.edu/hst/cos/documents/handbooks/current/ch08.Acquisitions09.html
A centering error of 0.1 arcsec would cause an wavelength shift around 15 km/s, so in our case the shift in dispersion direction must be well below that.
The throughput of the COS primary science aperture (PSA) is almost flat with in 0.4 arcsec, so this small change in centering cannot explain why the count rate a factor of three lower. (See the link above and http://www.stsci.edu/hst/cos/documents/isrs/ISR2010_09.pdf).
One might think that the centering worked in dispersion direction, but not in cross-dispersion direction. However, it seems very unlikely, that we can center in dispersion direction to within 0.05 arcsec, but are off by > 1 arcsec (as required to explain the flux drop) in cross-dispersion direction.
Furthermore, we see in the spectrum that the lower count rate occurs only in the the continuum and in the stellar lines, but not in $H_2$. Assuming that the reduction in flux is due to a centering error > 1 arcsec, then this means that the $H_2$ emission is specially extended by at least that much. Herczeg et al. (2002, http://iopscience.iop.org/0004-637X/572/1/310/fulltext/55255.text.html#sc3.2) find in their STIS observations that the $H_2$ is not extended more than 0.05 arsec. (They look at the line profile perpendicular to the dispersion in their article, so the limit valid for an extension in cross-dispersion direction, but if the flux was really extended by 1 arcsec, that would have lead to a decreased spectral resolution as well.) Given all we know about the excitation of $H_2$ in the UV it is impossible that the $H_2$ emission is this wide.
In [25]:
hist = plt.hist(evhigh['XCORR'][(evhigh['DQ'] == 512) & (evhigh['PHA']==0)])
hist = plt.hist(evlow['XCORR'][(evlow['DQ'] == 512) & (evlow['PHA']==0)])
The data quality flag 512 marks events with an unusual PHA. Most, but not all of those events are background events, e.g. cosmics, which are independent of the count rate. Thus, in orbits with a lower stellar flux, the fraction of cosmics is higher and in regions of the detector where no stellar signal is observed (last bin in the histrogram) the absolut number is almost identical. However, some source events (the ISR cited above says <2%) are also flagged. This explains that the number of flagged evens in spectral regions with a strong signal (e.g. X = 6000) in the histogram is about three times as strong in the orbit with high stellar fluxes as in the orbit with low stellar fluxes.
The orbits with the lower fluxes are also shifted in wavelength compared to the average. Instrumental effects or pointing error could explain both facts at the same time. Yet, there is no hint of either of those. Indeed, the wavelength calibration is good to 6 km/s in all observations, which is better than typical in COS observations. We conclude that the causal relation in the other way around: It is not a poiting offset that causes a lower flux, but the real, lower stellar flux, causes a slightly less accurate poiting. This change is so slight, that it will not effect the flux scale, but can cause the 6 km/s shift in wavelength; again, this value is well within the expected uncertainty for COS. So, we can trust the flux calibration and just shift the wavelength to make the feature match.
We know from Herzceg et al. (2002,2004) that the $H_2$ lines are emitted from the upper layers of the disk. In TW Hya we see the disk face-on, so those lines don't have a line-of-sight velocity and if they had, this velocity should not change at all. Thus, we cross-correlate several regions with many $H_2$ lines to measure the relative wavelength offset between the different spectra.
In [193]:
# most of those are regions with well separated H2 lines, the last one looks different
# (the lines are not well separated)
# but I did not find any atomic lines in CHIANTI, so it's probably mostly H2, too.
wranges = np.array([[1385,1415], [1430, 1458], [1458, 1475], [1458,1510], [1510,1530],
[1585,1605], [1605,1625], [1654,1660]]) * u.Angstrom
shifts = np.zeros((wranges.shape[0], len(fuv)))
for i, w in zip(range(wranges.shape[0]), wranges):
shifts[i] = spectrum.spectrum.xcorr(fuv, w)
In [194]:
print shifts.mean(axis=0)
print shifts.std(axis=0)
In [28]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(wranges.mean(axis=1).value, shifts, 'o')
plt.xlabel('cental wavelength [$\AA$]')
plt.ylabel('wavelength shift [km/s]')
Out[28]:
In [179]:
fig = plt.figure(figsize = (12,3))
fig = plot_spectra(nuv)
temp = plt.xlim([2530, 2590])
temp = plt.xlim([2650, 2700])
temp = plt.xlim([2770, 2820])
temp = plt.ylim([0,3e-1])
In [30]:
# Repeat for NUV
wrangesNUV = np.array([[2560,2570], [2575,2580], [2667, 2680],
[2790, 2793]]) * u.Angstrom
shiftsNUV = np.zeros((wrangesNUV.shape[0], len(nuv)))
for i, w in zip(range(wrangesNUV.shape[0]), wrangesNUV):
shiftsNUV[i] = spectrum.spectrum.xcorr(nuv, w)
In [31]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(wrangesNUV.mean(axis=1).value, shiftsNUV, 'o')
plt.xlabel('cental wavelength [$\AA$]')
plt.ylabel('wavelength shift [km/s]')
plt.ylim([-50,+50])
Out[31]:
Use Gaussfits of lines (H_2 or ISM - Mg II). The center of the line is more accurate than 1 pixel. Here is a list of lines to try:
In [139]:
fig = plt.figure(figsize = (12,3))
fig =plot_spectra(nuv)
temp = plt.xlim([2575, 2580])
#temp = plt.xlim([2650, 2700])
#temp = plt.xlim([2770, 2820])
temp = plt.ylim([0,.3])
#temp = plt.xlim([2796,2796.5])
#temp = plt.xlim([2793,2800])
#temp = plt.xlim([2802,2804])
In [33]:
from utils import LineModel, LineAbsModel
In [35]:
line = nuv[5].slice_disp([2537.*u.AA, 2542*u.AA])
g_init = LineModel(amplitude=np.mean(line.flux.value), mean=np.mean(line.disp.value),
stddev=.35, const=np.mean(line.flux.value))
f2 = fitting.LevMarLSQFitter()
g = f2(g_init, line.disp.value, line.flux.value)
In [136]:
def fitline(spectrum, rest, bounds):
'''Fit the position and shape of a line
To-Do:
- bounds = 1 elem -> make symmetric interval
- allow default model to be set or at least control input to model
- This could be an object that interacts with a spectrum object
- This could be filili (make it part of spectrum?)
- have way g, f2 are still available in object to help quality control and debug
'''
line = spectrum.slice_disp(bounds)
g_init = LineModel(amplitude=np.mean(line.flux.value)*1e12, mean=np.mean(line.disp.value),
stddev=.35, const=np.mean(line.flux.value))
f2 = fitting.SLSQPLSQFitter()
g = f2(g_init, line.disp.value, line.flux.value*1e12)
return (g.mean * line.disp.unit).to(u.km/u.s, equivalencies=u.doppler_optical(rest))
def spectrafitline(spectra, rest):
'''make this interface similar to xcorr'''
shift = np.zeros((len(spectra))) * u.km/u.s
for i in range(len(spectra)):
shift[i] = fitline(spectra[i], rest, [rest-2*u.AA, rest+2*u.AA])
return shift
In [140]:
shift2539 = spectrafitline(nuv, 2539.7628 * u.AA)
shift2578 = spectrafitline(nuv, 2576.9 * u.AA)
In [141]:
shift2539, shift2578
Out[141]:
In [125]:
def lineshift(spec, region, **kwargs):
# normalize flux to about one, otherwise fit won't work
line = spec.slice_disp(region)
g_init = LineAbsModel(amplitude=np.mean(line.flux.value)*3*1e12,
stddev=.8, const=.1,
amplitude2=np.mean(line.flux.value)*2.5*1e12, stddev2=0.05, **kwargs)
f2 = fitting.SLSQPLSQFitter()
g = f2(g_init, line.disp.value, line.flux.value*1e12) # Is there a way to ignore stdout?
return g.mean, g.mean2
MgII2796shift = np.zeros((len(nuv), 2))
for i, spec in enumerate(nuv):
MgII2796shift[i,:] = np.array(lineshift(spec, [2793*u.AA,2800*u.AA],
mean=2796.55, mean2=2796.25)).reshape(-1)
MgII2803shift = np.zeros((len(nuv), 2))
for i, spec in enumerate(nuv):
MgII2803shift[i,:] = np.array(lineshift(spec, [2800*u.AA,2807*u.AA],
mean=2803.7, mean2=2803.45)).reshape(-1)
#print (g.mean2 * line.disp.unit).to(u.km/u.s, equivalencies=u.doppler_optical(2796.3521*u.AA))
#print g_init
#print g
In [126]:
print (MgII2796shift[:,:] * line.disp.unit).to(u.km/u.s,
equivalencies=u.doppler_optical(2796.3521*u.AA))
In [127]:
print (MgII2803shift[:,:] * line.disp.unit).to(u.km/u.s,
equivalencies=u.doppler_optical(2803.5310*u.AA))
Here we fit the Mg II lines in the NUV. We use a model of a Gaussian emission line and a Gaussian absorption line. Presumably, the absorbtion is due to the ISM and should be exactly at the rest wavelenght. Since Mg II is a doublet, we perform the fit twice and compare the results. They all agree within 1 km/s. There is a larger variability between the peaks of the emission component, but that's not as reliable, because the line is not always Gaussian, while the absorption line is fairly narrow and the position of the minimum should be well determined.
Now there are two tasks open:
In [42]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [42]:
In [42]:
In [43]:
# Idea 1
ISMwave = [1526.7090, 2796.3521, 2803.5310]
resISM[1,:]/resISM[2:] #-> both Mg II fits agree to <7%
#some 1526.7090 are not fitted, or don't even show abs.
#line is much weaker.
# Si II 1526.7090 probably is not an ISM abs line, because the other member of the doublet
# (Si II 1533.4320) has no abs, although it is stronger. The abs. is probably something else - H2?
resISM[0,:]/resISM[2:] #-> little agreeement with NUV
#-> Let's see if relative shifts are compatible with the H_2 lines shifts
# Molat R(12) 3-10 is 1615.439
# http://amrel.obspm.fr/molat/index.php
h2lines=[1386.519, 1387.366, 1396.221, 1402.648, 1431.014, 1442.867, 1446.121, 1453.101,
1555.891, 1556.869]
#-> H2 and ISM on NUV move in parallel.
#-> ISM on FUV is not present on every obs and, where present, it is similar to ISM NUV and H2,
# but it's not clear if it is identical.
#This is an absolute shift that should make everything right
nuvshift = resISM[1:,:].mean(axis=0)
# This is a relative shift, that should put all FUVs on top of each other.
# Absolute shift can only be done, when enough H2 lines are identified.
# TBD: There seems to be a shift <-> lambda relation
# Maybe the grating parameters are not quite right?
# Watch out for that and come back after more H2 lines are securely identified.
# So far, I am not sure if there is a systemtic difference or
# if just some of the lines are misidentified.
fuvshift = (b1-b2).mean(axis=1)
for band, bandshift in zip([fuv, nuv],[fuvshift, nuvshift]):
for tab, shift in zip(band, bandshift):
tab.WAVELENGTH -=shift
plot_spectra(fuv, nuv)
#taken from stackoverflow
y_formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)
plt.gca().xaxis.set_major_formatter(y_formatter)
I played around a lot with different input parameters and tried to see what happens. From this I clearly see that the fitter is just not very good. Some parameters are not changed at all from their starting value (not always the same parameters, that depends on the starting value). If I want to use this, the fitters have to improve or I have to work out how the sets the fit parameters (required accuaracy, step size etc.) or normalize the flux before I fit it or something like that and hope that that is numerically more stable.
For now, I just see that the velocity shift is only of the order <6 km/s in the FUV and <20 km/s in the NUV. If I don't apply any RV shift for now (fix that after AAS and CS) I don't make much of an error. For posters that is fine.
In [186]:
spectra = []
for f,n in zip(fuv, nuv):
disp = np.hstack([f.disp.value, n.disp.to(f.disp.unit).value]) * f.disp.unit
spectra.append(spectrum.spectrum.coadd_simple([f,n], dispersion=disp, bounds_error=False))
In [191]:
tab.flux.shape
Out[191]:
In [190]:
# The line fluxes (I tired different H_2) line differ only by 10% between observations, the continuum by a factor 3. Maybe that is not a flux cal error, but is real?
for tab in fuv:
cor = np.correlate(fuv[0].flux[0], tab.flux[0],mode = 'full')
print tab.dispersion[0][1000], np.argmin(np.abs(tab.dispersion[0]-1575.)), len(tab.flux[0]), np.argmax(cor)
In [ ]:
def clip_expand(data, n=5):
'''mask at least n values around each masked value
Use this to prevent the wings of the H_2 emission lines to contribute'''
mask = np.array(data.mask, dtype=np.float)
kernel = np.ones(2*n)
data.mask = np.convolve(mask, kernel, mode='same') > 0
return data
In [ ]:
from astropy.stats import sigma_clip
line=spectra[9].slice_disp([1630*u.AA, 1680*u.AA])
plt.plot(line.disp.value, line.flux.value)
fclip = sigma_clip(line.flux.value, iters=None)
plt.plot(line.disp.value, fclip)
plt.plot(line.disp.value, clip_expand(fclip,n=20))
In [ ]:
cliprange = np.array([[1400,1450], [1450,1500], [1500, 1550], [1580,1630], [1630,1680],
[1680, 1730],
[2560, 2580], [2650, 2700], [2777, 2792]]) * u.Angstrom
In [ ]:
cont_flux = np.zeros((10,9))
for i, spec in enumerate(spectra):
for j in range(len(cliprange)):
line = spec.slice_disp(cliprange[j,:])
clipped = sigma_clip(line.flux.value, iters=None)
cont_flux[i,j] = np.ma.median(clip_expand(clipped, n=20))
In [ ]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(cliprange.mean(axis=1).value, cont_flux.T, 'o')
plt.xlabel('central wavelength [$\AA$]')
plt.ylabel('flux [erg/s/cm$^2$]')
In [ ]:
from scipy.optimize import curve_fit
import astropy.constants as const
def planck(T, wave):
'''Planck function
Parameters
----------
Input units are hardcoded, because scipy's curvefit chops them off
Returns
-------
spectral radiance per unit projected area of emitting surface, per unit solid angle,
per Ang
'''
wave = wave * u.AA
T = T * u.K
factor = 2.*const.h*const.c**2. / wave**5.
expfactor = np.exp(const.h*const.c / (wave*const.k_B*T)).decompose() - 1.
return factor.to(u.erg/u.AA/u.second/u.cm**2) / expfactor
def fitfunc(wave, amp, T):
# 4 pi to go from per unit solit angle to full sphere
return amp * 4 * np.pi* planck(T, wave).value
fitres= np.zeros((cont_flux.shape[0],2))
for i in range(len(spectra)):
a,b = curve_fit(fitfunc, cliprange.mean(axis=1).value, cont_flux[i,:], p0=[1e-21, 2e5])
fitres[i,:] = a
In [ ]:
# Assuming R_* = 0.7
plt.plot((fitres[:,0] * (57*u.pc)**2 / (0.7*const.R_sun)**2).decompose().value)
plt.ylabel('Filling factor of stellar surface')
plt.xlabel('Number of spectrum')
In [ ]:
plt.plot(fitres[:,1])
plt.ylabel('Black-body temperature [K]')
plt.xlabel('Number of spectrum')
In [ ]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(cliprange.mean(axis=1).value, cont_flux.T, 'o')
plt.xlabel('central wavelength [$\AA$]')
plt.ylabel('flux [erg/s/cm$^2$]')
ax.set_color_cycle(colorlist)
for i in range(10):
ax.plot(cliprange.mean(axis=1).value,
fitfunc(cliprange.mean(axis=1).value, fitres[i,0].T, fitres[i,1]))
The plot above shows black-body fits to the "continuum" measurements that are performed using sigma clipping. All temperatures are around 12000 K, which indicates that my numbers are not totally bogus. However, the fits are not particualrly good. That might just indicate that there are other continuum features in the data (see Kevin France's work) or that I did not subtract the emission lines well enough yet.
I think the temperature scale is around 12000 K and the filling factor is around 0.1%, but the fits are not good enough to believe the $\Delta T$ or $\Delta$ amplitude. Not surprisingly, the spectra in the low state have the smallest filling factors, but also relatively high temperatures. It almost seems that more accretion cools down the chromosphere.
In [ ]:
CIV = np.array([1548.19, 1550.775]) * u.Angstrom
def eq(wave):
return {'equivalencies': u.doppler_optical(wave)}
label = ['t = 0', '1.5 h', ' 2 h', ' 1 d', ' 3 d', ' 5 d', ' 8 d', '16 d', '20 d', '26 d']
def plot_doublet(spectra, waves, title='', bin_up=1):
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111, axisbg='0.95')
for s, c, l in zip(spectra, colorlist, label):
spec = s.bin_up(bin_up)
linespec = spec.slice_rv([-550*u.km/u.s, 1500*u.km/u.s], waves[0])
line = plt.plot(linespec.disp.to(u.km/u.s, **eq(waves[0])).value,
linespec.flux.value *1e12, color = c, label=l, lw=2)
ax.set_xlabel(r'v [km/s]', fontsize='x-large')
ax.set_ylabel(r'flux $\left[10^{-12} \frac{\rm{erg}}{\rm{s\;cm}^2 \AA}\right]$', fontsize='xx-large')
ax.plot([0,0], ax.get_ylim(), 'k:', lw=2)
ax.plot(waves[[1,1]].to(u.km/u.second, **eq(waves[0])).value, ax.get_ylim(), 'k:', lw=2)
ax.set_xlim([-550, 1500])
ax.set_ylim([0, 1.6])
ax.legend()
ax.tick_params(labelsize='xx-large')
ax.set_title(title, fontsize='xx-large')
return fig
fig = plot_doublet(fuv, CIV, 'C IV doublet')
fig.subplots_adjust(right=0.96, bottom=0.15)
#fig.savefig(plotpath+'CIV.png')
plt.savefig(posterdir + 'CIV.png', dpi=600, max_dpi=None, facecolor='0.8', edgecolor='none')
In [ ]:
def linerate(spectrum, waves, bin_up=1):
# Make a copy
spec = spectrum.bin_up(bin_up)
contregion = spec.slice_rv([-5000*u.km/u.s,-500*u.km/u.s], waves[0])
continuum = sigma_clip(contregion.flux.value, iters=1)
spec[spec.fluxname] -= np.ma.mean(continuum)
shifted = spec[np.ones((len(spec)), dtype=np.bool)]
shifted.shift_rv(-waves[1].to(u.km/u.s, **eq(waves[0])))
return (spec.flux / shifted.interpol(spec.disp, bounds_error=False).flux).decompose()
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111, axisbg='0.95')
for s, c, l in zip(fuv, colorlist, label):
s = COStools.spectrum.COSspectrum(s, dispersion=s.dispersion)
s = s.bin_up(2)
line = plt.plot(s.disp.to(u.km/u.s, **eq(CIV[0])).value,
linerate(s, CIV).value, color = c, label=l, lw=2)
ax.set_xlabel(r'v [km/s]', fontsize='xx-large')
ax.set_ylabel(r'flux ratio', fontsize='xx-large')
ax.plot([0,0], ax.get_ylim(), 'k:', lw=2)
ax.plot([-1000,1000], [2,2], 'k:', lw=4)
#ax.set_xlim([-500, 1550])
ax.set_xlim([-50, 550])
ax.set_ylim([0,4])
ax.legend()
ax.tick_params(labelsize='xx-large')
ax.grid()
ax.bar(175,4,50, fc='0.8', ec='none', zorder=2.1,alpha=.7)
ax.text(185, 2, r'blend with $H_2$ line', rotation='vertical', va='center', fontsize='xx-large')
ax.set_title('C IV doublet line ratio', fontsize='xx-large')
fig.subplots_adjust(right=0.98, bottom=0.15)
plt.savefig(posterdir + 'CIV_ratio.png', dpi=600, max_dpi=None, facecolor='0.8', edgecolor='none')
In [ ]:
MgII = np.array([2796.3521, 2803.5310]) * u.Angstrom
fig = plot_doublet(nuv, MgII, 'Mg II doublet')
fig.subplots_adjust(right=0.97)
plt.savefig(posterdir + 'MgII.png', dpi=600, max_dpi=None, facecolor='0.8', edgecolor='none')
In [ ]:
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
for s, c, l in zip(nuv, colorlist, label):
s = COStools.spectrum.COSspectrum(s, dispersion=s.dispersion)
s = s.bin_up(2)
line = plt.plot(s.disp.to(u.km/u.s, **eq(MgII[0])).value,
linerate(s, MgII).value, color = c, label=l, lw=2)
ax.set_xlabel(r'v [km/s]', fontsize='xx-large')
ax.set_ylabel(r'flux ratio', fontsize='xx-large')
ax.plot([0,0], ax.get_ylim(), 'k:', lw=2)
ax.plot([-1000,1000], [2,2], 'k:', lw=4)
#ax.set_xlim([-500, 1550])
ax.set_xlim([-300, 300])
ax.set_ylim([0,4])
ax.legend()
ax.tick_params(labelsize='xx-large')
ax.grid()
#fig.savefig(plotpath+'optdepth.png')
In [ ]:
SiII = np.array([1526.707, 1533.431]) * u.Angstrom
fig = plot_doublet(fuv, SiII, 'Si II doublet - not as strong. H2 contamination?', bin_up=4)
ax = fig.axes[0]
ax.set_xlim([-200,1900])
ax.set_ylim([0,0.2])
#fig.savefig(plotpath+'CIV.png')
In [ ]:
HeII = np.array([1640.4771]) * u.Angstrom
def plot_line(spectra, waves, title, bin_up=1):
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111,axisbg='0.95')
for spec, c, l in zip(spectra, colorlist, label):
linespec = spec.slice_rv([-550*u.km/u.s, 1500*u.km/u.s], waves[0])
linespec = linespec.bin_up(bin_up)
line = ax.plot(linespec.disp.to(u.km/u.s, **eq(waves[0])).value,
linespec.flux.value *1e12, color = c, label=l, lw=2)
ax.set_xlabel(r'v [km/s]', fontsize='x-large')
ax.set_ylabel(r'flux $\left[10^{-12} \frac{\rm{erg}}{\rm{s\;cm}^2 \AA}\right]$', fontsize='xx-large')
ax.plot([0,0], ax.get_ylim(), 'k:', lw=2)
ax.set_xlim([-550, 750])
ax.set_ylim([0, 2.6])
ax.legend()
ax.tick_params(labelsize='xx-large')
ax.set_title(title, fontsize='xx-large')
return fig
fig = plot_line(fuv, HeII, 'He II')
fig.subplots_adjust(right=0.98, bottom=0.15)
plt.savefig(posterdir + 'HeII.png', dpi=600, max_dpi=None, facecolor='0.8', edgecolor='none')
In [ ]:
fig = plot_line(fuv, HeII, 'He II', bin_up=4)
fig.axes[0].set_ylim([0, .4])
In [ ]:
fig = plot_line(fuv, [CIV[1]], 'C IV blue', bin_up=2)
fig.axes[0].set_ylim([0, 1])
In [ ]:
SiIV = np.array([1393.757]) * u.Angstrom
fig = plot_line(fuv, SiIV, '2 * H_2 + Si IV', bin_up=2)
fig.axes[0].set_ylim([0, .6])
In [ ]:
linespec = fuv[5]
line = plt.plot(linespec.disp.to(u.km/u.s, **eq(CIV[1])).value,
linespec.flux.value *1e12, color = 'k', label=l, lw=2)
line = plt.plot(linespec.disp.to(u.km/u.s, **eq(HeII[0])).value,
linespec.flux.value *1e12*1.4, color = 'r', label=l, lw=2)
plt.xlim([-550, 750])
In [ ]:
CI = np.array([1656.3]) * u.Angstrom # Multiplet - only strongest line listed for quicklook
fig = plot_line(fuv, CI, 'C I', bin_up=2)
fig.axes[0].set_ylim([0, .6])
The C I lines look stronger at 26 d than ever. The ratios also change slightly between exposures (e.g. the first 2 vs the later three lines in 1d<->8d), but I am not sure that is statistically significant.
In [ ]:
SV = np.array([1501.776]) * u.Angstrom # Multiplet - only strongest line listed for quicklook
fig = plot_line(fuv, SV, 'S V ? or just some other H_2 line?', bin_up=2)
fig.axes[0].set_ylim([0, .3])
In [ ]:
#line = fuv[0].slice_disp([1498*u.AA,1502*u.AA])
#line = fuv[0].slice_rv([-1000*u.km/u.s, 2500*u.km/u.s], HeII[0])
line = nuv[0].slice_rv([-1600*u.km/u.s, 1300*u.km/u.s], MgII[0])
#contregion = spec.slice_rv([-5000*u.km/u.s,-500*u.km/u.s], waves[0])
continuum = sigma_clip(line.flux.value, iters=None, cenfunc=np.ma.median)
#spec[spec.fluxname] -= np.ma.mean(continuum)
continuum.mask.sum()
# something is persistent here and it might not be cigma_clip
# sigma_clip should use ma versions.
# sigma_clip should improve documentation.
In [ ]:
plt.plot(line.disp.value, line.flux.value)
plt.plot(line.disp.value[~continuum.mask], line.flux.value[~continuum.mask])
In [ ]:
# Now is fully automatic, but that requires either tuning kwargs or have a large section of
# good continuum
def find_continuum(spectrum, iters=None, **kwargs):
# This works well for isoalted lines. Might want alternative ways for longer spectra.
# Maybe look for column "continuum" in table?
# add flux_norm and flux_subtracted as properties?
continuum = sigma_clip(spectrum.flux.value, iters=iters, **kwargs)
return np.ma.median(continuum)*spectrum.flux.unit, np.std(continuum)*spectrum.flux.unit
def sum_flux(spectrum):
return (spectrum.flux[:-1] * np.diff(spectrum.disp)).sum()
def line_flux(spectrum, **kwargs):
return sum_flux(spectrum) - find_continuum(spectrum, **kwargs)[0] * np.abs(spectrum.disp[0] - spectrum.disp[-1])
def EW(spectrum, **kwargs):
# does not really make sense for longer spectra...
return -line_flux(spectrum, **kwargs) / find_continuum(spectrum, **kwargs)[0]
In [ ]:
# Here is the exact same thing, except that line region and background region can be provided as
# separate inputs. I try this out, too, and see how it works and then put one of them into my spectrum
# class.
# I almost don't need the sigma clip here any longer if I hand select continuum regions.
# Other idea is continuum column
def find_continuum(cont_spec, iters=None, **kwargs):
# This works well for isoalted lines. Might want alternative ways for longer spectra.
# Maybe look for column "continuum" in table?
# add flux_norm and flux_subtracted as properties?
continuum = sigma_clip(cont_spec.flux.value, iters=iters, **kwargs)
return np.ma.median(continuum) * cont_spec.flux.unit, np.std(continuum)*cont_spec.flux.unit
def sum_flux(spectrum):
return (spectrum.flux[:-1] * np.diff(spectrum.disp)).sum()
def line_flux(spectrum, cont_spec, **kwargs):
return sum_flux(spectrum) - find_continuum(cont_spec, **kwargs)[0] * np.abs(spectrum.disp[0] - spectrum.disp[-1])
def EW(spectrum, cont_spec, **kwargs):
# does not really make sense for longer spectra...
return -line_flux(spectrum, cont_spec, **kwargs) / find_continuum(cont_spec, **kwargs)[0]
In [ ]:
(line.flux[:-1] * np.diff(line.disp)).sum()
In [ ]:
find_continuum(line, cenfunc=np.ma.median)
In [ ]:
line_flux(line, cenfunc=np.ma.median)
In [ ]:
EW(line, cenfunc=np.ma.median)
In [ ]:
ewtab = astropy.table.Table()
fluxtab = astropy.table.Table()
#for n, r in zip(['C IV doublet', 'He II', 'C I multiplet', 'Mg II doublet', 'H_2 P(3) 0-5',
# 'H_2 P(2) 0-6'],
# [[1540*u.AA, 1555*u.AA], [1635.*u.AA, 1645.*u.AA], [1650.*u.AA, 1660.*u.AA],
# [2790*u.AA, 2810*u.AA], [1402*u.AA, 1404*u.AA],[1459.2*u.AA, 1461*u.AA]]):
for n, r, c in zip(['C IV doublet', 'He II', 'C I multiplet', 'Mg II doublet', 'H$_2$ P(3) 0-5',
'H_2 P(2) 0-6'],
[[1545.5*u.AA, 1554*u.AA], [1639.*u.AA, 1644.*u.AA], [1655.*u.AA, 1660.*u.AA],
[2793*u.AA, 2807*u.AA], [1402*u.AA, 1404*u.AA],[1459.5*u.AA, 1461*u.AA]],
[[1540*u.AA, 1546*u.AA], [1635.*u.AA, 1639.*u.AA], [1650.*u.AA, 1656.*u.AA],
[2785*u.AA, 2795*u.AA], [1404*u.AA, 1406*u.AA],[1458*u.AA, 1460*u.AA]]):
ew = np.zeros(10) * u.AA
flux = np.zeros(10) * u.erg/u.cm**2/u.second
for i in range(10):
line = nf[i].slice_disp(r)
cont = nf[i].slice_disp(c)
ew[i] = EW(line, cont, iters=None, cenfunc=np.ma.median)
flux[i] = line_flux(line, cont, iters=None, cenfunc=np.ma.median)
ewtab.add_column(astropy.table.Column(np.array(ew), n, unit=ew.unit))
fluxtab.add_column(astropy.table.Column(np.array(flux), n, unit=ew.unit))
In [ ]:
line=spectra[7].slice_disp([1459.5*u.AA, 1461.*u.AA])
plt.plot(line.disp.value, line.flux.value)
fclip = sigma_clip(line.flux.value, iters=None, cenfunc=np.ma.median)
plt.plot(line.disp.value, fclip)
In [ ]:
for col in ewtab.colnames:
plt.plot(ewtab[col]/np.min(ewtab[col]), label=col)
In [ ]:
for col in fluxtab.colnames:
plt.plot(fluxtab[col]/np.max(fluxtab[col]), label=col)
In [ ]:
temp = plt.plot(cont_flux)
In [ ]:
smarts = astropy.table.Table.read('/data/hguenther/TWHya/SMARTSlc.dat', format='ascii')
filters = ['B','V','R','I','J','H']
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(111)
for i, filt in enumerate(filters):
ax.errorbar(tab['JD']-2455650, tab[filt]+i, tab[filt+'_err'], label=filt)
ax.set_xlabel('JD - 2455650')
ax.set_ylabel('mag [arb. zero point]')
ax.legend()
ax.set_xlim([0,80])
ax.set_ylim([6,-1])
fig.subplots_adjust(top = .96, bottom = .16, left=.15, right=.95)
#plt.draw()
#fig.savefig('/data/hguenther/Dropbox/my_articles/TW_Hya_UV/SMARTSlc.png', bb_inches='tight')
In [ ]:
smarts = astropy.table.Table.read('/data/hguenther/TWHya/SMARTSlc.dat', format='ascii')
fig = plt.figure(figsize=(5,4))
ax = fig.add_axes([.2, .13, 0.75, .75], axisbg='0.95')
ax.errorbar(smarts['JD']-2455650, -smarts['B']/2.5-11., tab['B_err']/2.5, label='B (scaled)', elinewidth=2, capthick=2)
ax.errorbar(smarts['JD']-2455650, -smarts['R']/2.5-11.2, tab['R_err']/2.5, label='R (scaled)', elinewidth=2, capthick=2)
for col in ['C IV doublet', 'He II', 'C I multiplet', 'H$_2$ P(3) 0-5']:
ax.plot(fuvtime.jd-2455650, np.log10(fluxtab[col]), 'o-', label=col)
col='Mg II doublet'
ax.plot(nuvtime.jd-2455650, np.log10(fluxtab[col]), 'o-', label=col)
ax.plot(fuvtime.jd-2455650, np.log10(cont_flux[:,0]*10), 'o-', label='FUV cont\n(scaled)')
# comment NUV cont because plot looks messy with too many lines
#ax.plot(fuvtime.jd-2455650, np.log10(cont_flux[:,-1]*10), 'o-', label='NUV cont\n(scaled)')
ax.set_ylabel(r'log(flux) $\left[\frac{\rm{erg}}{\rm{s\;cm}^2}\right]$', fontsize='x-large')
ax.set_xlabel('JD - 2455650', fontsize='x-large')
ax.tick_params(labelsize='large')
ax.set_title('Lightcurves\n(continuum, B, R scaled for clarity)', fontsize='x-large')
ax.legend(loc='lower right')
ax.set_xlim([10,60])
ax.set_ylim([-13, -10.75])
plt.savefig(posterdir + 'lc.png', dpi=600, max_dpi=None, facecolor='0.8', edgecolor='none')
In [ ]:
'''Results from the analysis above
General
-------
tuvb :
All low FAP have also low p(lc = const),
but some low p as high FAP
Individual
----------
lbl207d0q_rawtag_b.fits : Unusual in the CIV acorr
lc: has emission event 200 < t < 300 (_a, _b and CIV), with increase in count rate 6-8 %
extra flux seems concentrated in peak of lines (H2 and CIV). This this true? Or does it just look like that, because cont. is more bins but lower numbers per bin?
'''
# plot all lc in subplots
for i, spec in enumerate(tuvb['file']):
t = atpy.Table(str(spec), hdu = 1)
# cut off wavelcal lamp
ypos = ev.tabmax(t, 'RAWY', bin = 10)plot_spectra(*split)
t = t.where((t.RAWY > ypos-40) & (t.RAWY < ypos+40))
xpos = ev.tabmax(t, 'RAWX', bin = 100)
#t = t.where((t.RAWX > xpos-250) & (t.RAWX < xpos+500))
ax = plt.subplot(6,5,i+1)
hist, bins = np.histogram(t.TIME/0.032, bins = np.arange(-.5,max(t.TIME/0.032), 120))
h = ax.plot(bins[1:], hist)
'''
plotting fuv_b reveals some strong trends in lcs, often spanning all 3 lcs per obs
In CIV only the strongest changes are visible.
That could have different reasons:
-they are driven by the cont.
-by the H2 lines
-they are due to instrumental effects, e.g. the changing FP-POS (howver, that should mean changing flux between obs and not during obs
-> acutally extract spectra with the pipeline on time scales of minutes and see what I find, if I integrate that over large (but consistent) energy ranges/ the H2 lines, the cont etc.
Timing To-Do:
simulations, analytical number on max periodicity
timing in minute-month scales
extract spectra with pipeline (see notes above)
write that down
Spectral To-Do:
Mg II should be 2:1 (same conf as C IV), but it is not
find Lamzin profiles (in HH on some disk?) and check what Anthoony did
-> Can that explain the line profiles? (No, if I remember correctly)
'''
m1 = np.zeros_like(n)
s1 = np.zeros_like(n)
for i in range(len(n)):
dist = scipy.stats.poisson(mean[i]).rvs(n[i])
m1[i] = np.mean(dist)
s1[i] = scipy.var(detrend_poly.detrend_poly(dist, deg = 3))
'''
Compare extraction of time split data with original processing
-> Looks good. No wavelength shift between individual files
-> can apply the wavelength shift determined from the summed spectrum
-> differences in cont and line flux are marginal.
-> Need to measure by fitting lines + cont
'''
import glob
import atpy
import sys
sys.path.append('/data/hguenther/Dropbox/my_articles/TW_Hya_UV/python')
import HSTplotdefs as pd
base = 'lbl207d0q'
ta = atpy.Table(pd.datapath + base + '_x1d.fits')
splitlist = glob.glob('/media/MAX/moritz/processing/HST/TW_Hya/calcos/'+base+'*x1d.fits')
split = []
for f in splitlist:
tab = atpy.Table(f)
split.append(tab)
plot_spectra(split)
for a in split:
a.FLUX[0,:].sum()
splitf = np.zeros((9,2,16384))
for i in range(len(split)):
splitf[i,:,:] = split[i].FLUX
Next step: Interpol SMARTs lightcurve and look at B<->CIV and similar
\acknowledgements Support for program GO-12315.01-A was provided by NASA through a grant from the Space Telescope Science Institute. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. %Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts.
This research made use of Astropy, a community-developed core Python package for Astronomy \citep{2013arXiv1307.6212T}.