I now calcaulte the branchign rations in Abgrall myself. Use that instead of Kavin France table, because I calcualte them for every line, evne those that are not in his table.
In [1]:
import glob
import sys
from copy import deepcopy
import re
import os
import shutil
from collections import OrderedDict
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, MaskedColumn
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
import spectrum.coadd
from spectrum.spectrum import Spectrum
# Dirty hack to import stuff that's in the same directory as this notebook.
sys.path.append('/melkor/d1/guenther/projects/TWHya/')
import code as TWHya
from H2 import H2
from utils import LineModel, LineAbsModel # file in same directory as this notebook.
%matplotlib inline
In [2]:
base_data_path = '/melkor/d1/guenther/projects/TWHya/'
datapath = base_data_path + 'COS/'
datastis = base_data_path + 'STIS/'
#plotdir = '/melkor/d1/guenther/Dropbox/my_talks/15_ESTEC/'
#plt.style.use(('/melkor/d1/guenther/soft/python/mplstylelib/highdpi.mplstyle'))
In [3]:
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']
labellist = ['t = 0', '1.5 h', ' 2 h', ' 1 d', ' 3 d', ' 5 d', ' 8 d', '16 d', '20 d', '26 d']
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 [4]:
filelist = glob.glob(datapath + '*sum.fits')
filelist.sort()
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 [5]:
len(filelist)
Out[5]:
In [6]:
plot_spectra(nuv)
In [ ]:
plot_spectra(nuv)
%%\documentclass[12pt,preprint]{aastex}
%\documentclass[manuscript]{aastex}
%% preprint2 produces a double-column, single-spaced document:
%\documentclass[preprint2]{aastex} \documentclass{emulateapj}
%% Sometimes a paper's abstract is too long to fit on the %% title page in preprint2 mode. When that is the case, %% use the longabstract style option. %% \documentclass[preprint2,longabstract]{aastex}
\usepackage{natbib} \citestyle{aa} %% If you wish, you may supply running head information, although %% this information may be modified by the editorial offices. %% The left head contains a list of authors, %% usually a maximum of three (otherwise use et al.). The right %% head is a modified title of up to roughly 44 characters. %% Running heads will not print in the manuscript style.
\shorttitle{UV emission from TW Hya} \shortauthors{G\"unther et al.}
\begin{document} %% Use \author, \affil, and the \and command to format %% author and affiliation information. %% Note that \email has replaced the old \authoremail command %% from AASTeX v4.0. You can use \email to mark an email address %% anywhere in the paper, not just in the front matter. %% As in the title, use \ to force line breaks. \title{UV emission from TW Hya}
\author{H.~M.~G\"unther} \affil{Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138} \email{hguenther@cfa.harvard.edu}
\begin{abstract} Abstract \end{abstract}\keywords{circumstellar matter -- infrared: stars -- Stars: formation -- Stars: pre-main sequence -- X-rays: stars}
In this section, we describe how the data was taken and how we reduced it.
TW~Hya was observed with \emph{HST}/COS for ten orbits in program ID 12315 in 2011. Table~\ref{tab:obsCOS} gives a detailed observation log. The setup of all orbits is almost identical. The target is aquired in the NUV with the G285M grating using one peak-up exposure in cross-dispersion direction and another one in dispersion direction to center TW~Hya properly. There are four science exposures, with an exposure time of about 10~min each. First, a G285M exposure (central wavelength 2676~\AA{}) is taken. The NUV channel provides a resolution of around 20,000 and consists of three detectors, each of which covers about 40~\AA{} of a spectrum. In the setting chosen those spectral reagions are centered on 2566, 2675, and 2795~\AA{}. The remaining three exposures are taken with the G160M FUV grating centered on 1577~\AA{}. The spectral resolution is similar, but the FUV channel consists of only two chips which cover roughly 1385-1555~AA{} and 1580-1750~\AA{}. To reduce the fixed-pattern noise on the detector, each of the three FUV exposures is taken with a different grating offset position (\texttt{FPPOS=2,3,4}). These three exposures are combined in the pipeline extraction process and table~\ref{tab:obsCOS} reports their summed exposure time. All data were taken in time-tag mode, where the arrival time of individual photons is recorded.
Since the time scale of the variability in the hot ion lines was not known before the observations, the cadence is chosen to cover a range of time scales. The first three observations were scheduled in consequtive orbits. Therefore, no target aquisition was required in orbit 2 and 3, which results in slightly longer exposure times (table~\ref{tab:obsCOS}). The remaining observations where spread out over the visibility window of TW~Hya covering time scales between one day and one month.
We retrieved the data from the Mikulski Archive for Space Telescopes (MAST), where the default COS pipeline is run on the data. We further processed the data using custom python routines. All code used in the analysis is available at https://github.com/hamogu/TWHya.
In [7]:
# output->LaTeX
COSobs = Table([Column(name='date', data=[n.meta['DATE-OBS'] for n in nuv]),
Column(name='time', data=[n.meta['TIME-OBS'] for n in nuv]),
Column(name='NUV-ID', data=['\\dataset{{ADS/sa.HST\\#{0}}}'.format(n.meta['EXPNAME']) for n in nuv]),
Column(name='NUV exp.time', data=[n.meta['EXPTIME'] for n in nuv], format='%4.0f'),
Column(name='FUV-ID', data=['\\dataset{{ADS/sa.HST\\#{0}}}'.format(n.meta['EXPNAME']) for n in fuv]),
Column(name='FUV exp.time', data=[n.meta['EXPTIME'] for n in fuv], format='%4.0f')
])
latexdict = deepcopy(astropy.io.ascii.latexdicts['AA'])
latexdict['tabletype'] = 'deluxetable*'
latexdict['units'] = {'NUV exp.time':'s', 'FUV ext.time':'s'}
COSobs.write(sys.stdout, format='ascii.aastex', latexdict=latexdict,
caption=r"\label{tab:obsCOS}Observations log of \emph{HST}/COS observations")
In [6]:
def xcorr_order(n, stisdata, cos, xcorrlim):
'''Cross-correlate one order of a STIS spectrum with COS data.
Parameters
----------
n : int
STIS order number
stisdata : `astropy.table.Table`
STIS data read from a ``*x1d.fits`` file
cos : spectrum or list of spectra
COS data
xcorrlim : list
Min and max of the wavelength interval (in Ang) of the STIS spectra that is used for cross-correlation.
``None`` defaults to 1 Ang from the edge of the order.
'''
stis = Spectrum({'WAVE': stisdata['WAVELENGTH'][n,:], 'FLUX': stisdata['FLUX'][n,:],
'ERROR': stisdata['ERROR'][n,:]},
dispersion='WAVE', uncertainty='ERROR')
stis['WAVE'].unit = u.Angstrom # unit in fits file is "Angstroems" which is not recognized
start = min(stis.disp) + 1.*u.AA if xcorrlim[0] is None else xcorrlim[0]
stop = max(stis.disp) - 1.*u.AA if xcorrlim[1] is None else xcorrlim[1]
base = stis.slice_disp([start, stop])
shifts = spectrum.spectrum.xcorr(base, cos, np.arange(-20,20)*u.km/u.s)
return shifts
In [7]:
stisdata = Table.read(datastis+'o59d01030_x1d.fits', hdu=1)
# I have no idea why this is read in as masked, since no value is actually masked?
#assert np.all([0 == data[c].mask.sum() for c in data.colnames])
stisdata = stisdata.filled()
# None means start 1 ang from the edges
# These ranges are selected by hand (see above) to avoid ion emission lines
xcorrlimits = {
2: [1654.*u.AA, 1660.*u.AA],
4: [None, 1625*u.AA],
5: [None, None],
6: [1585.*u.AA, None],
10: [None, 1530.*u.AA],
11: [None, None],
12: [None, None],
13: [None, None],
14: [None, None],
15: [None, None],
16: [1430*u.AA, None],
18: [None, None],
19: [None, None]
}
shifts = np.zeros((len(xcorrlimits), len(fuv))) * u.km/u.second
for i, n in enumerate(xcorrlimits.keys()):
shifts[i,:] = xcorr_order(n, stisdata, fuv, xcorrlimits[n])
In [8]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(shifts, 'o')
plt.xlabel('cental wavelength [$\AA$]')
plt.ylabel('wavelength shift [km/s]')
# The last two plots are at the lowest wavelength. There is clearly a systematic here.
# For the purposes of this paper, we only need to make sure that things are consistent, so we will use
# the first ten values.
print shifts[:-2,:].mean(axis=0)
print np.std(shifts[:,:], axis=0)
The wavelength scale in \emph{HST}/COS spectroscopy is usually calibrated with an internal calibration lamp, that flashes on and off during the observation on a part of the detector where it does not overlap with the science data. According to the COS instrument handbook, this calibration is good to about 15~km~s$^{-1}$ for the medium resolution gratings we used\footnote{http://www.stsci.edu/hst/cos/documents/handbooks/current/ch05.COS_Spectroscopy02.html}. We improve this wavelength calibration using features in the spectra. The FUV data contains hundreds of $H_2$ emission lines, which are observed to show a consistent velocity shift of $13.55\pm0.10$~km~s$^{-1}$ in \emph{HST}/STIS \citep{2002ApJ...572..310H}. In turn, this shift is fully consistent with the radial velocity of TW~Hya \citep{2002ApJ...572..310H}. We cross-correlate each of our ten FUV spectra with the STIS/E140M spectrum used in \cite{2002ApJ...572..310H} (\dataset{ADS/sa.HST#o59d01030}) in spectral regions that are dominated by $H_2$ lines (1385 to 1415, 1430 to 1458, 1458 to 1510, 1510 to 1530, 1585 to 1605, 1605 to 1625, and 1654 to 1660 \AA{}) and measure the relative velocity shift. The maximal difference between two spectra is $<8$~km~s$^{-1}$, indicating a better than expected wavelength stability for \emph{HST}/COS \citep[similar to the findings of][]{2013ApJS..207....1A}; for each individual spectrum, the standard deviation between the eight spectral regions is $<4$~km~s$^{-1}$. We shift each spectrum such that the $H_2$ lines are at rest. In the FUV channel, the accuracy is thus limited by the \emph{HST}/STIS calibration from \citet{2002ApJ...572..310H} to $3$~km~s$^{-1}$, assuming that the $H_2$ lines are still at rest with respect to TW~Hya. This is a resonable assumption, since they orignate in the upper layer of the accretion disk and we see this disk face-on.
In [9]:
stisdata = Table.read(datastis+'o59d01020_x1d.fits', hdu=1)
# I have no idea why this is read in as masked, since no value is actually masked?
#assert np.all([0 == data[c].mask.sum() for c in data.colnames])
stisdata = stisdata.filled()
# None means start 1 ang from the edges
# These ranges are selected by hand (see above) to avoid ion emission lines
xcorrlimits = [
(3, [2790.*u.AA, 2793.*u.AA]),
(6, [2667.*u.AA, 2680.*u.AA]),
(9, [2575.*u.AA, 2580.*u.AA]),
(9, [None, 2570.*u.AA])
]
shiftsnuv = np.zeros((len(xcorrlimits), len(nuv))) * u.km/u.second
for i, xcorrlim in enumerate(xcorrlimits):
shiftsnuv[i,:] = xcorr_order(xcorrlim[0], stisdata, nuv, xcorrlim[1])
In [10]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(shiftsnuv, 'o')
plt.xlabel('cental wavelength [$\AA$]')
plt.ylabel('wavelength shift [km/s]')
# The last two plots are at the lowest wavelength. There is clearly a systematic here.
# For the purposes of this paper, we only need to make sure that things are consistent, so we will use
# the first twn values.
print shiftsnuv.mean(axis=0)
print np.median(shiftsnuv, axis=0)
print np.std(shiftsnuv, axis=0)
In [11]:
shifts.value
Out[11]:
In [12]:
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?
print g.mean, g.mean2, g.amplitude, g.amplitude2, g.stddev, g.stddev2
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)
In [13]:
print (MgII2796shift[:,:] * spec.disp.unit).to(u.km/u.s,
equivalencies=u.doppler_optical(2796.3521*u.AA))
In [14]:
print (MgII2803shift[:,:] * spec.disp.unit).to(u.km/u.s,
equivalencies=u.doppler_optical(2803.5310*u.AA))
Unfortunately, the COS/NUV spectra cover a much smaller wavelength range and there are much fewer $H_2$ lines. We fit the sharp absorption line within the wider \ion{Mg}{2} emission lines at 2796.3521~\AA{} and 2803.5310~\AA{}. We assume that the absorption is due to interstellar material at rest and we shift the NUV spectra to the rest frame of TW~Hya using the values found here. The difference between the fitted values for both lines in the \ion{Mg}{2} doublet is less than 1~km~s$^{-1}$; the absolute uncertaity is about 4~km~s$^{-1}$, mostly due to the uncertainty in the radial velocity of TW~Hya.
We also correlate our NUV data with the STIS/E230M spectrum \dataset{ADS/sa.HST#o59d01020}, but we find a larger scatter up to 20~km~s$^{-1}$ between the different regions used for cross-correlation within a single spectrum. The difference between the values found from the \ion{Mg}{2} absorption lines and this method is also of the order 10-20~km~s$^{-1}$, presumably because there are only very few features in the noisy spectra used for cross-correlation.
In [15]:
# Units are lost in the fitting, so add km/s here again
#for spec, rv in zip(fuv, shifts.mean(axis=0)):
# spec.shift_rv(rv * u.km / u.s)
NUVshifts = np.zeros_like(MgII2796shift)
NUVshifts[:, 0] = (MgII2796shift[:,1] * spec.disp.unit).to(u.km/u.s,
equivalencies=u.doppler_optical(2796.3521*u.AA))
NUVshifts[:, 1] = (MgII2803shift[:,1] * spec.disp.unit).to(u.km/u.s,
equivalencies=u.doppler_optical(2803.5310*u.AA))
In [16]:
NUVshifts.mean(axis=1), np.median(shiftsnuv, axis=0)-13.55*u.km/u.s
Out[16]:
In [17]:
# Do the actual shift of the data
dv = shifts[:-2,:].mean(axis=0)
for i in range(10):
fuv[i].shift_rv(dv[i] - 13.55*u.km/u.s)
In [18]:
# plot a line after the shift to check that I did not screw up the signs.
plot_spectra(fuv)
plt.xlim([1454.5,1455.5])
# observed wavelength in Herczeg et al 2002
plt.plot([1454.892, 1454.892], [0,1], 'r')
plt.plot([1455.038, 1455.038], [0,1], 'r')
# Real wavelengths in Abgrall
plt.plot([1454.971, 1454.971], [0,1],'k')
plt.plot([1454.829, 1454.829], [0,1], 'k')
plt.ylim([0,.5])
Out[18]:
In [19]:
# Shift the NUV data.
dv = NUVshifts.mean(axis=1)*u.km/u.s
# use - dv instead of +dv here, because the Mg II fit gives the sign the other way around compared to xcorr
for i in range(10):
nuv[i].shift_rv(-dv[i] - 13.55*u.km/u.s)
In [20]:
plot_spectra(nuv)
plt.plot([2796.3521, 2796.3521],[0,.2]) # Line position from CHIANTI, which is 13.55 km/s from TW Hya at rest
plt.xlim([2795,2798])
Out[20]:
Photometric monitoring was performed with the ANDICAM instrument \citep{2003SPIE.4841..827D} on the SMARTS/CTIO 1.3m telescope operated by the SMARTS consortium. Exposures were taken nightly in queue mode. ANDICAM uses a Fairchild 447 2048x2048 CCD with 15-micron pixels for observations in KPNO-recipe Johnson-Kron-Cousins $BVRI$ filters. The field-of-view is about 6\arcsec{} on each side. TW~Hya was roughly centered in this field. $YJHK$ images are taken with a Rockwell 1024x1024 HgCdTe ``Hawaii'' Array with 18-micron pixels. The FOV is about 2.4\arcsec{} on the side. The $JHK$ filters are standard CIT/CTIO filters and the $Y$ filter is a 1.05-micron central bandpass filter. The ANDICAM pipeline subtracts the overscan bias and a zero frame before flatfielding the CCD images. Three exposures are taken in $BVRI$ each, with exposure times of 20, 8, 4 and 2~s, respectively. The IR images are corrected for bad pixels, cosmic-ray rejected and the six dither positions (exposure time 4~s each) are co-added. We perform aperture photometry on the resulting CCD and IR images, using an aperture radius of four times the FWHM of the PSF. In each case, the four brightest stars in the field are used for relative photometry. We discard the $K$ and $Y$ data due to the low signal of the comparison stars.
The pipeline processed data for this observation is available at \textbf{put in dataverse, where put my lightvurve? Dataverse or electronic table in article?}
In [9]:
#Merge FUV and NUV of each observation into one array for easier plotting
cos = []
for f, n in zip(fuv, nuv):
disp = np.hstack([f.disp.value, n.disp.value]) * f.disp.unit
cos.append(spectrum.coadd.coadd_simple([f,n], dispersion=disp, bounds_error=False))
cos[-1].meta = n.meta
In [10]:
plt.figure(figsize=(20,4))
plot_spectra(cos)
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 [11]:
fig = plt.figure(figsize = (15,3))
fig = plot_spectra(nuv)
temp = plt.xlim([2530, 2590])
temp = plt.xlim([2650, 2700])
temp = plt.xlim([2790, 2810])
temp = plt.ylim([0,3e-1])
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:
FUV atomic lines. Often with ISM abs -> absolute wavecal possible
What are the other NUV lines, e.g.
In [ ]:
In [ ]:
In [ ]:
In [24]:
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 [25]:
from astropy.stats import sigma_clip
line=cos[9].slice_disp([2660*u.AA, 2700*u.AA])
#line=cos[9].slice_disp([1500*u.AA, 1630*u.AA])
plt.plot(line.disp.value, line.flux.value)
fclip = sigma_clip(line.flux.value, iters=3, sigma=2)
plt.plot(line.disp.value, fclip)
plt.plot(line.disp.value, clip_expand(fclip,n=20))
plt.ylim([0,2e-13])
Out[25]:
In [26]:
cliprange = np.array([[1400,1450], [1450,1500], [1500, 1550], [1580,1630], [1630,1680],
[1680, 1730],
[2560, 2580], [2650, 2700], [2777, 2792]]) * u.Angstrom
In [27]:
cont_flux = np.zeros((10,9))
for i, spec in enumerate(cos):
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))
#cont_flux[i,j] = np.percentile(clip_expand(clipped, n=20), 5)
In [28]:
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$]')
Out[28]:
In [29]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(cliprange.mean(axis=1).value, cont_flux.T-cont_flux.T[:,6][:,None], 'o')
plt.xlabel('central wavelength [$\AA$]')
plt.ylabel('flux [erg/s/cm$^2$]')
Out[29]:
In [30]:
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(cos)):
# subtract the lowest flux as a template of non-accreting state
a,b = curve_fit(fitfunc, cliprange.mean(axis=1).value, cont_flux[i,:] - cont_flux[6,:], p0=[1e-21, 2e5])
fitres[i,:] = a
In [31]:
# Assuming R_* = 0.7
fitres[:, 0] = (fitres[:,0] * (57*u.pc)**2 / (0.7*const.R_sun)**2).decompose()
plt.plot(fitres[:,0])
plt.ylabel('Filling factor of stellar surface')
plt.xlabel('Number of spectrum')
Out[31]:
In [32]:
# Set those to nan where the filling factor = 0
ind = fitres[:,0] < 0.0001
fitres[:,1][ind] = np.nan
plt.plot(fitres[:,1],'o')
plt.ylabel('Black-body temperature [K]')
plt.xlabel('Number of spectrum')
Out[32]:
In [33]:
plt.plot(fitres[:,0], fitres[:,1])
Out[33]:
In [34]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_color_cycle(colorlist)
lines = ax.plot(cliprange.mean(axis=1).value, cont_flux.T-cont_flux.T[:,6][:,None], '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], 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.
The fit is much better is I subtract a base value (e.g. the continuum in Obs 6). There are still some systematics (cp. the middle group of three points in the plot above). Still I find an anticorrelation between area and temperature. I doubt that that's physicla. It might just mean that those parameters are strongly correlated in the fit and I don't learn much above the "flux".
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.
Mal die Bereiche von Kavin France (Christians email) ausprobieren. Das sind Bereiche von dnen wir wissen, dass sie (fast) Linien frei sind. Ich bin mir sicher, dass das nicht viel anders (ich glaube, dass sigma clipping funktioniert), aber kann man ja mal testen.
In [35]:
# This is the line list that Christian gave me for comparison
H2progressions = Table.read('../H2lines/H2_model.dat', format='ascii.commented_header', header_start=3)
In [36]:
H2prorg = H2progressions.group_by('codeu')
In [37]:
H2prorg.groups[1][:-10]
Out[37]:
In [901]:
# Let's see if I can use the Abgrall line list directly. Reduce the risk of errors in the propagation
# of derived files.
# If this is the same, I'll use the routine I already wrote.
data = np.loadtxt(H2.H2.Abgrall93file, skiprows=3, dtype=[('vu','i4'), ('Ju','i4'), ('vl','i4'), ('Jl','i4'),
('A','f4'), ('wavenumber','f8')])
Abgrall93 = Table(data)
Abgrall93.add_column(Column(name='wave', data=1./Abgrall93['wavenumber'] * 1e8, unit=u.AA))
# Add string columns useful for output to LaTeX table etc.
out = []
for d in data:
if d['Ju'] - d['Jl'] == -1:
letter = 'P'
elif d['Ju'] - d['Jl'] == 1:
letter = 'R'
else:
raise NotImplementedError
code = '{0}({1}) {2}-{3}'.format(letter, d['Jl'], d['vu'], d['vl'])
out.append(code)
Abgrall93.add_column(Column(name='speccode', data=out))
out = []
for d in data:
out.append('Ju:{0} vu:{1}'.format(d['Ju'], d['vu']))
Abgrall93.add_column(Column(name='upperlevel', data=out))
Abgrall93 = Abgrall93.group_by('upperlevel')
Abgrall93.add_column(Column(name='branch', dtype=float, length=len(Abgrall93)))
for g in Abgrall93.groups:
g['branch'] = g['A'] / g['A'].sum()
At this point I have two options with different pros and cons:
1) BFG: make Spectra for each progression and fit shift, width and norm
2) Use summer student's work and use fluxes of already identified lines.
Implement 2 for new and see how far I get.
In [466]:
from H2 import sherpacode
import H2
import H2.sherpacode
reload(H2.sherpacode)
from sherpa import models
from sherpa.data import Data1D
from sherpa.astro import models as astromodels
sys.path.append('/melkor/d1/guenther/soft/python/filili/')
from filili.shmodelshelper import copy_pars
# Make one instance of every type of model that we need
# and set sensible defaults for the numbers.
constbase = models.Const1D('baseconst')
constbase.c0 = 5e-14
linebase = astromodels.Lorentz1D('linebase')
linebase.fwhm = 0.07
linebase.fwhm.min = .04
linebase.fwhm.max = 1.
linebase.ampl.max = 5e-12
linebase.ampl = 2e-13
linebase.ampl.min = 0
abslinebase = linebase.__class__('abslinebase')
copy_pars(linebase, abslinebase)
abslinebase.ampl.min = -2e-12
abslinebase.ampl = -2e-13
abslinebase.ampl.max = 0
H2linebase = linebase.__class__('H2linebase')
H2linebase.fwhm.min = .04
H2linebase.fwhm.max = 0.09
H2linebase.ampl.max = 1e-12
H2linebase.ampl = 2e-13
H2linebase.ampl.min = 0
H2linebase.fwhm.val = 0.05681
H2linebase.fwhm.frozen = True
In [467]:
def gabrielcsv2dict(gabtable):
'''Convert Gabriels csv tables with prelininary fits results into nested dicts and lists.
During this time as a summer intern, Gabriel performed preliminary fits of the emission lines
in the COS spectra. In particular, he visually compared the existing STIS data with our COS observations
and used that as guidance to decide where we can split the spectra into small regions with just a handful of
lines which can be fit independently, to improve the speed and the numerical stability of the fits.
He then examined the fits by hand, checking for fit artifacts
'''
gg = gabtable.group_by('Region')
regions = OrderedDict()
for g in gg.groups:
H2linelist = []
linelist = []
abslinelist = []
for l in g:
if l['line'].startswith('H2'):
safename = re.sub('[^0-9a-zA-Z]', '_', l['line'])
H2linelist.append({'name': safename, 'pos.val': l['wave_in']})
elif l['line'].startswith('?ab'):
safename = 'abs' + l['line'].replace('.', '_')
abslinelist.append({'name': safename, 'pos.val': l['wave_in']})
else:
safename = 'em' + l['line'].replace('.', '_')
linelist.append({'name': safename, 'pos.val': l['wave_in']})
fmodellist = [[{'c0.val': g['const1d'][0]}]]
basemodels = [constbase]
if H2linelist:
fmodellist.append(H2linelist)
basemodels.append(H2linebase)
if linelist:
fmodellist.append(linelist)
basemodels.append(linebase)
if abslinelist:
fmodellist.append(abslinelist)
basemodels.append(abslinebase)
regions[g['Region'][0]] = {'range': [g['start'][0], g['stop'][0]],
'fmodellist': fmodellist, 'basemodels': basemodels}
return regions
gabregion = gabrielcsv2dict(gabtabsmaster)
In [469]:
reload(H2)
Abgraldict = H2.H2.read_Abgrall93()
In [490]:
import filili.low_fit
sys.path.append('/melkor/d1/guenther/soft/python/COSlsf/')
from COSlsf import empG160M, tabNUV
class COSModelMaker(filili.low_fit.ModelMaker):
'''Make a model doe out COS spectra.
This adds three features to its base class by overriding finalize_model:
- The wavelength for H2 lines are set from the Abgrall et al. (1993) line list.
- The wavelength of several H2 lines in each region is coupled.
- The entire model is warped into the appropriate LSF fro mthe COSlsf module.
'''
def finalize_model(self, modellist):
model = super(COSModelMaker, self).finalize_model(modellist)
# find which sublist has the H2 lines.
for sublist in modellist:
if "H2" in sublist[0].name:
for l in sublist:
# Set to Abgrall line list + TW Hya stellar velocity
code = l.name.split('_')
wave = Abgraldict['{2}-{3} {0}({1})'.format(code[1], code[2], code[4], code[5])]
l.pos.val = wave # We calibrated the wavelength, so that H2 is at rest.
l.pos.max = wave + 0.05
l.pos.min = wave - 0.05
# freeze difference
filili.low_fit.constant_difference(sublist, 'pos')
# look at some line that has a position to decide if this FUV or NUV
if hasattr(sublist[0], 'pos'):
samplewave = sublist[0].pos.val
if samplewave < 2000:
return empG160M(model)
else:
return tabNUV(model)
In [491]:
from sherpa.data import Data1D
import filili
reload(filili)
import filili.shmodelshelper
reload(filili.shmodelshelper)
import filili.low_fit
reload(filili.low_fit)
modelmaker = COSModelMaker()
fitter = filili.low_fit.Fitter()
from sherpa import optmethods
#fitter.fit_settings['method'] = optmethods.MonCar()
fitter.fit_settings['method'].config['epsfcn'] = 2.2e-16
In [627]:
# While many of the fits work really well without intervention,
# there are cases where the fit "jumps out" of the solution we are
# looking for, e.g. an emission line becomes so wide that it adds to the continuum instead of fitting
# a weak emission line we care for.
# Here, I have some code to interactively modify the fit parameters for a specific region. My goal is to define starting
# conditions for the fit that ensure that each model component has the same meaning for each epoch,
# e.g. prevent lines from shifting into other features to the left or right.
gabmaster = Table.read('../H2lines/master.csv', format='ascii.csv')
gabregion = gabrielcsv2dict(gabmaster)
# Note sure this feature is physical, but the continuum on the left and right side of the line differ.
# Fitting a line here works.
gabregion['6']['fmodellist'][2][0]['pos.frozen'] = True
# region 7: FWHM fluxtuates widely, most likely it's just not well constrained.
# Might profit from fitting all profiles jointly with coupled FWHM.
# However, no H2 in region 7 and lines are unidentified -> low priority for now
gabregion['10']['fmodellist'][2][0]['fwhm.frozen'] = True
# region 11: width of abs feature varies. Feature moves. No good start value for all fits.
# -> refit manually later.
gabregion['12']['fmodellist'][2][0]['fwhm.frozen'] = True
# reg 14 - > revisit, do not use P (18) 3-4
gabregion['14']['fmodellist'][2][0]['fwhm.frozen'] = True
gabregion['14']['fmodellist'][2][1]['fwhm.frozen'] = True
gabregion['16']['fmodellist'][2][0]['pos.val'] = 1429.95
gabregion['16']['fmodellist'][2][1]['pos.val'] = 1430.15
gabregion['16']['fmodellist'][2][0]['pos.max'] = gabregion['16']['fmodellist'][2][0]['pos.val'] + 0.05
gabregion['16']['fmodellist'][2][0]['pos.min'] = gabregion['16']['fmodellist'][2][0]['pos.val'] - 0.05
gabregion['16']['fmodellist'][2][1]['pos.max'] = gabregion['16']['fmodellist'][2][1]['pos.val'] + 0.05
gabregion['16']['fmodellist'][2][1]['pos.min'] = gabregion['16']['fmodellist'][2][1]['pos.val'] - 0.05
gabregion['18']['fmodellist'][1][0]['fwhm.max'] = 0.1
gabregion['18']['fmodellist'][1][1]['fwhm.max'] = 0.1
gabregion['18']['fmodellist'][1][0]['pos.val'] = 1432.6
gabregion['18']['fmodellist'][1][0]['pos.max'] = gabregion['18']['fmodellist'][1][0]['pos.val'] + 0.1
gabregion['18']['fmodellist'][1][0]['pos.min'] = gabregion['18']['fmodellist'][1][0]['pos.val'] - 0.1
# region 20b: Variable, but no H2 -> low priority
gabregion['38']['fmodellist'][2][2]['fwhm.frozen'] = True
# region 46: C I line is crazy. Don't trust P(8) 4-8, but P(11) 1-6 is OK.
# Look at C I line, because we might learn something there.
gabregion['58']['fmodellist'][2][0]['fwhm.max'] = 0.2
gabregion['58']['fmodellist'][2][0]['fwhm.val'] = 0.1
gabregion['58']['fmodellist'][3][0]['fwhm.max'] = 0.1
gabregion['58']['fmodellist'][3][0]['fwhm.val'] = 0.05
gabregion['58']['fmodellist'][3][0]['pos.val'] = 1526.7
gabregion['86']['fmodellist'][1][0]['pos.val'] = 1600.9
gabregion['89']['fmodellist'][2][0]['fwhm.max'] = 0.1
gabregion['89']['fmodellist'][2][0]['fwhm.val'] = 0.05
gabregion['89']['fmodellist'][2][1]['fwhm.max'] = 0.1
gabregion['89']['fmodellist'][2][1]['fwhm.val'] = 0.05
gabregion['89']['fmodellist'][1][0]['pos.max'] = 1608.5
gabregion['89']['fmodellist'][1][0]['pos.min'] = 1608.35
# TODO: NUV regions. Since they don't have H2, they are low priority for now.
# Add lines here that were not part of the Gabriel list, but that France et al (2012) use.
gabregion['france1'] = {'basemodels': [constbase, H2linebase],
'fmodellist': [[{'c0.val': 9.8400024825599998e-14}],
[{'name': 'H2_P_13__2_8', 'pos.val': 1588.79}]],
'range': [1588., 1589.]}
# Not included in our setup
#gabregion['france2'] = {'basemodels': [constbase, H2linebase],
# 'fmodellist': [[{'c0.val': 9.8400024825599998e-14}],
# [{'name': 'H2_P_2__0_2', 'pos.val': 1221.95}]],
# 'range': [1221.5, 1222.5]}
gabregion['france3'] = {'basemodels': [constbase, H2linebase],
'fmodellist': [[{'c0.val': 9.8400024825599998e-14}],
[{'name': 'H2_R_12__4_8', 'pos.val': 1509.45}]],
'range': [1509, 1510.]}
gabregion['france4'] = {'basemodels': [constbase, H2linebase, linebase],
'fmodellist': [[{'c0.val': 9.8400024825599998e-14}],
[{'name': 'H2_P_14__3_9', 'pos.val': 1608.33}],
[{'name': 'em?1608_15', 'pos.val': 1608.15, 'fwhm.frozen': True}]],
'range': [1608., 1609.]}
gabregion['france5'] = {'basemodels': [constbase, H2linebase],
'fmodellist': [[{'c0.val': 9.8400024825599998e-14}],
[{'name': 'H2_P_5__4_9', 'pos.val': 1526.54}]],
'range': [1526., 1527.]}
#gabregion['france6'] = {'basemodels': [constbase, H2linebase],
# 'fmodellist': [[{'c0.val': 9.8400024825599998e-14}],
# [{'name': 'H2_R_14__4_6', 'pos.val': 1415.33}]],
# 'range': [1415., 1416.]}
reg 1: We fit 3-5~R(12) and 1-5~P(5) in this region. The first of those lines appreas consistently shifted towards the blue, possibly due to an unidentified blend.
reg 2: This region contains two $H_2$ lines (0-5~R(0) and 0-5~R(1)) that are blended with the wide \ion{Si}{4} line. The line fluxes are fit well, but both $H_2$ lines are shifted in position. The line shift might be real or an artifact of an asymmetric \ion{Si}{4} line, that we approximate using a Lorentian line profile.
reg 3: Line 0-5~R(2) always apprears slightly blue-shifted, most notable in spectrum 5, 6 and 10.
reg 5: This a clean region of the spectrum with two lines that blend together only in the wings. Following \citet{2002ApJ...572..310H}, we identify these lines as 0-5 P(2) and 2-5 R(11). While 2-5 R(11) is found at the listed wavelength, 0-5~P(2) is shifted about 0.05~\AA{} towards the blue with respect to the theoretical wavelength. \textbf{What does this mean?}
reg 15: 4-7~P(5) is a weak line on the red wing of a broader and stronger unidentified line; thus the fits are not reliable. We do not consider 4-7~P(5) any further.
reg 69: 2-8~R(11) is fit slightly red-shifted. The line is isolated and well fit, but found very close to the edge of the detector, where the wavelengths calibration might not be as accurate any longer.
reg 77: This region contains 4-10~R(12) at 1586.81~\AA. Some spectra show an absorbtion feature next to it around 1567.0~\AA{}, in others so sub-continuum absorption is seen, but the emission lines is too broad for an $H_2$ line. This other, unidentified feature can also be seen in the STIS spectra of \citet{2002ApJ...572..310H}. Since we cannot resolve it in the COS spectra, we do not consider the 4-10~R(12) line any further.
reg 89: The region is dominated by two unidentified features. The 5-12~R(0) line is seen only as an extension of the line wing. Without identifications (and thus fixed line positions) for the other features the 5-12~R(0) line cannot be fit reliably. Thus, we do not consider it any further.
france3: In our data, the 4-8~R(12) line is relatively weak and in several spectra the fit shifts the centroid considerably. The theoretical wavelength appreas in a weak feature, that might be wider than an $H_2$ line or contain a blend; the low signal in this region is insufficient to discern these options. We thus do not consider the 4-8~R(12) line any further.
france5: This is a more complex region consisting of either one broad line with an absorbtion feature in the center or two narrow lines. \citet{2012ApJ...756..171F} identify the left component as 4-9~P(5) and this seems resonable in some of our epochs. Most of the time, however, the observed feature is narrower than the other $H_2$ lines and the line shape indicates that absorption on the right side of the line takes away some of the flux. Thus, we do not consider 4-9~P(5) any further.
france6: \citet{2012ApJ...756..171F} suggest that the 4-6~R(12) line is relatively isolated. In our data, we do not see this line at all.
In [636]:
for r in reps: print r.results['test']['parvals'][1]
In [608]:
# Compare listed with fitted H2 wavelength
print Abgraldict['2-8 P(13)'], Abgraldict['0-2 P(2)'], Abgraldict['4-8 R(12)'], Abgraldict['3-9 P(14)'], Abgraldict['4-9 P(5)'], Abgraldict['4-6 R(12)']
In [671]:
### Code for debugging individual regions in the list interactively
selectregion = OrderedDict(test=gabregion['36'])
print selectregion
#selectregion['test']['fmodellist'][2][2]['fwhm.frozen'] = True
#selectregion['test']['fmodellist'][2][1]['pos.val'] = 1608.3
# Make a separate reporter for each region
reps = [filili.low_fit.SherpaReporter(plot_path='/melkor/d1/guenther/projects/TWHya/plots/{0}_'.format(i))
for i in range(10)]
for i, c in enumerate(cos):
master = filili.low_fit.Master(modelmaker=modelmaker, fitter=fitter, fitreporter=reps[i], confreporter=reps[i])
data = Data1D('cos0', c['WAVELENGTH'], c['FLUX'], c['ERROR'])
master.loop_regions(data, selectregion)
In [628]:
# Make a separate reporter for each COS spectrum
reporters = [filili.low_fit.SherpaReporter(plot_path='/melkor/d1/guenther/projects/TWHya/plots/{0}_'.format(i))
for i in range(10)]
for i, c in enumerate(cos):
master = filili.low_fit.Master(modelmaker=modelmaker, fitter=fitter, fitreporter=reporters[i], confreporter=reporters[i])
data = Data1D('cos0', c['WAVELENGTH'], c['FLUX'], c['ERROR'])
master.loop_regions(data, gabregion)
In [629]:
# make a second directory, where the images are ordered such that the same wavelength are next to each other
im_to_copy = glob.glob('/melkor/d1/guenther/projects/TWHya/plots/*_*.png')
for i in im_to_copy:
filename = os.path.basename(i)
n1, n2 = filename.split('_')
shutil.copy(i, '/melkor/d1/guenther/projects/TWHya/plots2/{0}_{1}.png'.format(n2.split('.')[0], n1))
In [497]:
def get_model_parts(model):
if hasattr(model, 'parts'):
for p in model.parts:
for parts in get_model_parts(p):
yield parts
else:
yield model
list(get_model_parts(master.last_model))
Out[497]:
In [498]:
def tabfromsherpareporter(reporter):
# Get a list of all existing model names in this reporter
modelnames = set(n.split('.')[0] for k in reporter.results.keys() for n in reporter.results[k]['parnames'])
parnames = set(n.split('.')[1] for k in reporter.results.keys() for n in reporter.results[k]['parnames'])
# Turn the dict of stuff into a useful table
pnames = []
pvals = []
pmins = []
pmaxes = []
for n, reg in reporter.results.iteritems():
pnames.extend(reg['parnames'])
pvals.extend(reg['parvals'])
pmins.extend(reg['parmins'])
pmaxes.extend(reg['parmaxes'])
# make astropy table
fittab = astropy.table.Table([list(modelnames)], names=['mdlname'], masked=True)
for c in parnames:
for n in ['val', 'min', 'max']:
fittab.add_column(MaskedColumn(name='{0}.{1}'.format(c, n), dtype=float, length=len(fittab),
mask=True, format='4.2e'))
for i, n in enumerate(pnames):
modname = n.split('.')[0]
pname = n.split('.')[1]
ind = (fittab['mdlname'] == modname).nonzero()[0][0]
for v, c in zip([pvals[i], pmins[i], pmaxes[i]], ['val', 'min', 'max']):
fittab[ind]['{0}.{1}'.format(pname, c)] = v
# Add columns for upper and lower level for H2 lines
for n in ['Ju', 'Jl', 'vu', 'vl']:
fittab.add_column(MaskedColumn(name=n, dtype=np.int, length=len(fittab), mask=True))
for l in fittab:
if l['mdlname'][0:2]=='H2':
Ju, Jl, vu, vl = H2.H2.string2numbers('{0} {1}({2}) {4}-{5}'.format(*l['mdlname'].split('_')))
l['Ju'] = Ju
l['Jl'] = Jl
l['vu'] = vu
l['vl'] = vl
return fittab
In [630]:
fittabs = [tabfromsherpareporter(r) for r in reporters]
In [631]:
fittabs[0]
Out[631]:
In [637]:
# cannot have missing values in merging, so filter to H2 lines only first.
H2tabs = []
for fittab in fittabs:
H2tab = fittab[~fittab['Ju'].mask]
H2tab.remove_columns(['c0.val', 'c0.min', 'c0.max'])
H2tab = astropy.table.join(H2tab, Abgrall93)
H2tab.sort('wave')
H2tab = H2tab.group_by(('Ju','vu'))
H2tabs.append(H2tab)
# remove lines which we know to be unreliable from the fit because of blends.
for n in ['H2_P_8__4_8', 'H2_P_18__3_4', 'H2_P_5__4_7', 'H2_R_12__4_10', 'H2_R_0__5_12', 'H2_P_5__4_9', 'H2_R_12__4_8']:
H2tab.remove_row((H2tab['mdlname'] == n).nonzero()[0][0])
In [638]:
# Values that are high or low indiate that either something wend wrong in the fit or the COS wavelength scale
# needs non-linear corrections.
for H2tab in H2tabs:
plt.plot(H2tab['ampl.val'], H2tab['pos.val'] - H2tab['wave'], '.')
In [639]:
# Find out which lines are not fitted at the right wavelength. Most likely, they don't represent a good fit
# and thus we also should not trust the fitted amplitude. We will inspect those regions again by and and adjust the model
# that is used iteratively.
ind = (np.abs(H2tab['pos.val'] - H2tab['wave']) > 0.04).nonzero()[0]
for i in ind:
for reg in gabregion:
if gabregion[reg]['range'][0] < H2tab['wave'][i] and gabregion[reg]['range'][1] > H2tab['wave'][i]:
print reg, gabregion[reg]['range'], H2tab['mdlname'][i]
In [641]:
# Plot all progressions for one spectrum
H2tab = H2tabs[4]
colors=['r', 'g','b','c','m','y','orange','maroon','0.7', 'pink', 'navy', 'brown', 'salmon','lime', 'r', 'g', 'b', 'c', 'k']
i = 0
H2tab = H2tab.group_by(('Ju','vu'))
for g in H2tab.groups[:].groups:
if len(g) >= 5:
c = colors[i]
i += 1
plotnorm = np.mean(g['A'] * g['wavenumber'] / g['ampl.val'])
# plotnorm=1
plt.plot(g['wave'], g['A'] * g['wavenumber'] / g['ampl.val'] / plotnorm, '-o', color=c)
plt.errorbar(g['wave'],
g['A'] * g['wavenumber'] / g['ampl.val'] / plotnorm,
# Error propagation here! to make that correct
g['A'] * g['wavenumber'] / g['ampl.val']**2 * g['ampl.max'] / plotnorm, color=c)
plt.ylim(0, 2)
#plt.ylim([0, 2e27])
plt.xlabel(r'wavelength [$\AA$]')
plt.ylabel(r'$\frac{A f}{ampl}$ [normalized]')
Out[641]:
In [568]:
# Plot one progression for all spectra and look for changes over time. Pick a good one!
# Check how many lines we have in each progression:
H2tab = H2tabs[0]
H2tab = H2tab.group_by(('Ju','vu'))
for g in H2tab.groups[:].groups:
print 'Progression Ju {0} vu {1} has {2} lines'.format(g['Ju'][0], g['vu'][0], len(g))
In [594]:
Ju = 10
vu = 1
for i, H2tab in enumerate(H2tabs):
g = H2tab[(H2tab['Ju'] == Ju) & (H2tab['vu'] == vu)]
#plotnorm = np.mean(g['A'] * g['wavenumber'] / g['ampl.val'])
plotnorm=1
plt.plot(g['wave'], g['ampl.val'] / (g['A'] * g['wavenumber']) / plotnorm, '-o', color=colorlist[i])
plt.errorbar(g['wave'],
g['ampl.val'] / (g['A'] * g['wavenumber']) / plotnorm,
# Error propagation here! to make that correct
g['ampl.max'] / (g['A'] * g['wavenumber']) / plotnorm, color=colorlist[i])
#plt.ylim(0, 2)
#plt.ylim([0, 5e26])
plt.xlabel(r'wavelength [$\AA$]')
plt.ylabel(r'$\frac{ampl}{A f}$')
Out[594]:
The plot above shows the fitted amplitude of each $H_2$ line in a progression, plotted for each COS spectrum. The energy fluxes are converted to photon fluxes by deviding by the fequency of each transition and by the A value of this transtion. In the absence of blends, fit errors and optical depth, the normalized fluxes would fall on a straight line, where the value of the line is given by the strength of the pumping.
In [605]:
# Compare my line selection with the line selection from France et al. (2012)
francelines = Table.read('../H2lines/France2012tab2.txt', format='ascii', comment='#')
p = re.compile("\((\d+)-(\d+)\)(\w)\((\d+)\)")
mdlname = []
for l in francelines['lineid']:
m = p.match(l)
mdlname.append('H2_{0}_{1}__{2}_{3}'.format(m.group(3), m.group(4), m.group(1), m.group(2)))
francelines['mdlname'] = mdlname
In [642]:
print len(francelines)
H2tab = H2tabs[0]
print len(H2tab)
set1 = set(francelines['mdlname'].data)
set2 = set(H2tab['mdlname'].data)
print set1 - set2
print '----'
print set2 - set1
In [675]:
import mpld3
# Plot all progressions for one spectrum
H2tab = astropy.table.join(H2tabs[1], francelines)
H2tab.sort('wave')
colors=['r', 'g','b','c','m','y','orange','maroon','0.7', 'pink', 'navy', 'brown', 'salmon','lime', 'r', 'g', 'b', 'c', 'k']
i = 0
H2tab = H2tab.group_by(('Ju','vu'))
fig, ax = plt.subplots()
ax.grid(True, alpha=0.3)
for g in H2tab.groups[:].groups:
if len(g) >= 1:
c = colors[i]
i += 1
plotnorm = np.mean(g['A'] * g['wavenumber'] / g['ampl.val'])
# plotnorm=1
line = ax.plot(g['wave'], g['A'] * g['wavenumber'] / g['ampl.val'] / plotnorm, '-o', color=c)
tooltip = mpld3.plugins.PointLabelTooltip(line[0], labels=list(g['mdlname']))
mpld3.plugins.connect(fig, tooltip)
ax.errorbar(g['wave'],
g['A'] * g['wavenumber'] / g['ampl.val'] / plotnorm,
# Error propagation here! to make that correct
g['A'] * g['wavenumber'] / g['ampl.val']**2 * g['ampl.max'] / plotnorm, color=c)
#ax.aylim(0, 2)
#plt.ylim([0, 2e27])
ax.set_xlabel(r'wavelength [Ang]')
ax.set_ylabel(r'Af / ampl')
mpld3.plugins.connect(fig, mpld3.plugins.MousePosition(fontsize=14, fmt='.5g'))
mpld3.display()
Out[675]:
\textbf{For now, we restric ourselves to the line list in \citet{2012ApJ...756..171F}. We can always add more progression later, but those are relatively clean and the branching ratios are already conveniently calcualted in their table.}
Explain how I fut
Explaine how we chose "good" lines
In [680]:
# The following code deomsntrates that Sherpa has a (known) bug, where the "ampl" of a Lorentz function
# really is the "integrated flux".
linetest = astromodels.Lorentz1D('linetest')
linetest.fwhm = 0.02
linetest.ampl = 2e-13
linetest.pos = 1.
x = np.arange(0., 2., 0.01)
y = linetest(x)
print y.sum() * 0.01 # flux in each bin * binwidth
In [690]:
### TODO: Should there be a devision by frequecy to get from erg to photons here?
# Schindhelm does not do that. Maybe it's included in his definiton of the branching ratio?
def progression_from_linelist(H2tab):
'''
Parameters
----------
H2tab : `astropy.table.Table`
Table of measured H2 line fluxes. This table should onlt contain "known good" lines, i.e.
blended lines, lines affected by extra absorbtion etc. need t obe removed before calling this function.
Each line must have the following columns:
- 'Ju' : J of upper level
- 'vu' : v of upper level (those to numbers are used to group lines into progressions)
- 'wave' : wavelength (for sorting)
- 'flux' : dereddened line flux
- 'branch' : the branching ratio for this line
- 'lambda_pump' : pumping wavelength for this transition (will be used for ordering)
TBD : This is a property of the progression, not of the individual line in the input.
It would be better to input this as a table of "Ju, vu, pumping". I can change that later.
For now, I just use the France et al (2012) table as input, which already has this column.
Returns
-------
H2prog : `astropy.table.Table`
Table of progressions, where each progression has an estiamted total flux and an uncertainty.
Typically, the measurement uncertainties of the individual lines are small and the total uncertainties
are dominted by the dereddening or line blends. Thus, the uncertainty is simply evaluated as the standard
deviation of the total progression flux extimated fomr each line.
If there is only one line in a progression, the error field will be masked.
'''
H2tab = H2tab.group_by(('Ju','vu'))
H2prog = Table({'Ju': H2tab.groups.keys['Ju'], 'vu': H2tab.groups.keys['vu']})
H2prog.add_column(MaskedColumn(name='flux', dtype=float, mask=True, length=len(H2prog)))
H2prog.add_column(MaskedColumn(name='e_flux', dtype=float, mask=True, length=len(H2prog)))
H2prog.add_column(MaskedColumn(name='lambda_pump', dtype=float, mask=True, length=len(H2prog)))
for i, group in enumerate(H2tab.groups):
temp = group['flux'] / group['branch']
H2prog['flux'][i] = np.mean(temp)
if len(group) > 1:
H2prog['e_flux'][i] = np.std(temp)
H2prog['lambda_pump'][i] = group['lambda_pump'][0]
H2prog.sort('lambda_pump')
return H2prog
In [755]:
H2tab
Out[755]:
In [785]:
H2progs = []
for t in H2tabs:
t['flux'] = t['ampl.val'] # In Sherpa, ampl.val == integrated flux for Lorentzian lines
t = astropy.table.join(t, francelines)
H2progs.append(progression_from_linelist(t))
In [786]:
H2progs[0]
Out[786]:
Explain the progression from line list. Can refer to Schindhelm 2012 and Herzceg 2004, use there notation and be bried here.
In [775]:
# Read and parse Tab 2 fomr Herczeg et al. 2002 for comparison
herczeg02 = Table.read('../H2lines/Herczeg2002tab2.txt', comment='#', format='ascii.tab')
for n in ['Ju', 'Jl', 'vu', 'vl']:
herczeg02.add_column(Column(name=n, dtype=int, length=len(herczeg02)))
for i, n in enumerate(herczeg02['lineid']):
ju,jl,vu,vl = H2.H2.string2numbers(n)
herczeg02['Ju'][i] = ju
herczeg02['Jl'][i] = jl
herczeg02['vu'][i] = vu
herczeg02['vl'][i] = vl
herczeg02.add_column(Column(name='flux', dtype=float, length=len(herczeg02)))
for i, n in enumerate(herczeg02['data']):
herczeg02['flux'][i] = float(n.split('(')[0]) * 1e-15
# Restrict to France et al 2012 list, because we do the same thing for the COS lines
# Before: lineid in francelist is spelled differently - is does not contain any spaces, but brackets
francenotation = re.compile('\((?P<vu>[0-9]+)-(?P<vl>[0-9]+)\)(?P<dJ>[PR])\((?P<Jl>[0-9]+)\)')
for i, n in enumerate(francelines['lineid']):
m = francenotation.match(n)
francelines['lineid'][i] = '{0}-{1} {2}({3})'.format(m.group('vu'), m.group('vl'), m.group('dJ'), m.group('Jl'))
t = astropy.table.join(herczeg02, francelines)
H2prog_Herczeg02 = progression_from_linelist(t)
In [777]:
H2prog_Herczeg02
Out[777]:
In [926]:
# normalize this to the numbers that Greg observed ?
# tabnorm = H2prog_Herczeg02.copy()
tabnorm = deepcopy(H2progs[0])
tabnorm.rename_column('flux', 'norm')
tabnorm.remove_columns(['e_flux', 'lambda_pump'])
fig = plt.figure()
axl = fig.add_subplot(111)
for i, g in enumerate(H2progs):
# Make an intermediate table to make sure that all columns match up right
g = astropy.table.join(g, tabnorm)
g.sort('lambda_pump')
axl.plot(g['lambda_pump'], g['flux']/g['norm'], '-o', color=colorlist[i])
axl.errorbar(g['lambda_pump'], g['flux']/g['norm'], g['e_flux']/g['norm'], color=colorlist[i])
axl.set_xlabel(r'wavelength [$\AA$]')
axl.set_ylabel('flux [normalized]')
axl.get_xaxis().get_major_formatter().set_useOffset(False)
axl.set_ylim([0.6, 1.3])
fig.savefig(plotdir + 'changeinlya')
The observed fluxes in the different $H_2$ progressions can be used to reconstructed the Ly$\alpha$ profile that is incident on the molecual hydrogen, given some assumptions on the geometry. \citet{2004ApJ...607..369H} present such a calcualtion based on their STIS data of TW~Hya and the Monte-Carlo radiative transfer model of \citet{2002ApJ...575.1057W}. For each progression, the pumping depends on the incident flux at the wavelength of the pumping transition, the opacity of the pumping transition (which in turn depends on the properties of the hydrogen molecule, the geometry of the gas around TW~Hya and the temperature column density of the hydrogen gas). Thus, given a geometry, e.g. a slab of constant temperature, the observations determine the Ly$\alpha$ flux as seen by the hydrogen at different wavelength. Using this method \citet{2004ApJ...607..369H} find that a broad, Gaussian Ly$\alpha$ line, which is absorbed by a single Voigt profile blueshifted 90~km~s$^{-1}$ from line center presents a plausible model for the flux that is incident on the inner disk region where warm ($T=2500$~K), molecular hydrogen is located. They model the absorption feature as a thermal gas with a column density of $\log N(H~I) = 18.7$. More complex scenarios would also include non-thermal level populations in the $H_2$ as expected from X-ray and UV radiation of the central star \citep{2007ApJ...661..334N} and a temperature distribution in the molecular hydrogen.
In the context of variability, we are most interested in how flux and shape of the Ly$\alpha$ line changes over time. Figure~\ref{fig:changeinlya} shows the flux observed in different progressions, ordered by the wavelengths of the pumping transitions. The flux is normalized to the flux observed in the first COS spectrum; as discussed above, the error bars are calculated as the standard deviation of the calcaulted flux in a progression when several lines are available. The figure shows that the flux changes between observations by roughyl 20\% and that essentially all progressions can be described with just a single scaling factor. This indicates that only the total flux of the Ly$\alpha$ line changes in the monitoring period, but not the properties of the molecular hydrogen that is getting excited. If, e.g., the absorbtion between the star and the inner disk edge were to increase significantly, only transitions that are located close to the narrow absorption trough \citep[see Fig.~14 in][]{2004ApJ...607..369H} would be effected, but we see the same scaling at 1214.8~\AA{} as in the other transitions; if the column density of the layer of molecular hydrogen increased significantly, the fluxes in progressions that are optically thin would scale proportionally, while the fluxes in progressions with considerable opacity would not. Figure~\ref{fig:changeinlya} contains transition from levels that are close to the ground state such as 1-0 and others that are only weakly populated such as 16-3. Yet, all progressions shown in that figure are compatible with constant scaling factors.
In [673]:
# find which regions contains a line:
for reg in gabregion:
if gabregion[reg]['range'][0] < 1405 and gabregion[reg]['range'][1] > 1395:
print reg, gabregion[reg]['range']
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [859]:
CIV = np.array([1548.19, 1550.775]) * u.Angstrom
def eq(wave):
return {'equivalencies': u.doppler_optical(wave)}
def plot_doublet(spectra, waves, title='', bin_up=1):
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
for s, c, l in zip(spectra, colorlist, labellist):
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)
ax.set_xlabel(r'v [km/s]')
ax.set_ylabel(r'flux $\left[10^{-12} \frac{\rm{erg}}{\rm{s\;cm}^2 \AA}\right]$')
ax.plot([0,0], ax.get_ylim(), 'k:')
ax.plot(waves[[1,1]].to(u.km/u.second, **eq(waves[0])).value, ax.get_ylim(), 'k:')
ax.set_xlim([-550, 1500])
ax.set_ylim([0, 1.6])
ax.legend()
ax.set_title(title)
return fig
fig = plot_doublet(fuv, CIV, 'C IV doublet')
fig.subplots_adjust(right=0.96, bottom=0.15)
plt.savefig(plotdir + 'CIV')
In [868]:
# The C IV line is blended with several H2 lines, but the only one we can fit separately is 1-8 R(3) at 1547.40
basepoly = sherpa.models.Polynom1D()
selectregion = OrderedDict(test={'basemodels': [basepoly, H2linebase],
'fmodellist': [[{'c0.val': 5e-14, 'c1.frozen': False, 'c2.frozen': False}],
[{'name': 'H2_R_3__1_8', 'pos.val': 1547.40}]],
'range': [1546.8, 1547.8]})
print selectregion
reps = [filili.low_fit.SherpaReporter(plot_path='/melkor/d1/guenther/projects/TWHya/plots/CIVcont_{0}_'.format(i))
for i in range(10)]
for i, c in enumerate(cos):
master = filili.low_fit.Master(modelmaker=modelmaker, fitter=fitter, fitreporter=reps[i], confreporter=reps[i])
data = Data1D('cos0', c['WAVELENGTH'], c['FLUX'], c['ERROR'])
master.loop_regions(data, selectregion)
In [925]:
H2R318fluxes = [r.results['test']['parvals'][4] for r in reps]
H2R318errors = [r.results['test']['parmins'][4] for r in reps]
H2R318 = Abgrall93[Abgrall93['speccode'] == 'R(3) 1-8']
# Check again: Need units of photons / s or erg / s for this calculation?
scale = 1e14
for i, H2prog in enumerate(H2progs):
H2prog41 = H2prog[(H2prog['Ju'] == 4) & (H2prog['vu'] == 1)]
plt.bar(left=i-0.2, width=0.4, bottom=(H2prog41['flux']-H2prog41['e_flux'])*H2R318['branch']*scale,
height=2*H2prog41['e_flux']*H2R318['branch']*scale,
alpha=0.5, linewidth=0, color=colorlist[i],
label='Predicted from other lines in same progression')
plt.errorbar(i, H2R318fluxes[i]*scale, H2R318errors[i]*scale, color=colorlist[i], fmt='o', label='observed')
plt.xlim(-0.25, 9.25)
ticks = plt.xticks(np.arange(10), labellist)
plt.ylabel(r'flux $\left[10^{-14} \frac{\rm{erg}}{\rm{s\;cm}^2}\right]$')
plt.title('$H_2$ 1-8 R(3)')
ax=plt.gca()
handles, labels = ax.get_legend_handles_labels()
plt.legend(handles=handles[:2], numpoints=1, loc='lower left')
plt.ylim(3.7,7.2)
plt.savefig(plotdir + 'H2R318abs')
In [904]:
H2prog41
Out[904]:
In [886]:
H2R318fluxes = [r.results['test']['parvals'][4] for r in reps]
H2R318errors = [r.results['test']['parmins'][4] for r in reps]
H2R318 = Abgrall93[Abgrall93['speccode'] == 'R(3) 1-8']
Ju = 4
vu = 1
for i, H2tab in enumerate(H2tabs):
g = H2tab[(H2tab['Ju'] == Ju) & (H2tab['vu'] == vu)]
#plotnorm = np.mean(g['A'] * g['wavenumber'] / g['ampl.val'])
plotnorm=1
plt.plot(g['wave'], g['ampl.val'] / (g['A'] * g['wavenumber']) / plotnorm, '-', color=colorlist[i])
plt.errorbar(g['wave'],
g['ampl.val'] / (g['A'] * g['wavenumber']) / plotnorm,
# Error propagation here! to make that correct
g['ampl.max'] / (g['A'] * g['wavenumber']) / plotnorm, color=colorlist[i])
plt.errorbar(H2R318['wave']+2*i, H2R318fluxes[i] / (H2R318['A'] * H2R318['wavenumber']) / plotnorm,
H2R318errors[i] / (H2R318['A'] * H2R318['wavenumber']) / plotnorm,
color=colorlist[i])
#plt.ylim(0, 2)
#plt.ylim([0, 5e26])
plt.xlabel(r'wavelength [$\AA$]')
plt.ylabel(r'$\frac{ampl}{A f}$')
Out[886]:
In [870]:
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)
for s, c, l in zip(fuv, colorlist, label):
s = spectrum.spectrum.Spectrum(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]')
ax.set_ylabel(r'flux ratio')
ax.plot([0,0], ax.get_ylim(), 'k:')
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.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')
fig.subplots_adjust(right=0.98, bottom=0.15)
plt.savefig(plotdir + 'CIV_ratio')
In [915]:
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(plotdir + 'MgII')
In [701]:
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
for s, c, l in zip(nuv, colorlist, label):
s = spectrum.spectrum.Spectrum(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)
ax.set_xlabel(r'v [km/s]')
ax.set_ylabel(r'flux ratio'
ax.plot([0,0], ax.get_ylim(), 'k:'
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.grid()
#fig.savefig(plotpath+'optdepth.png')
In [148]:
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')
Out[148]:
In [862]:
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)
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)
ax.set_xlabel(r'v [km/s]')
ax.set_ylabel(r'flux $\left[10^{-12} \frac{\rm{erg}}{\rm{s\;cm}^2 \AA}\right]$')
ax.plot([0,0], ax.get_ylim(), 'k:', lw=2)
ax.set_xlim([-550, 750])
ax.set_ylim([0, 2.6])
ax.legend()
ax.set_title(title)
return fig
fig = plot_line(fuv, HeII, 'He II')
fig.subplots_adjust(right=0.98, bottom=0.15)
plt.savefig(plotdir + 'HeII')
In [150]:
fig = plot_line(fuv, HeII, 'He II', bin_up=4)
fig.axes[0].set_ylim([0, .4])
Out[150]:
In [151]:
fig = plot_line(fuv, [CIV[1]], 'C IV blue', bin_up=2)
fig.axes[0].set_ylim([0, 1])
Out[151]:
In [152]:
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])
Out[152]:
In [153]:
linespec = fuv[0]
line = plt.plot(linespec.disp.to(u.km/u.s, **eq(CIV[0])).value,
linespec.flux.value *1e12*0.5, color = 'k', label=l, lw=2)
line = plt.plot(linespec.disp.to(u.km/u.s, **eq(CIV[1])).value,
linespec.flux.value *1e12, color = 'b', 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])
Out[153]:
In [154]:
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])
Out[154]:
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 [155]:
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])
Out[155]:
In [804]:
#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.
Out[804]:
In [805]:
plt.plot(line.disp.value, line.flux.value)
plt.plot(line.disp.value[~continuum.mask], line.flux.value[~continuum.mask])
Out[805]:
In [806]:
# 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 [818]:
# 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)
# use filled here, because Quantity does not play nicely with masked array
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 [819]:
(line.flux[:-1] * np.diff(line.disp)).sum()
Out[819]:
In [821]:
find_continuum(line, cenfunc=np.ma.median)
Out[821]:
In [824]:
line_flux(line, line, cenfunc=np.ma.median)
Out[824]:
In [822]:
EW(line, line, cenfunc=np.ma.median)
Out[822]:
In [825]:
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 = cos[i].slice_disp(r)
cont = cos[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 [826]:
line=cos[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)
Out[826]:
In [827]:
for col in ewtab.colnames:
plt.plot(ewtab[col]/np.min(ewtab[col]), label=col)
In [828]:
for col in fluxtab.colnames:
plt.plot(fluxtab[col]/np.max(fluxtab[col]), label=col)
In [829]:
temp = plt.plot(cont_flux)
In [830]:
smarts = Table.read('../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(smarts['JD']-2455650, smarts[filt]+i, smarts[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)
fig.savefig(plotdir + 'SMARTSlc.png', bb_inches='tight')
In [831]:
fuvtime = np.array([f.meta['EXPSTRTJ'] for f in fuv])
nuvtime = np.array([f.meta['EXPSTRTJ'] for f in nuv])
In [914]:
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., smarts['B_err']/2.5, label='B (scaled)', elinewidth=2, capthick=2)
ax.errorbar(smarts['JD']-2455650, -smarts['R']/2.5-11.2, smarts['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', 'Mg II doublet']:
ax.plot(fuvtime-2455650, np.log10(fluxtab[col]), 'o-', label=col)
ax.plot(fuvtime-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,70])
ax.set_ylim([-13, -10.75])
plt.savefig(plotdir + 'lc')
In [833]:
for col in ['He II', 'C I multiplet', 'H$_2$ P(3) 0-5', 'Mg II doublet']:
print col, scipy.stats.pearsonr(fluxtab['C IV doublet'], fluxtab[col])
plt.plot(fluxtab['C IV doublet'], fluxtab[col]/np.mean(fluxtab[col]), 'o-', label=col)
plt.plot(fluxtab['C IV doublet'], cont_flux[:,3]/np.mean(cont_flux[:,3]), 'o-', label='FUV cont')
print 'FUV cont', scipy.stats.pearsonr(fluxtab['C IV doublet'], cont_flux[:,3])
bflux = scipy.interpolate.interp1d(smarts['JD'], smarts['B'], kind='nearest')
rflux = scipy.interpolate.interp1d(smarts['JD'], smarts['R'], kind='nearest')
b_interpol = 10**(-bflux(fuvtime)/2.5)
r_interpol = 10**(-rflux(fuvtime)/2.5)
plt.plot(fluxtab['C IV doublet'], b_interpol/np.mean(b_interpol), 'o-', label='B')
print 'B', scipy.stats.pearsonr(fluxtab['C IV doublet'], b_interpol)
plt.plot(fluxtab['C IV doublet'], r_interpol/np.mean(r_interpol), 'o-', label='R')
print 'R', scipy.stats.pearsonr(fluxtab['C IV doublet'], r_interpol)
plt.xlim([1e-12,6e-12])
plt.legend()
Out[833]:
In [834]:
for band in 'BVRIJH':
flux = scipy.interpolate.interp1d(smarts['JD'], smarts[band], kind='nearest')
b_interpol = 10**(-flux(fuvtime)/2.5)
plt.plot(fluxtab['C IV doublet'], b_interpol/np.nanmean(b_interpol), 'o-', label=band)
print band, scipy.stats.pearsonr(fluxtab['C IV doublet'], b_interpol)
plt.legend(loc='upper left')
Out[834]:
Everything is pretty well correlated with each other and with C IV (which I put on the x-axis). B band and FUV continuum scatter most in Fig 1. The second Fig shows that even with the scatter, the correlation is strongest for B and weakest for B-R. The correlations looks stronger in B, but it also has more scatter, so the numbers don't work out (J and H have missing values noted with nan. I could fill those, but that would not make a fair comparison, because then the number of points would vary). So, while not true in any strict numbers, we can see that the correlation is strongest in short wavelenght and falls of the longerwavelengths. Makes sense. Still, in the end everything seems to be correlated with each other (even if I cannot trust the exact numbers.)
Idea: This indicates that all the existing optical monitoring with all the different periods they find really measure periodicty in the accretion (or as we would say from YSOVAR) in the inner disk which drives the accretion.
I don't think that's been written like this so far.
In [710]:
'''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}.
In [ ]:
# DRAFT for an automatic STIS spectrum READER
# This is all nive an good, but for now it's probably better to do cross-correlation with hand selected,
# single orders only.
# So, come back to this later.
hdus = fits.open(datastis+'o59d01030_x1d.fits')
print hdus.info()
data = Table.read(datastis+'o59d01030_x1d.fits', hdu=1)
data['WAVELENGTH'].shape
wavearr = data['WAVELENGTH']
wavelist = [wavearr[i,:] for i in range(wavearr.shape[0]-1,-1,-1)]
mins = np.array([w[0] for w in wavelist])
maxs = np.array(w[-1] for w in wavelist])
ind_overlap = maxs[:-1] > mins[1:]
dispersion = spectrum.coadd.wave_little_interpol(wavelist[ind_overlap])
# take the interpolated dispersion in the beginning and then hstack the later orders.
dispersion = np.hstack([dispersion, wavelist[])
spectralist = [Spectrum({'WAVE': data['WAVELENGTH'][-i,:], 'FLUX': data['FLUX'][-i,:], 'ERROR': data['ERROR'][-i,:]},
dispersion='WAVE', uncertainty='ERROR')
]
STISspec = spectrum.coadd.coadd_errorweighted(spectralist, dispersion)
In [835]:
plt.style.available
Out[835]:
In [836]:
import matplotlib as mpl
mpl.get_configdir()
Out[836]:
In [ ]: