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)