diff --git a/tutorials/UVES/UVES.ipynb b/tutorials/UVES/UVES.ipynb new file mode 100644 index 000000000..745b66821 --- /dev/null +++ b/tutorials/UVES/UVES.ipynb @@ -0,0 +1,1427 @@ +{ + "metadata": { + "astropy-tutorials": { + "author": "Hans Moritz G\u00fcnther , Miguel de Val-Borro ", + "date": "August 2014", + "link_name": "Analyzing UVES Spectroscopy with Astropy", + "name": "", + "published": true, + "description" : "This tutorial demonstrates the data analysis of MN Lup obtained with the UVES spectrograph on the VLT." + }, + "signature": "sha256:1f74796dfdd5b13d06d2a9bd9913a15fe6aa475e7967f8bd969ddcf17fb1e825" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "heading", + "level": 1, + "metadata": {}, + "source": [ + "Analyzing UVES Spectroscopy with Astropy" + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Acknowledgments" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "License: This tutorial is made available under the terms of the\n", + "[Creative Commons Attribution 3.0 License](http://creativecommons.org/licenses/by/3.0/)\n", + "\n", + "Authors: Moritz Guenther, Miguel de Val-Borro\n", + "\n", + "Copyright: 2013, 2014 Moritz Guenther, Miguel de Val-Borro\n", + "\n", + "Based on observations made with ESO Telescopes at the La Silla Paranal Observatory\n", + "under program ID 087.V0991(A)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This tutorial follows my real live data analysis of MN Lup and the code developed\n", + "below is taken (with only minor modifications) from the code that I used to\n", + "actually prepare the publication. The plots that will be developed below\n", + "appear in very similar form in the article published in\n", + "[ApJ, 771, 1, 70](http://adsabs.harvard.edu/abs/2013ApJ...771...70G).\n", + "\n", + "The examples below depend on each other and the plots in the last section make\n", + "use of things calculated in the earlier sections. Thus, if you need to restart\n", + "your python session in the course of this tutorial, please execute all the code\n", + "again.\n", + "\n", + "Also, this tutorial works best if you follow along and execute the code shown\n", + "on your own computer. The page is already quite long and I do not include the\n", + "output you would see on your screen in this document." + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Before you proceed" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Please download the tar file and extract\n", + "the content by executing the python code below.\n", + "Executing this block will download and extract a tar file with data\n", + "necessary for this tutorial. By default it will run in whatever\n", + "directory you have placed this notebook. If that's some place\n", + "you don't want to fill with data files, change the\n", + "location_for_data_file variable." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import tarfile\n", + "from astropy.utils.data import download_file\n", + "url = 'http://data.astropy.org/tutorials/UVES/data_UVES.tar.gz'\n", + "f = tarfile.open(download_file(url, cache=True), mode='r|*')\n", + "location_for_data_file = '.' # CHANGE TO WHEREVER YOU WANT THE DATA TO BE EXTRACTED\n", + "f.extractall(path=location_for_data_file)" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Scientific background" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this tutorial we analyze data from MN Lup, a T Tauri star in the\n", + "Taurus-Auriga star forming region located at a distance of about 140 pc. MN Lup\n", + "has been observed simultaneously with XMM-Newton and the UVES spectrograph on the\n", + "VLT. MN Lup is suspected to be a classical T Tauri star, that is accreting mass\n", + "from a circumstellar disk. MN Lup has been Doppler imaged by\n", + "[Strassmeier et al. 2005](http://adsabs.harvard.edu/abs/2005A%26A...440.1105S)\n", + "with a very similar UVES setup and those authors claim an rotationally modulated\n", + "accretion spot.\n", + "\n", + "In the X-ray data we find moderate indications for accretion. In this\n", + "tutorial we analyze (some of) the UVES data to search for rotationally modulated\n", + "features in the emission line profiles, which could be due to an accretion spot\n", + "on the stellar surface." + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Reading the data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A previous astropy tutorial already covered\n", + "[handling FITS files](../FITS-header.html) and WCS transformations, so the explanation here\n", + "is only very brief. Check the [astropy documentation](http://docs.astropy.org)\n", + "or the other two tutorials for more details:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from glob import glob\n", + "\n", + "import numpy as np\n", + "\n", + "from astropy.wcs import WCS\n", + "from astropy.io import fits\n", + "\n", + "# glob searches through directories similar to the Unix shell\n", + "filelist = glob('UVES/*.fits')\n", + "# sort alphabetically - given the way the filenames are\n", + "# this also sorts in time\n", + "filelist.sort()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Read the first fits file in the list and check what is in there:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "sp = fits.open(filelist[0])\n", + "sp.info()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the data is given as the primary image and all other info is\n", + "part of the primary header. So, we can extract the WCS from that header\n", + "to get the wavelength coordinate.\n", + "If you see warnings about a non-standard RADECSYS, don't worry\n", + "about this - the WCS will still work, it just doesn't\n", + "fully conform to the WCS standard." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "header = sp[0].header\n", + "\n", + "wcs = WCS(header)\n", + "#make index array\n", + "index = np.arange(header['NAXIS1'])\n", + "\n", + "wavelength = wcs.wcs_pix2world(index[:,np.newaxis], 0)\n", + "wavelength.shape\n", + "#Ahh, this has the wrong dimension. So we flatten it.\n", + "wavelength = wavelength.flatten()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The flux is contained in the primary image." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "flux = sp[0].data" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Making code reusable as a function" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we do not want to repeat this process for every single file by hand,\n", + "so let us define a function that takes the filename as input and returns\n", + "the wavelength and flux arrays and the time of the observation.\n", + "In python, functions are created with the ``def`` statements.\n", + "All lines that have an indentation level below the `def` statement are part\n", + "of the function. Functions can (but do not have to) return values using\n", + "the ``return`` statement.\n", + "\n", + "If a function ``func`` is contained in a file called ``spectra_utils.py`` in\n", + "the current directory, then this file can be imported into a python session in\n", + "order to use the function `func` with the following command\n", + "\n", + "```import spectra_utils\n", + "a = spectra_utils.func(param1, param2, ...)```\n", + "\n", + "Alternatively, you can import just one (or a few) of many different functions\n", + "that are defined in your file ``spectra_utils.py``:\n", + "\n", + "```from spectra_utils import func\n", + "a = func(param1, param2, ...)```\n", + "\n", + "You will recognize that python does not make a difference between modules that come\n", + "with python (e.g. `glob`), external modules (e.g. `numpy` or `astropy`) and modules\n", + "that you write yourself. The syntax to import those modules or functions\n", + "is the same in all cases, provided that the directory where your module is\n", + "defined is in the search path [more about python modules and the search path](http://docs.python.org/2/tutorial/modules.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once you used ``import spectra_utils`` python will not monitor the source file.\n", + "If you change the source code of ``func`` in the file, you need to\n", + "``reload(spectra_utils)`` to load the new version of ``func``.\n", + "\n", + "So, after all this discussion, we can now define a function that automates the\n", + "loading of a single spectrum using the commands we developed above. Even if\n", + "this function is fairly short, we still add some documentation to the header,\n", + "so that we can look up what parameters it needs when we come back to this\n", + "project a while later. Personally, I comment every function that is longer\n", + "than two lines." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def read_spec(filename):\n", + " '''Read a UVES spectrum from the ESO pipeline\n", + "\n", + " Parameters\n", + " ----------\n", + " filename : string\n", + " name of the fits file with the data\n", + "\n", + " Returns\n", + " -------\n", + " wavelength : np.ndarray\n", + " wavelength (in Ang)\n", + " flux : np.ndarray\n", + " flux (in erg/s/cm**2)\n", + " date_obs : string\n", + " time of observation\n", + " '''\n", + " sp = fits.open(filename)\n", + " header = sp[0].header\n", + "\n", + " wcs = WCS(header)\n", + "#make index array\n", + " index = np.arange(header['NAXIS1'])\n", + "\n", + " wavelength = wcs.wcs_pix2world(index[:,np.newaxis], 0)\n", + " wavelength = wavelength.flatten()\n", + " flux = sp[0].data\n", + "\n", + " date_obs = header['Date-OBS']\n", + " return wavelength, flux, date_obs" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Exercise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Try to find out how you can read the help for this function from the\n", + "command line." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "help(read_spec)\n", + "# or\n", + "read_spec?\n", + "# In the IPython notebook, the easiest way to see the help for a function is to type read_spec" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Exercise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The dataset of UVES spectra should have been taken using all the same setup.\n", + "Write a function that returns the exposure time (``EXPTIME``),\n", + "the wavelength zero point\n", + "(``CRVAL1``), and the arm used (UVES has a red and a blue arm - see keyword\n", + "``HIERARCH ESO INS PATH``). Then check that all exposures have the same\n", + "setup." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def read_setup(filename):\n", + " '''Get setup for UVES spectrum from the ESO pipeline\n", + "\n", + " Parameters\n", + " ----------\n", + " filename : string\n", + " name of the fits file with the data\n", + "\n", + " Returns\n", + " -------\n", + " exposure_time : float\n", + " wavelength_zero_point : float\n", + " optical_arm : string\n", + " '''\n", + " sp = fits.open(filelist[0])\n", + " header = sp[0].header\n", + "\n", + " return header['EXPTIME'], header['CRVAL1'], header['HIERARCH ESO INS PATH']\n", + "\n", + "# Let's just print the setup on the screen\n", + "# We'll see if it's all the same.\n", + "for f in filelist:\n", + " print(read_setup(f))" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The UVES pipeline that was used to reduce the data that we use in the this example\n", + "employs a fixed wavelength grid (see exercise above),\n", + "thus the ``wavelength`` is the same for all spectra.\n", + "This makes it easy to define an array that can hold the fluxes of all\n", + "observations. Then, we loop over the list of all filenames and fill this array\n", + "with data." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "flux = np.zeros((len(filelist), len(wavelength)))\n", + "# date comes as string with 23 characters (dtype = 'S23')\n", + "date = np.zeros((len(filelist)), dtype = 'S23')\n", + "\n", + "for i, fname in enumerate(filelist):\n", + " w, f, date_obs = read_spec(fname)\n", + " flux[i,:] = f\n", + " date[i] = date_obs" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 1, + "metadata": {}, + "source": [ + "Units and constants in astropy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Often, one has to keep track of the units for certain values. Was the wavelength\n", + "given in Angstrom or in nm? In X-ray observations, a common unit of wavelength is\n", + "keV. How many nm is 0.65 keV?\n", + "[`astropy.units`](http://docs.astropy.org/en/stable/units/index.html)\n", + "offers a framework that can take\n", + "care of this book-keeping and propagate the units through many (but not all)\n", + "mathematical operations (e.g. addition, division, multiplication).\n", + "Furthermore,\n", + "[`astropy.constants`](http://docs.astropy.org/en/stable/constants/index.html) supplies the values of\n", + "many physical and astronomical constants.\n", + "The easiest way to attach a unit to a number is by multiplication." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import astropy.units as u\n", + "from astropy.constants.si import c, G, M_sun, R_sun\n", + "\n", + "wavelength = wavelength * u.AA\n", + "\n", + "# Let's define some constants we need for the exercises further down\n", + "# Again, we multiply the value with a unit here\n", + "heliocentric = -23. * u.km/u.s\n", + "v_rad = -4.77 * u.km / u.s # Strassmeier et al. (2005)\n", + "R_MN_Lup = 0.9 * R_sun # Strassmeier et al. (2005)\n", + "M_MN_Lup = 0.6 * M_sun # Strassmeier et al. (2005)\n", + "vsini = 74.6 * u.km / u.s # Strassmeier et al. (2005)\n", + "period = 0.439 * u.day # Strassmeier et al. (2005)\n", + "\n", + "inclination = 45. * u.degree # Strassmeier et al. (2005)\n", + "# All numpy trigonometric functions expect the input in radian.\n", + "# So far, astropy does not know this, so we need to convert the\n", + "# angle manually\n", + "incl = inclination.to(u.radian)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we can use those variables in our calculations. MN Lup is a T Tauri\n", + "star (TTS), which is possibly surrounded by an accretion disk. In the spectra\n", + "we will be looking for signatures of accretion. We expect those accretion\n", + "signatures to appear close to the free-fall velocity v that a mass m reaches, when\n", + "it hits the stellar surface. We can calculate the infall speed using simple\n", + "energy conservation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$E_{kin} = E_{grav}$$\n", + "$$\\frac{1}{2} m v^2 = G \\frac{m M_*}{R_*}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So, let us calculate the free-fall velocity for MN Lup." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "v_accr = (2.* G * M_MN_Lup/R_MN_Lup)**0.5 \n", + "print(v_accr)\n", + "# Maybe astronomers prefer it in the traditional cgs system?\n", + "print(v_accr.cgs)\n", + "# Or in some really obscure unit?\n", + "from astropy.units import imperial\n", + "print(v_accr.to(imperial.yd / u.hour))" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "How does the accretion velocity relate to the rotational velocity?" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "v_rot = vsini / np.sin(incl)\n", + "v_accr / v_rot" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Oh, what is that? The seconds are gone, but ``astropy.quantity`` objects keep\n", + "their different length units unless told otherwise." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "(v_accr / v_rot).decompose()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The reason for this is that it is not uncommon to use different length units in\n", + "a single constant, e.g. the Hubble constant is commonly given in \"km/ (s * Mpc)\".\n", + "\"km\" and \"Mpc\" are both units of length, but generally you do *not* want to\n", + "shorten this to \"1/s\".\n", + "\n", + "We can now use the ``astropy.units`` mechanism to correct the wavelength\n", + "scale to the heliocentric velocity scale." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$\\lambda_{heliocentric} = \\lambda_{bariocentric} * (1 + \\frac{v_{helio}}{c})$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Naively, we could try:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "wavelength = wavelength * (1. + heliocentric/c)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "However, this fails, because ``heliocentric/c`` is in units of \"km/m\" and ``1.``\n", + "is just a number. From the notation above, it is not clear what we actually want.\n", + "Do we ask for the value of ``heliocentric/c + 1.`` or do we want to simplify the\n", + "units of ``heliocentric/c`` and after that add ``1.``?\n", + "There are several ways to make the instruction precise,\n", + "but one is to explicitly add ``u.dimensionless_unscaled`` to ``1.``\n", + "to tell astropy that this number is dimensionless and does not carry any scaling." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "wavelength = wavelength * (1. * u.dimensionless_unscaled+ heliocentric/c)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I want to mention one more feature here (check out\n", + "[`astropy.units`](http://docs.astropy.org/en/stable/units/index.html) for\n", + "more): The ability to convert the spectral axis to frequencies or energies.\n", + "Normally, a unit of length is not equivalent to a unit of energy or to a\n", + "frequency, but this conversion makes sense for the wavelength of a spectrum.\n", + "This is how it can be done:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "wavelength.to(u.keV, equivalencies=u.spectral())\n", + "wavelength.to(u.Hz, equivalencies=u.spectral())" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Exercise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Spectroscopically, MN Lup is classified as spectral type M0 V, thus\n", + "the gravitational acceleration on the surface $\\log(g)$\n", + "should be comparable to the sun.\n", + "(For non-stellar astronomers: Conventionally, all values are given\n", + "in the cgs system. The value for the sun is $\\log(g) = 4.4$.)\n", + "\n", + "Calculate $\\log(g)$ for MN Lup with the values for the mass\n", + "and radius given above. Those values were determined from\n", + "evolutionary tracks. Check if the $\\log(g)$ is consistent\n", + "with the value expected from spectroscopy." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The values from evolutionary tracks are indeed consistent with the\n", + "spectroscopically estimated surface gravity." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print(np.log10((G*M_MN_Lup/R_MN_Lup**2)/u.cm*u.second**2))" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Exercise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Write a function that turns a wavelength scale into a velocity scale.\n", + "We want to input a wavelengths array and the rest wavelength of a spectral\n", + "line. We need this function later to show the red- and blueshift of the\n", + "spectrum relative to the the Ca II H line. Use the following definition\n", + "to make sure the that code below can use it later.\n", + "You can test if you function works by calculating the a Dopplershift\n", + "of the following wavelengths relative to $H_\\alpha$." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "waveclosetoHa = np.array([6562.,6563,6565.]) * u.AA" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I get -132, -86 and +5 km/s." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# This function uses the Doppler equivalency between wavelength and velocity\n", + "import astropy.units as u\n", + "def wave2doppler(w, w0):\n", + " w0_equiv = u.doppler_optical(w0)\n", + " w_equiv = w.to(u.km/u.s, equivalencies=w0_equiv)\n", + " return w_equiv\n", + "\n", + "print(wave2doppler(waveclosetoHa, 656.489 * u.nm).to(u.km/u.s))" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Exercise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Write a function that takes a wavelength array and the rest wavelength of\n", + "a spectral line as input, turns it into a Doppler shift (you can use\n", + "the function from the last exercise),\n", + "subtracts the radial velocity of MN Lup (4.77 km/s) and expresses\n", + "the resulting velocity in units of vsini.\n", + "We need this function later to show the red- and blueshift of the\n", + "spectrum relative to the Ca II H line. Use the following definition\n", + "to make sure the that code below can use it later." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def w2vsini(wavelength_array, wavelength_line):\n", + " # .. replace this with your implementation ..\n", + " return array_of_shifts_in_vsini" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def w2vsini(w, w0):\n", + " v = wave2doppler(w, w0) - 4.77 * u.km/u.s\n", + " return v / vsini" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Converting times" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[`astropy.time`](http://docs.astropy.org/en/stable/time/index.html)\n", + "provides methods to convert times and dates between different\n", + "systems and formats. Since the ESO fits headers already contain the time of the\n", + "observation in different systems, we could just read the keyword in the time\n", + "system we like, but we will use ``astropy.time`` to make this conversion here.\n", + "``astropy.time.Time`` will parse many common input formats (strings, floats), but\n", + "unless the format is unambiguous the format needs to be specified (e.g. a number\n", + "could mean JD or MJD or year). Also, the time system needs to be given (e.g. UTC).\n", + "Below are several examples, initialized from different header keywords." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from astropy.time import Time\n", + "t1 = Time(header['MJD-Obs'], format = 'mjd', scale = 'utc')\n", + "t2 = Time(header['Date-Obs'], scale = 'utc')" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Times can be expressed in different formats." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "t1\n", + "t1.isot\n", + "t2" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "or be converted to a different time system." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "t1.tt" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Times can also be initialized from arrays and we can calculate time differences." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "obs_times = Time(date, scale = 'utc')\n", + "delta_t = obs_times - Time(date[0], scale = 'utc')" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we want to express the time difference between the individual spectra of\n", + "MN Lup in rotational periods. While the unit of ``delta_t`` is days, unfortunately\n", + "``astropy.time.Time`` and ``astropy.units.Quantity`` objects do not work together\n", + "yet, so we will have to convert from one to the other explicitly." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "delta_p = delta_t.value * u.day / period" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Normalize the flux to the local continuum" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we want to look at the time evolution of a single specific\n", + "emission line in the spectrum. In order to estimate the equivalent width\n", + "or make reasonable plots we need to normalize the flux to the local continuum.\n", + "In this specific case the emission line is bright and the continuum can be\n", + "described reasonably by a second-order polynomial." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So, we define two regions left and right of the emission line, where we fit the\n", + "polynomial. Looking at the figure, ``[3925*u.AA, 3930*u.AA]`` and\n", + "``[3938*u.AA, 3945*u.AA]`` seem right for that. Then, we normalize the flux by\n", + "this polynomial.\n", + "\n", + "The following function will do that:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def region_around_line(w, flux, cont):\n", + " '''cut out and normalize flux around a line\n", + "\n", + " Parameters\n", + " ----------\n", + " w : 1 dim np.ndarray\n", + " array of wavelengths\n", + " flux : np.ndarray of shape (N, len(w))\n", + " array of flux values for different spectra in the series\n", + " cont : list of lists\n", + " wavelengths for continuum normalization [[low1,up1],[low2, up2]]\n", + " that described two areas on both sides of the line\n", + " '''\n", + "#index is true in the region where we fit the polynomial\n", + " indcont = ((w > cont[0][0]) & (w < cont[0][1])) |((w > cont[1][0]) & (w < cont[1][1]))\n", + "#index of the region we want to return\n", + " indrange = (w > cont[0][0]) & (w < cont[1][1])\n", + "# make a flux array of shape\n", + "# (number of spectra, number of points in indrange)\n", + " f = np.zeros((flux.shape[0], indrange.sum()))\n", + " for i in range(flux.shape[0]):\n", + "# fit polynomial of second order to the continuum region\n", + " linecoeff = np.polyfit(w[indcont], flux[i, indcont],2)\n", + "# divide the flux by the polynomial and put the result in our\n", + "# new flux array\n", + " f[i,:] = flux[i,indrange]/np.polyval(linecoeff, w[indrange])\n", + " return w[indrange], f\n", + "\n", + "wcaII, fcaII = region_around_line(wavelength, flux,\n", + " [[3925*u.AA, 3930*u.AA],[3938*u.AA, 3945*u.AA]])" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Publication ready output" + ] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Tables" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will calculate the equivalent width in Angstroms of the emission line\n", + "for the first spectrum." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "ew = fcaII[0,:] - 1.\n", + "ew = ew[:-1] * np.diff(wcaII.to(u.AA).value)\n", + "print(ew.sum())" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using ``numpy`` array notation we can actually process all spectra at once." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "delta_lam = np.diff(wcaII.to(u.AA).value)\n", + "ew = np.sum((fcaII - 1.)[:,:-1] * delta_lam[np.newaxis, :], axis=1)" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, we want to generate a LaTeX table of the observation times, period\n", + "and equivalent width that we can directly paste into our manuscript. To do so,\n", + "we first collect all the columns and make an ``astropy.table.Table`` object. (Please\n", + "check [`astropy.table`](http://docs.astropy.org/en/stable/table/index.html)\n", + "or `tabular-data` for more\n", + "details on ``Table``). So, here is the code:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from astropy.table import Column, Table\n", + "import astropy.io.ascii as ascii\n", + "\n", + "datecol = Column(name = 'Obs Date', data = date)\n", + "pcol = Column(name = 'phase', data = delta_p, format = '{:.1f}')\n", + "ewcol = Column(name = 'EW', data = ew, format = '{:.1f}', unit = '\\\\AA')\n", + "tab = Table((datecol, pcol, ewcol))\n", + "# latexdicts['AA'] contains the style specifics for A&A (\\hline etc.)\n", + "tab.write('EWtab.tex', latexdict = ascii.latexdicts['AA'])" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will make two plots. The plotting is done with\n", + "[`matplotlib`](http://matplotlib.org), and does not involve Astropy itself.\n", + "Plotting is introduced in `plotting-and-images` and more details on\n", + "plotting can be found there. When in doubt, use the search engine of your choice\n", + "and ask the internet. Here, I mainly want to illustrate that Astropy can be\n", + "used in real-live data analysis.\n", + "Thus, I do not explain every step in the plotting in detail.\n", + "The plots we produce below appear in very\n", + "similar form in Guenther et al. 2013 (ApJ, 771, 70)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In both cases we want the x-axis to show the Doppler shift expressed in units\n", + "of the rotational velocity. In this way, features that are rotationally\n", + "modulated will stick out between -1 and +1." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "x = w2vsini(wcaII, 393.366 * u.nm).decompose()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we will show the line profile." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import matplotlib.pyplot as plt\n", + "# set reasonable figsize for 1-column figures\n", + "fig = plt.figure(figsize = (4,3))\n", + "ax = fig.add_subplot(1,1,1)\n", + "ax.plot(x, fcaII[0,:])\n", + "ax.set_xlim([-3,+3])\n", + "ax.set_xlabel('line shift [v sin(i)]')\n", + "ax.set_ylabel('flux')\n", + "ax.set_title('Ca II H line in MN Lup')\n", + "# when using this interface, we need to explicitly call the draw routine\n", + "plt.draw()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Exercise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The plot above shows only a single spectrum. Plot all spectra into a single\n", + "plot and introduce a sensible offset between them, so that we can follow\n", + "the time evolution of the line." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are clearly several ways to produce a well looking plot. Here is one\n", + "way." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "yshift = np.arange((fcaII.shape[0])) * 0.5\n", + "#shift the second night up by a little more\n", + "yshift[:] += 1.5\n", + "yshift[13:] += 1\n", + "\n", + "fig = plt.figure(figsize = (4,3))\n", + "ax = fig.add_subplot(1,1,1)\n", + "\n", + "for i in range(25):\n", + "ax.plot(x, fcaII[i,:]+yshift[i], 'k')\n", + "\n", + "#separately show the mean line profile in a different color\n", + "ax.plot(x, np.mean(fcaII, axis =0), 'b')\n", + "ax.set_xlim([-2.5,+2.5])\n", + "ax.set_xlabel('line shift [$v \\\\sin i$]')\n", + "ax.set_ylabel('flux')\n", + "ax.set_title('Ca II H line in MN Lup')\n", + "fig.subplots_adjust(bottom = 0.15)\n", + "plt.draw()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we will make a more advanced plot. For each spectrum we calculate\n", + "the difference to the mean flux." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fmean = np.mean(fcaII, axis=0)\n", + "fdiff = fcaII - fmean[np.newaxis,:]" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the following simple plot, we can already see features moving through the line.\n", + "However, the axis scales are not right, the gap between both nights is not visible\n", + "and there is no proper labeling." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "fig = plt.figure(figsize = (4,3))\n", + "ax = fig.add_subplot(1,1,1)\n", + "im = ax.imshow(fdiff, aspect = \"auto\", origin = 'lower')" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the following, we will plot the spectra from both nights separately.\n", + "Also, we will pass the ``extent`` keyword to ``ax.imshow`` which takes care\n", + "of the axis." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "ind1 = delta_p < 1 * u.dimensionless_unscaled\n", + "ind2 = delta_p > 1 * u.dimensionless_unscaled\n", + "\n", + "fig = plt.figure(figsize = (4,3))\n", + "ax = fig.add_subplot(1,1,1)\n", + "\n", + "for ind in [ind1, ind2]:\n", + "im = ax.imshow(fdiff[ind,:], extent = (np.min(x), np.max(x), np.min(delta_p[ind]), np.max(delta_p[ind])), aspect = \"auto\", origin = 'lower')\n", + "\n", + "ax.set_ylim([np.min(delta_p), np.max(delta_p)])\n", + "ax.set_xlim([-1.9,1.9])\n", + "plt.draw()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, this plot is already much better, but there are still some things that can be\n", + "improved:\n", + "\n", + "* Introduce an offset on the y-axis to reduce the amount of white space.\n", + "* Strictly speaking, the image shown is not quite the right scale because the\n", + "``extent`` keyword gives the edges of the image shown, while ``x`` and\n", + "``delta_p`` contain the bin mid-points.\n", + "* Use a gray scale instead of color to save publication charges.\n", + "* Add labels to the axis.\n", + "\n", + "The following code addresses these points." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# shift a little for plotting purposes\n", + "pplot = delta_p.copy().value\n", + "pplot[ind2] -= 1.5\n", + "# image goes from x1 to x2, but really x1 should be middle of first pixel\n", + "delta_t = np.median(np.diff(delta_p))/2.\n", + "delta_x = np.median(np.diff(x))/2.\n", + "# imshow does the normalization for plotting really well, but here I do it\n", + "# by hand to ensure it goes -1,+1 (that makes color bar look good)\n", + "fdiff = fdiff / np.max(np.abs(fdiff))\n", + "\n", + "fig = plt.figure(figsize = (4,3))\n", + "ax = fig.add_subplot(1,1,1)\n", + "\n", + "for ind in [ind1, ind2]:\n", + "im = ax.imshow(fdiff[ind,:],\n", + "extent = (np.min(x)-delta_x, np.max(x)+delta_x,\n", + "np.min(pplot[ind])-delta_t, np.max(pplot[ind])+delta_t),\n", + "aspect = \"auto\", origin = 'lower', cmap = plt.cm.Greys_r)\n", + "\n", + "ax.set_ylim([np.min(pplot)-delta_t, np.max(pplot)+delta_t])\n", + "ax.set_xlim([-1.9,1.9])\n", + "ax.set_xlabel('vel in $v\\\\sin i$')\n", + "ax.xaxis.set_major_locator(plt.MaxNLocator(4))\n", + "\n", + "def pplot(y, pos):\n", + "'The two args are the value and tick position'\n", + "'Function to make tick labels look good.'\n", + "if y < 0.5:\n", + "yreal = y\n", + "else:\n", + "yreal = y + 1.5\n", + "return yreal\n", + "\n", + "formatter = plt.FuncFormatter(pplot)\n", + "ax.yaxis.set_major_formatter(formatter)\n", + "ax.set_ylabel('period')\n", + "fig.subplots_adjust(left = 0.15, bottom = 0.15, right = 0.99, top = 0.99)\n", + "plt.draw()" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Exercise" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Understand the code for the last plot. Some of the commands used are\n", + "already pretty advanced stuff. Remember, any internet search engine can be\n", + "your friend." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Clearly, I did not develop this code for scratch.\n", + "The [matplotlib gallery](http://matplotlib.org/gallery.html) is my\n", + "preferred place to look for plotting solutions." + ] + }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Contributing to Astropy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Astropy](http://astropy.org) is an open-source and community-developed\n", + "Python package, which means that is only as good as the contribution of the\n", + "astronomical community. Clearly, there will always people who have more fun writing\n", + "code and others who have more fun using it. However, if you find a bug and do not\n", + "report it, then it is unlikely to be fixed. If you wish for a specific feature,\n", + "then you can either implement it and contribute it or at least fill in a feature\n", + "request.\n", + "\n", + "If you want to get help or discuss issues with other Astropy users, you can\n", + "sign up for the [astropy mailing list](http://mail.scipy.org/mailman/listinfo/astropy).\n", + "Alternatively, the [astropy-dev](http://groups.google.com/group/astropy-dev) list is where you should go to\n", + "discuss more technical aspects of Astropy with the developers.\n", + "\n", + "If you have come across something that you believe is a bug, please open a\n", + "ticket in the Astropy [issue tracker](http://github.com/astropy/astropy/issues), and we will look into it\n", + "promptly.\n", + "\n", + "Please try to include an example that demonstrates the issue and will allow the\n", + "developers to reproduce and fix the problem. If you are seeing a crash\n", + "then frequently it will help to include the full Python stack trace as well as\n", + "information about your operating system (e.g. MacOSX version or Linux version).\n", + "\n", + "Here is a practical example.\n", + "Above we calculated the free-fall\n", + "velocity onto MN Lup like this." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "v_accr = (2.*G *M_MN_Lup / R_MN_Lup)**0.5" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Mathematically, the following statement is equivalent:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "v_accr2 = np.sqrt((2.*G * M_MN_Lup/R_MN_Lup))" + ], + "language": "python", + "metadata": {}, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "However, this raises a warning and automatically converts the quantity object\n", + "to an numpy array and the unit is lost. If you believe that is a bug that should\n", + "be fixed, you might chose to report it in the [issue tracker](http://github.com/astropy/astropy/issues).\n", + "(But please check if somebody else has reported the same thing before, so we do\n", + "not clutter the issue tracker needlessly.)" + ] + } + ], + "metadata": {} + } + ] +} diff --git a/tutorials/UVES/UVES.py b/tutorials/UVES/UVES.py deleted file mode 100644 index d1ef88f17..000000000 --- a/tutorials/UVES/UVES.py +++ /dev/null @@ -1,242 +0,0 @@ -from glob import glob - -import numpy as np -import matplotlib.pyplot as plt - -from astropy.wcs import WCS -from astropy.io import fits - -filelist = glob('data/*.fits') -filelist.sort() -sp = fits.open(filelist[0]) -header = sp[0].header - -wcs = WCS(header) -index = np.arange(header['NAXIS1']) -wavelength = wcs.wcs_pix2world(index[:,np.newaxis], 0) -wavelength = wavelength.flatten() - - -flux = sp[0].data - - - -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 - -flux = np.zeros((len(filelist), len(wavelength))) -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 - - -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 excersices 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) -incl = inclination.to(u.radian) - -wavelength = wavelength * (1. * u.dimensionless_unscaled+ heliocentric/c) - - -def wave2doppler(w, w0): - doppler = ((w-w0)/w0 * c) - return doppler - - - -def w2vsini(w, w0): - v = wave2doppler(w, w0) - 4.77 * u.km/u.s - return v / vsini - -from astropy.time import Time -obs_times = Time(date, scale = 'utc') -delta_t = obs_times - Time(date[0], scale = 'utc') -delta_p = delta_t.val * u.day / period - - - -fig = plt.figure() -ax = fig.add_subplot(1,1,1) -ax.plot(wavelength, flux[0,:]) -ax.set_xlim([3900, 4000]) -ax.set_ylim([0,2000]) -ax.set_xlabel('wavelength [$\AA$]') -ax.set_ylabel('flux') -fig.savefig('CaII.png') - -def region_around_line(w, flux, cont): - '''cut out and normalize flux around a line - - Parameters - ---------- - w : 1 dim np.ndarray - array of wanvelenghts - flux : np.ndarray of shape (N, len(w)) - array of flux values for different spectra in the series - cont : list of lists - wavelengths for continuum normalization [[low1,up1],[low2, up2]] - that described two areas on both sides of the line - ''' - #index is true in the region where we fit the polynomial - indcont = ((w > cont[0][0]) & (w < cont[0][1])) | ((w > cont[1][0]) & (w < cont[1][1])) - #index of the region we want to return - indrange = (w > cont[0][0]) & (w < cont[1][1]) - # make a flux array of shape - # (nuber of spectra, number of pointsin indrange) - f = np.zeros((flux.shape[0], indrange.sum())) - for i in range(flux.shape[0]): - # fit polynom of second order to the continuum region - linecoeff = np.polyfit(w[indcont], flux[i, indcont],2) - # devide the flux by the polynom and put the result in our - # new flux array - f[i,:] = flux[i,indrange]/np.polyval(linecoeff, w[indrange]) - return w[indrange], f - -wcaII, fcaII = region_around_line(wavelength, flux, - [[3925*u.AA, 3930*u.AA],[3938*u.AA, 3945*u.AA]]) - - -x = w2vsini(wcaII, 393.366 * u.nm).decompose() - - -import matplotlib.pyplot as plt -# set reasonable figsize for 1-column figures -fig = plt.figure(figsize = (4,3)) -ax = fig.add_subplot(1,1,1) -ax.plot(x, fcaII[0,:]) -ax.set_xlim([-3,+3]) -ax.set_xlabel('line shift [v sin(i)]') -ax.set_ylabel('flux') -ax.set_title('Ca II H line in MN Lup') -# when using this interface, we need to explicitly call the draw routine -plt.draw() -fig.savefig('CaII-lines-one.png') - - -yshift = np.arange((fcaII.shape[0])) * 0.5 -yshift[:] += 1.5 -yshift[13:] += 1 - -fig = plt.figure(figsize = (4,3)) -ax = fig.add_subplot(1,1,1) - -for i in range(25): - ax.plot(x, fcaII[i,:]+yshift[i], 'k') - -ax.plot(x, np.mean(fcaII, axis =0), 'b') -ax.set_xlim([-2.5,+2.5]) -ax.set_xlabel('line shift [$v \\sin i$]') -ax.set_ylabel('flux') -ax.set_title('Ca II H line in MN Lup') -fig.subplots_adjust(bottom = 0.15) -plt.draw() -fig.savefig('CaII-lines-all.png') - - -fmean = np.mean(fcaII, axis=0) -fdiff = fcaII - fmean[np.newaxis,:] - - -fig = plt.figure(figsize = (4,3)) -ax = fig.add_subplot(1,1,1) -im = ax.imshow(fdiff, aspect = "auto", origin = 'lower') -fig.savefig('CaII-1.png') - -ind1 = delta_p < 1 * u.dimensionless_unscaled -ind2 = ~ind1 - -fig = plt.figure(figsize = (4,3)) -ax = fig.add_subplot(1,1,1) - -for ind in [ind1, ind2]: - im = ax.imshow(fdiff[ind,:], extent = (np.min(x), np.max(x), np.min(delta_p[ind]), np.max(delta_p[ind])), aspect = "auto", origin = 'lower') - -ax.set_ylim([np.min(delta_p), np.max(delta_p)]) -ax.set_xlim([-1.9,1.9]) -plt.draw() -fig.savefig('CaII-2.png') - - -# shift a little for plotting purposes -pplot = delta_p.copy().value -pplot[ind2] -= 1.5 -# image goes from x1 to x2, but really x1 should be middle of first pixel -delta_t = np.median(np.diff(delta_p))/2. -delta_x = np.median(np.diff(x))/2. -# imshow does the normalization for plotting really well, but here I do it -# by hand to ensure it goes -1,+1 (that makes color bar look good) -fdiff = fdiff / np.max(np.abs(fdiff)) -fig = plt.figure(figsize = (4,3)) -ax = fig.add_subplot(1,1,1) - -for ind in [ind1, ind2]: - im = ax.imshow(fdiff[ind,:], - extent = (np.min(x)-delta_x, np.max(x)+delta_x, - np.min(pplot[ind])-delta_t, np.max(pplot[ind])+delta_t), - aspect = "auto", origin = 'lower', cmap = plt.cm.Greys_r) - -ax.set_ylim([np.min(pplot)-delta_t, np.max(pplot)+delta_t]) -ax.set_xlim([-1.9,1.9]) -ax.set_xlabel('vel in $v\\sin i$') -ax.xaxis.set_major_locator(plt.MaxNLocator(4)) - -def pplot(y, pos): - 'The two args are the value and tick position' - 'Function to make tick labels look good.' - if y < 0.5: - yreal = y - else: - yreal = y + 1.5 - return yreal - - -formatter = plt.FuncFormatter(pplot) -ax.yaxis.set_major_formatter(formatter) -ax.set_ylabel('period') -fig.subplots_adjust(left = 0.15, bottom = 0.15, right = 0.99, top = 0.99) -plt.draw() -fig.savefig('CaII-3.png') - diff --git a/tutorials/UVES/UVES.rst b/tutorials/UVES/UVES.rst deleted file mode 100644 index f8b6af499..000000000 --- a/tutorials/UVES/UVES.rst +++ /dev/null @@ -1,841 +0,0 @@ -Analyzing UVES Spectroscopy with Astropy -======================================= - -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. - -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:: - - import urllib2, tarfile - url = 'http://python4astronomers.github.io/_downloads/astropy_UVES.tar.gz' - tarfile.open(fileobj=urllib2.urlopen(url), mode='r|*').extractall() - cd UVES - ls - -Then start up IPython with the ``--pylab`` option to enable easy plotting:: - - ipython --pylab - -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(filelist[0]) - 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 -