In this notebook we're going to be exploring our galaxy by looking at stars around the Sun.
\n", "\n", "We are going to do this using data from two surveys, the Sloan Digital Sky Survey (SDSS) and Gaia. SDSS measured the spectra of millions of stars. The SDSS data reduction pipeline took these spectra and measured some important quantities about each individual star. This includes information about the composition of each star as well as the star's velocity towards or away from the Earth. Gaia measured the positions and motions of billions of stars, giving us lots of information about where an individual star is in the galaxy and how it is moving.
\n", "\n", "## Imports\n", "This tutorial requires three commonly used python packages\n", "- _numpy_ for efficienct math with arrays\n", "- _matplotlib_ for plotting\n", "- _astropy_ for reading fits table\n", "\n", "Missing packages can be installed via `pip install numpy matplotlib astropy` \n", "\n", "Let's start by importing some libraries. The warnings import is entirely optional, but suppresses some annoying formatting errors that may appear on some computers.
" ] }, { "cell_type": "code", "execution_count": null, "id": "4afa1463", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import astropy.coordinates as coord\n", "import astropy.units as u\n", "from astropy.io import fits\n", "from astropy.coordinates import SkyCoord\n", "from astropy.table import Table\n", "\n", "# for plots\n", "from matplotlib import rc\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "markdown", "id": "4d969d3e", "metadata": {}, "source": [ "### Loading data from a fits file\n", "\n", "We can now import our data file. This combines the parameters measured by SDSS and Gaia for lotes of different stars. It loads the data from our input file into a table.
\n", "\n", "**Where**: you can download the data --> https://dr19.sdss.org/sas/dr19/spectro/astra/0.6.0/summary/astraAllStarASPCAP-0.6.0.fits.gz\n", "\n", "**OR** if you are running this notebook on sciServer, the summary file is available at /home/idies/workspace/sdss_sas/dr19/spectro/astra/0.6.0/summary/\n", "\n", "Below, the default loadpath is set to sciServer. If you download it locally, change loadpath as needed, e.g., \"/Users/myuser/Downloads/\"\n", "\n", "Next, we'll use `astropy.fits` to load the file and choose the extension we want. In this case, BOSS data is in extension 1, but we want APOGEE data in extension 2." ] }, { "cell_type": "code", "execution_count": null, "id": "031d8a61", "metadata": {}, "outputs": [], "source": [ "filename='astraAllStarASPCAP-0.6.0.fits.gz'\n", "#tb = fits.open(filename)\n", "data = Table.read(filename, format='fits',hdu=2)" ] }, { "cell_type": "markdown", "id": "55cdfadd", "metadata": {}, "source": [ "Now lets have a look at the data columns
" ] }, { "cell_type": "code", "execution_count": null, "id": "e95a7cd5", "metadata": { "scrolled": true }, "outputs": [], "source": [ "data.colnames" ] }, { "cell_type": "markdown", "id": "d8348af4", "metadata": {}, "source": [ "That's a lot of columns. Luckily we will only be dealing with a few of them.
" ] }, { "cell_type": "markdown", "id": "d2c3f08e", "metadata": {}, "source": [ "**Exercise**, determine how many stars your table has data for." ] }, { "cell_type": "code", "execution_count": null, "id": "c8410ddd", "metadata": {}, "outputs": [], "source": [ "#enter your code here\n" ] }, { "cell_type": "markdown", "id": "0143f6ec", "metadata": {}, "source": [ "Now let's make a map of our data on the sky.
" ] }, { "cell_type": "code", "execution_count": null, "id": "5ca218ec", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(15,7.5))\n", "plt.scatter(data['ra'],data['dec'], color='black', marker='.', s=0.1,label=\"All Stars\" )\n", "plt.xlabel('ra', fontsize=15); plt.ylabel('dec', fontsize=15)\n", "plt.title('Sky plot', fontsize=15)\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "ac07c86a", "metadata": {}, "source": [ "What can you say about the distribution on the sky of the stars in your sample?
\n", "\n", "Gaia covered the whole sky but SDSS only has data from certain parts of the sky. This is a trade-off astronomers often face, cover as wide an area as possible or get as much detailed information on specific regions of the sky. This is a constant balancing act for astronomers designing surveys.
\n", "\n", "One thing you might notice is that near the celestial equator (declination=0) the stars appear grouped together in circular blobs. Each of these is one field observed by SDSS. Near the north and south celestial poles these blobs appear stretched in the ra direction. This isn't because the SDSS observations change near the poles, but because of a projection effect caused by making this plot. It's a bit like how Antarctica appears much larger than it really is on most flat world maps.
" ] }, { "cell_type": "markdown", "id": "479f5c3a", "metadata": {}, "source": [ "Let's start by doing a quick data quality cut on our sample.
" ] }, { "cell_type": "code", "execution_count": null, "id": "7297fee0", "metadata": {}, "outputs": [], "source": [ "mask_sp =(data['teff']>3e3)&(data['teff']<7e3) &(data['logg']<6)&(data['logg']>-1)&(data['snr']>50)\\\n", " &(np.abs(data['mg_h'])<2)&(np.abs(data['fe_h'])<2)&(data['plx'] > 0) & (data['plx']>10.0*data['e_plx'])" ] }, { "cell_type": "markdown", "id": "69f72af1", "metadata": {}, "source": [ "The bit of code above selects a subsample of our data. It selects only objects that meet certain criteria.
\n", "The effective temperature and the $\\log g$ in the star's atmosphere change the depth and shape of the lines in a stars spectrum. By looking at the depth and shape of the different lines in a star's spectrum, the SDSS pipeline can estimate the values of effective temperature and the $\\log g$ for each star. The same is true for parameters relating to the star's chemical composition.
\n", "\n", "**Exercise:** make a scatter plot of the data in the subsample we just selected. The plot should show effective temperature (column teff) on the x-axis and the $\\log g$ value (column logg) on the y-axis. Invert the y-axis so that larger values of log(g) are at the bottom of the plot. Set the symbol sizes to 0.1." ] }, { "cell_type": "code", "execution_count": null, "id": "f29310b3", "metadata": {}, "outputs": [], "source": [ "# code" ] }, { "cell_type": "markdown", "id": "04900b47", "metadata": {}, "source": [ "You should see a diagram with a few long thin blobs of points (a blob of points is technically called a locus of points). These blobs are stars at diffent stages in their evolution.
\n", "\n", "Let's consider giants and dwarfs. These are two different stages in a star's evolution. When a star is burning hydrogen it lies on the main sequence and is known as a dwarf star. Once it runs out of hydrogen fuel in its core it begins to evolve into a giant star. The star stays the same mass when it becomes a giant star, but it gets much bigger and cooler.
\n", "\n", "The acceleration due to gravity in a stars atmosphere depends on the star's mass $M$ and radius $R$:
\n", "$g=\\frac{GM}{R^2}$\n", "\n", "**Exercise:**\n", "Now let's add another parameter to this plot, the metallicity ($[Fe/H]$). This is the log of the ratio of how many iron atoms a star contains compared to many hydrogen atoms. Stars have lots more hydrogen than iron so the scale for $[Fe/H]$ is generally defined so that the Sun would have $[Fe/H]$ of zero. This means that a star with $[Fe/H]=-1$ has one tenth of the metallicity of the Sun. Stars with low $[Fe/H]$ values are said to be metal-poor and those with higher values are said to be metal-rich.
\n", "\n", "Add the following parameters to your plt.scatter command c=data[mask_sp]['fe_h'],vmin=-1.5,vmax=0.5 and add the following line to your plotting code plt.colorbar(label='[Fe/H]')
Have a look at the giant stars, do the metal poor giants have different effective temperatures and $\\log g$ values from metal-rich ones?
" ] }, { "cell_type": "markdown", "id": "abc01752", "metadata": {}, "source": [ "Now let's use the combined SDSS and Gaia data to look at where stars lie in the galaxy. We'll start by creating a SkyCoord object.
" ] }, { "cell_type": "code", "execution_count": null, "id": "b6d8e94b", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# use skycoord to input information\n", "c = SkyCoord(ra=data['ra'][mask_sp]*u.deg, dec=data['dec'][mask_sp]*u.deg, \n", " distance=(1/data['plx'][mask_sp])*u.kpc,\n", " pm_ra_cosdec=data['pmra'][mask_sp]*u.mas/u.yr, \n", " pm_dec=data['pmde'][mask_sp]*u.mas/u.yr,\n", " radial_velocity=data['v_rad'][mask_sp]*u.km/u.s)" ] }, { "cell_type": "markdown", "id": "f0813016", "metadata": {}, "source": [ "This SkyCoord object includes the sky coordinates and distance for each star. These three bits of information allow us to describe the star's position in 3D space. Additionally we have some information about how each star is moving through space. One bit of information is the radial velocity of each star measured by SDSS using the doppler shift in a star's spectral lines. We also have Gaia's proper motions (a star's motion across the sky) for each star which gives us information about the velocities of each star in the plane of the sky. These three bits of information, which allow us to describe the star's motion in 3D space, are also included in the SkyCoord object.
\n", "\n", "So far we've only interacted with each star's position on the sky (its R.A. and Dec.). But since we also have the distance to each star from Gaia's trigonometric parallax we can also start to think in terms of a star's position in space. One coordinate system we can use is called galactocentric coordinates. This is a cartesian coordinate system where $x$, $y$ and $z$ are all zero at the center or the galaxy. The $x$-axis points from the Sun to the Galactic Centre, the $y$-axis points roughly in the direction of galactic rotation at the Sun's position and the $z$-axis is perpendicular to the Milky Way's disk. Both the $x$ and $y$ axes lie in the plane of the Milky Way's disk. In this system the Sun has a negative $x$ position.
" ] }, { "cell_type": "markdown", "id": "561718d4", "metadata": {}, "source": [ "![\"A](\"https://github.com/nialldeacon/SDSS_EPO/blob/main/images/galactocentric_cartesian.png?raw=true\")
\n",
"Image credit: Niall Deacon
The galactocentric coordinate system with the Sun marked with a yellow symbol and the Galactic Center with a red symbol. The centre of the Galaxy's disk is the origin of the coordinate system. The arrows on the edge of the disk mark the direction of galactic rotation.
We can transform our SkyCoord object to the galactocentric coordinate system using the following code.
" ] }, { "cell_type": "code", "execution_count": null, "id": "c4cf7b75", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# transform to galactocentric coordinate\n", "galcen = c.transform_to(coord.Galactocentric(galcen_v_sun=[8, 254, 8] * u.km / u.s, galcen_distance=8.275 * u.kpc))\n" ] }, { "cell_type": "markdown", "id": "f2eb0afd", "metadata": {}, "source": [ "We can then grab the positions and velocities from the SkyCoord object.
" ] }, { "cell_type": "code", "execution_count": null, "id": "07e38db4", "metadata": {}, "outputs": [], "source": [ "x = galcen.x.value \n", "y = galcen.y.value\n", "z = galcen.z.value\n", "\n", "\n", "vx = galcen.v_x.value\n", "vy = galcen.v_y.value\n", "vz = galcen.v_z.value" ] }, { "cell_type": "markdown", "id": "964575cc", "metadata": {}, "source": [ "We can make one slight additional change to our coordinate system, going from a cartesian ($x$,$y$,$z$) system to a cylindrical system ($R$,$\\phi$,$z$) where R is the distance from the Galactic Centre in the plane of the Galaxy and $\\phi$ is the angle around the galactic disk in the direction of Galactic rotation. Here we define$\\phi$ for the Sun to be at $\\phi=180^{\\circ}$ as the Sun has a negative $x$ position in the coordinate system we previosuly defined.
\n", "\n", "![\"A](\"https://github.com/nialldeacon/SDSS_EPO/blob/main/images/galactocentric_cylindrical.png?raw=true\")
\n",
"Image credit: Niall Deacon
The galactocentric cylindrical coordinate system with the Sun marked with a yellow symbol and the Galactic Center with a red symbol. The centre of the Galaxy's disk is the origin of the coordinate system. The arrows on the edge of the disk mark the direction of galactic rotation.
We can transform our positions and velocities to the galactocentric cylindrical coordinate system using the following code (the $z$ coordinate stays the same).
" ] }, { "cell_type": "code", "execution_count": null, "id": "c35629ed", "metadata": {}, "outputs": [], "source": [ "phi = np.arctan2(y, x)\n", "R = np.sqrt(x**2+y**2)\n", "vphi = vx*np.sin(phi)-vy*np.cos(phi)\n", "vR = vx*np.cos(phi)+vy*np.sin(phi)" ] }, { "cell_type": "markdown", "id": "b89ef0f8", "metadata": {}, "source": [ "Let's slice our data up by metallicity. We already made a selection for objects with good quality SDSS and Gaia data mask_sp, now we can use this mask to divide this selection of objects with good quality data a little further.
We can define a slice in our data for the lowest metallicity stars with the following code
" ] }, { "cell_type": "markdown", "id": "618c41b2", "metadata": {}, "source": [ "```python\n", " mask_feh1 = (data['fe_h'][mask_sp]<-1)\n", "```" ] }, { "cell_type": "markdown", "id": "35083513", "metadata": {}, "source": [ "**Exercise:** split your sample up into four slices:\n", "Now let's make some plots of the distribution of the stars in our different metallicity selections.
\n", "\n", "**Exercise:** create histograms of the $z$ position coordinate for your four samples. Ensure that the range of $z$ values plotted is the same for all of the plots. You can add each histogram to the example code below." ] }, { "cell_type": "code", "execution_count": null, "id": "34d3dfca", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(25,40))\n", "plt.subplot(4,1,1)\n", "#plotting code here\n", "plt.subplot(4,1,2)\n", "#plotting code here\n", "plt.subplot(4,1,3)\n", "#plotting code here\n", "plt.subplot(4,1,4)\n", "#plotting code here" ] }, { "cell_type": "markdown", "id": "9131ca79", "metadata": {}, "source": [ "One Possible solution:" ] }, { "cell_type": "markdown", "id": "f794bac6", "metadata": {}, "source": [ "Firstly it's important to note that the underlying dataset we are using is a bit biased. It is not a complete survey of all stars in the Galaxy and will likely have lots more stars closer to the Sun as these are easier to detect.
\n", "\n", "**Exercise:** examine your plots. What happens to the distribution of stars in the $z$ direction? The $z$ coordinate measure the height above and below the plane of the disk of our Galaxy. What can you say about how spread out from the Galactic Disk the different metallicity samples are." ] }, { "cell_type": "markdown", "id": "d14f29d1", "metadata": {}, "source": [ "#answer" ] }, { "cell_type": "markdown", "id": "c2c7636a", "metadata": {}, "source": [ "Now let's look at how the stars in your different samples move in the galaxy.
\n", "\n", "Recall that we defined the positions and motions of each star in terms of three coordinates $R$, $\\phi$ and $z$.
\n", "\n",
"\n",
"**Exercise:** Think about a star orbiting the centre of the Galaxy with an orbital velocity of 200km/s in a circular orbit exactly in the plane of the Galaxy ($z=0$). What, roughly do you think velocities would be in the $R$, $\\phi$ and $z$ directions ($V_R$, $V_{\\phi}$ and $V_z$)?"
]
},
{
"cell_type": "markdown",
"id": "d59cae4f",
"metadata": {},
"source": [
"#answer"
]
},
{
"cell_type": "markdown",
"id": "2ca1a6e2",
"metadata": {},
"source": [
"**Exercise:** Now let's think about a star in an elliptical orbit in the disk. How would the values of $V_R$, $V_{\\phi}$ and $V_z$ look for this?"
]
},
{
"cell_type": "markdown",
"id": "8309cebe",
"metadata": {},
"source": [
"#answer"
]
},
{
"cell_type": "markdown",
"id": "3e334d45",
"metadata": {},
"source": [
"Now let's think about a star in a circular orbit which is tilted with respect to the disk so that it spends half its orbit above the disk and half below. How would the values of $V_R$, $V_{\\phi}$ and $V_z$ look for this?
" ] }, { "cell_type": "markdown", "id": "9f678a56", "metadata": {}, "source": [ "Let's see how the orbits of the stars in your different samples compare to the circular orbit in the disk you just considered.
\n", "\n", "**Exercise:** for each of your four samples, make three histograms, one each for $V_R$, $V_{\\phi}$ and $V_z$.