Astropy II: Analyzing UVES Spectroscopy ======================================== This tutorial follows my real live data analysis of MN Lup and the code developed below is taken (with only minor modifications) from the code that I used to actually prepare the publication. The plots that will be developed below appear in very similar form in the article published in ApJ (771, 70, 2013). The examples below depend on each other and the plots in the last section make use of things calculated in the earlier sections. Thus, if you need to restart your python session in the course of this tutorial, please execute all the code again. Also, this tutorial works best if you follow along and execute the code shown on your own computer. The page is already quite long and I do not include the output you would see on your screen in this document. .. toctree:: Before you proceed ------------------ Please download this :download:`this tar file <../downloads/astropy_UVES.tar.gz>` and extract the content, either by clicking on the link or by executing this python code:: from astropy.extern.six.moves.urllib import request import tarfile url = 'http://python4astronomers.github.io/_downloads/astropy_UVES.tar.gz' tarfile.open(fileobj=request.urlopen(url), mode='r|*').extractall() cd UVES ls Then start up IPython with the ``--matplotlib`` option and do the usual imports to enable interactive plotting:: $ ipython --matplotlib import numpy as np import matplotlib.pyplot as plt Acknowledgments --------------- :Authors: Moritz Guenther :Copyright: 2013 Moritz Guenther Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under program ID 087.V0991(A). Scientific background --------------------- In this tutorial we analyze data from MN Lup, a T Tauri star in the Taurus-Auriga star forming region located at a distance of about 140 pc. MN Lup has been observed simultaneously with XMM-Newton and the UVES spectrograph on the VLT. MN Lup is suspected to be a classical T Tauri star, that is accreting mass from a circumstellar disk. MN Lup has been Doppler imaged by `Strassmeier et al. (2005, A&A, 440, 1105) `_ with a very similar UVES setup and those authors claim an rotationally modulated accretion spot. In the X-ray data we find moderate indications for accretion. In this tutorial we analyze (some of) the UVES data to search for rotationally modulated features in the emission line profiles, which could be due to an accretion spot on the stellar surface. Reading the data ---------------- The previous astropy tutorial already covered :ref:`handling-fits-files` and :ref:`wcs-transformations`, so the explanation here is only very brief. Check the `astropy documentation `_ or the other two tutorials for more details:: from glob import glob import numpy as np from astropy.wcs import WCS from astropy.io import fits # glob searches through directories similar to the Unix shell filelist = glob('*.fits') # sort alphabetically - given the way the filenames are # this also sorts in time filelist.sort() Read the first fits file in the list and check what is in there:: sp = fits.open(filelist[0]) sp.info() We see that the data is given as the primary image and all other info is part of the primary header. So, we can extract the WCS from that header to get the wavelength coordinate:: header = sp[0].header wcs = WCS(header) #make index array index = np.arange(header['NAXIS1']) wavelength = wcs.wcs_pix2world(index[:,np.newaxis], 0) wavelength.shape #Ahh, this has the wrong dimension. So we flatten it. wavelength = wavelength.flatten() The flux is contained in the primary image:: flux = sp[0].data Making code reusable as a function ---------------------------------- Now, we do not want to repeat this process for every single file by hand, so let us define a function that takes the filename as input and returns the wavelength and flux arrays and the time of the observation. In python, functions are created with the ``def`` statements. All lines that have an indentation level below the `def` statement are part of the function. Functions can (but do not have to) return values using the ``return`` statement. If a function ``func`` is contained in a file called ``spectra_utils.py`` in the current directory, then this file can be imported into a python session in order to use the function `func` with the following command:: import spectra_utils a = spectra_utils.func(param1, param2, ...) Alternatively, you can import just one (or a few) of many different functions that are defined in your file ``spectra_utils.py``:: from spectra_utils import func a = func(param1, param2, ...) You will recognize that python does not make a difference between modules that come with python (e.g. `glob`), external modules (e.g. `numpy` or `astropy`) and modules that you write yourself. The syntax to import those modules or functions is the same in all cases, provided that the directory where your module is defined is in the search path (`more about python modules and the search path `_). .. note:: You can also define a function on the interactive interpreter by just typing it line by line. However, if you work in an interactive session, then the function ends as soon as there is a blank line. This makes it a little inconvenient to copy and paste more than the simplest functions into an interpreter. Fortunately, ipython offers two magic functions that take care of this: #) You can mark code in some editor, copy it to the clipboard and then type ``%paste`` in ipython. #) You type ``%cpaste`` in ipython you can paste code (e.g. with the mouse) and ipython will correct the leading number of white spaces and ignore empty lines. These two tricks can be used with any type of code, but are particularly useful to define functions. .. note:: Once you used ``import spectra_utils`` python will not monitor the source file. If you change the source code of ``func`` in the file, you need to ``reload(spectra_utils)`` to load the new version of ``func``. So, after all this discussion, we can now define a function that automates the loading of a single spectrum using the commands we developed above. Even if this function is fairly short, we still add some documentation to the header, so that we can look up what parameters it needs when we come back to this project a while later. Personally, I comment every function that is longer than two lines:: def read_spec(filename): '''Read a UVES spectrum from the ESO pipeline Parameters ---------- filename : string name of the fits file with the data Returns ------- wavelength : np.ndarray wavelength (in Ang) flux : np.ndarray flux (in erg/s/cm**2) date_obs : string time of observation ''' sp = fits.open(filename) header = sp[0].header wcs = WCS(header) #make index array index = np.arange(header['NAXIS1']) wavelength = wcs.wcs_pix2world(index[:,np.newaxis], 0) wavelength = wavelength.flatten() flux = sp[0].data date_obs = header['Date-OBS'] return wavelength, flux, date_obs .. admonition:: Exercise Try to find out how you can read the help for this function from the command line. .. raw:: html

Click to Show/Hide Solution

:: help(read_spec) #or read_spec? .. raw:: html
.. admonition:: Exercise The dataset of UVES spectra should have been taken using all the same setup. Write a function that returns the exposure time (``EXPTIME``), the wavelength zero point (``CRVAL1``), and the arm used (UVES has a red and a blue arm - see keyword ``HIERARCH ESO INS PATH``). Then check that all exposures have the same setup. .. raw:: html

Click to Show/Hide Solution

:: def read_setup(filename): '''Get setup for UVES spectrum from the ESO pipeline Parameters ---------- filename : string name of the fits file with the data Returns ------- exposure_time : float wavelength_zero_point : float optical_arm : string ''' sp = fits.open(filename) header = sp[0].header return header['EXPTIME'], header['CRVAL1'], header['HIERARCH ESO INS PATH'] # Let's just print the setup on the screen # We'll see if it's all the same. for f in filelist: print(read_setup(f)) .. raw:: html
The UVES pipeline that was used to reduce the data that we use in the this example employs a fixed wavelength grid (see exercise above), thus the ``wavelength`` is the same for all spectra. This makes it easy to define an array that can hold the fluxes of all observations. Then, we loop over the list of all filenames and fill this array with data:: flux = np.zeros((len(filelist), len(wavelength))) # date comes as string with 23 characters (dtype = 'S23') date = np.zeros((len(filelist)), dtype = 'S23') for i, fname in enumerate(filelist): w, f, date_obs = read_spec(fname) flux[i,:] = f date[i] = date_obs Units and constants in astropy ------------------------------ Often, one has to keep track of the units for certain values. Was the wavelength given in Angstrom or in nm? In X-ray observations, a common unit of wavelength is keV. How many nm is 0.65 keV? `astropy.units `_ offers a framework that can take care of this book-keeping and propagate the units through many (but not all) mathematical operations (e.g. addition, division, multiplication). Furthermore, `astropy.constants `_ supplies the values of many physical and astronomical constants. The easiest way to attach a unit to a number is by multiplication:: import astropy.units as u from astropy.constants.si import c, G, M_sun, R_sun wavelength = wavelength * u.AA # Let's define some constants we need for the exercises further down # Again, we multiply the value with a unit here heliocentric = -23. * u.km/u.s v_rad = -4.77 * u.km / u.s # Strassmeier et al. (2005) R_MN_Lup = 0.9 * R_sun # Strassmeier et al. (2005) M_MN_Lup = 0.6 * M_sun # Strassmeier et al. (2005) vsini = 74.6 * u.km / u.s # Strassmeier et al. (2005) period = 0.439 * u.day # Strassmeier et al. (2005) inclination = 45. * u.degree # Strassmeier et al. (2005) # All numpy trigonometric functions expect the input in radian. # So far, astropy does not know this, so we need to convert the # angle manually incl = inclination.to(u.radian) Now, we can use those variables in our calculations. MN Lup is a T Tauri star (TTS), which is possibly surrounded by an accretion disk. In the spectra we will be looking for signatures of accretion. We expect those accretion signatures to appear close to the free-fall velocity v that a mass m reaches, when it hits the stellar surface. We can calculate the infall speed using simple energy conservation: .. _formula-vaccr: .. math:: E_{kin} = E_{grav} \frac{1}{2} m v^2 = G \frac{m M_*}{R_*} So, let us calculate the free-fall velocity for MN Lup:: >>> v_accr = (2.* G * M_MN_Lup/R_MN_Lup)**0.5 >>> print(v_accr) 504469.027564 m / (s) >>> # Maybe astronomers prefer it in the traditional cgs system? >>> print(v_accr.cgs) 50446902.7564 cm / (s) >>> # Or in some really obscure unit? >>> print(v_accr.to(u.mile / u.hour)) 1128465.07598 mi / (h) How does the accretion velocity relate to the rotational velocity? :: v_rot = vsini / np.sin(incl) v_accr / v_rot Oh, what is that? The seconds are gone, but ``astropy.quantity`` objects keep their different length units unless told otherwise:: (v_accr / v_rot).decompose() The reason for this is that it is not uncommon to use different length units in a single constant, e.g. the Hubble constant is commonly given in "km/ (s * Mpc)". "km" and "Mpc" are both units of length, but generally you do *not* want to shorten this to "1/s". We can now use the ``astropy.units`` mechanism to correct the wavelength scale to the heliocentric velocity scale: .. math:: \lambda_{heliocentric} = \lambda_{bariocentric} * (1 + \frac{v_{helio}}{c}) Naively, we could try:: wavelength = wavelength * (1. + heliocentric/c) TypeError: unsupported operand type(s) for +: 'float' and 'Quantity' However, this fails, because ``heliocentric/c`` is in units of "km/m" and ``1.`` is just a number. From the notation above, it is not clear what we actually want. Do we ask for the value of ``heliocentric/c + 1.`` or do we want to simplify the units of ``heliocentric/c`` and after that add ``1.``? There are several ways to make the instruction precise, but one is to explicitly add ``u.dimensionless_unscaled`` to ``1.`` to tell astropy that this number is dimensionless and does not carry any scaling:: wavelength = wavelength * (1. * u.dimensionless_unscaled+ heliocentric/c) I want to mention one more feature here (check out `astropy.units `_ for more): The ability to convert the spectral axis to frequencies or energies. Normally, a unit of length is not equivalent to a unit of energy or to a frequency, but this conversion makes sense for the wavelength of a spectrum. This is how it can be done:: wavelength.to(u.keV, equivalencies=u.spectral()) wavelength.to(u.Hz, equivalencies=u.spectral()) .. admonition:: Exercise Spectroscopically, MN Lup is classified as spectral type M0 V, thus the gravitational acceleration on the surface :math:`\log(g)` should be comparable to the sun. (For non-stellar astronomers: Conventionally, all values are given in the cgs system. The value for the sun is :math:`\log(g) = 4.4`.) Calculate :math:`\log(g)` for MN Lup with the values for the mass and radius given above. Those values were determined from evolutionary tracks. Check if the :math:`\log(g)` is consistent with the value expected from spectroscopy. .. raw:: html

Click to Show/Hide Solution

The values from evolutionary tracks are indeed consistent with the spectroscopically estimated surface gravity:: >>> print(np.log10((G*M_MN_Lup/R_MN_Lup**2).cgs)) 4.3080943799180433 .. raw:: html
.. admonition:: Exercise Write a function that turns a wavelength scale into a velocity scale. We want to input a wavelengths array and the rest wavelength of a spectral line. We need this function later to show the red- and blueshift of the spectrum relative to the the Ca II H line. Use the following definition to make sure the that code below can use it later:: def wave2doppler(wavelength_array, wavelength_line): .. do something .. return array_of_dopplershifts You can test if you function works by calculating the a Dopplershift of the following wavelengths relative to :math:`H_\alpha`:: waveclosetoHa = np.array([6562.,6563,6565.]) * u.AA print(wave2doppler(waveclosetoHa, 656.489 * u.nm)) I get -132, -86 and +5 km/s. .. raw:: html

Click to Show/Hide Solution

:: def wave2doppler(w, w0): doppler = ((w-w0)/w0 * c) return doppler print(wave2doppler(waveclosetoHa, 656.489 * u.nm).decompose().to(u.km/u.s)) .. raw:: html
.. admonition:: Exercise Write a function that takes a wavelength array and the rest wavelength of a spectral line as input, turns it into a Doppler shift (you can use the function from the last exercise), subtracts the radial velocity of MN Lup (4.77 km/s) and expresses the resulting velocity in units of vsini. We need this function later to show the red- and blueshift of the spectrum relative to the Ca II H line. Use the following definition to make sure the that code below can use it later:: def w2vsini(wavelength_array, wavelength_line): .. do something .. return array_of_shifts_in_vsini .. raw:: html

Click to Show/Hide Solution

:: def w2vsini(w, w0): v = wave2doppler(w, w0) - 4.77 * u.km/u.s return v / vsini .. raw:: html
Converting times ---------------- `astropy.time `_ provides methods to convert times and dates between different systems and formats. Since the ESO fits headers already contain the time of the observation in different systems, we could just read the keyword in the time system we like, but we will use ``astropy.time`` to make this conversion here. ``astropy.time.Time`` will parse many common input formats (strings, floats), but unless the format is unambiguous the format needs to be specified (e.g. a number could mean JD or MJD or year). Also, the time system needs to be given (e.g. UTC). Below are several examples, initialized from different header keywords:: from astropy.time import Time t1 = Time(header['MJD-Obs'], format = 'mjd', scale = 'utc') t2 = Time(header['Date-Obs'], scale = 'utc') Times can be expressed in different formats:: >>> t1