Tutorial_10/Astropy_quantities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Astropy Quantities\n",
    "\n",
    "In astronomy (and most sciences) we typically have to deal with quantities which have some dimension (i.e. unit). The quantity is completely represented by not only a number but also a unit to go along with it. When we consider the same quantity in a different unit system, the quantity itself will not change, although the number used to represent it does change. While solving involved problems with several dimensional quantities, it may become cumbersome to keep track of the quantity if our program stores it as a single number with implied units. Astropy units library is a useful way to keep track of a dimensional quantity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy import units as u \n",
    "#The first step as usual is to import the required packages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can initialise a dimensional quantity as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 42.0 * u.m \n",
    "lengths = [1,2,3]*u.m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the numerical and unit part of q as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(q.value,q.unit) #Ofcourse, we can use q.value with any other command also\n",
    "\n",
    "print(lengths.value,q.unit)\n",
    "print(lengths[0].value , lengths[0].unit)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can multiply or divide quantities having different units but not add/subtract. Ofcourse, quantities with the same unit can be added as usual. Trying to add incompatible quantities will give unit conversion error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance2 =45 * u.m\n",
    "distance1 =5*u.m\n",
    "time = 2 * u.s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speed = (distance2 - distance1)/time\n",
    "print(speed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the units library to interconvert between different unit systems using the to() method. Note that after calling this method, the value stored in quantity.value and quantity.unit don't change. Another method for converting to cgs/SI is also shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(speed.to(u.km/u.hour))\n",
    "print(speed.value)\n",
    "print(speed.unit)\n",
    "print((1*u.Joule).cgs)\n",
    "print((10*u.cm).si)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also define new suitable units as required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mph = u.imperial.mile/u.hour #astropy.units.imperial is the module for fps units\n",
    "print(speed.to(mph))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A special type of unit in the astropy.units package is the dimensionless quantity. This can either be initialised directly or come out from units cancelling out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "refractive_index = 1.5 *u.dimensionless_unscaled\n",
    "print(refractive_index.value,refractive_index.unit) #No unit is printed, as the quantity is dimensionless\n",
    "\n",
    "distance_ratio1 = (1*u.m)/(2*u.m)\n",
    "print(distance_ratio1.value,distance_ratio1.unit)#This is a true unscaled dimensionless quantity\n",
    "\n",
    "distance_ratio2 = (1*u.km)/(2*u.m)\n",
    "print(distance_ratio2.value,distance_ratio2.unit)#This is also dimensionless, but value is scaled\n",
    "\n",
    "print(distance_ratio1.unit == refractive_index.unit)#both are dimensionless unscaled\n",
    "print(distance_ratio1.unit == distance_ratio2.unit)#both are dimensionless, but one is scaled, the other is not"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ofcourse, astropy being *__astro__*py, does support usual astronomical units like parsecs, solar masses, AU etc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 1*u.AU\n",
    "d = 1*u.parsec\n",
    "print(d/r)\n",
    "print(float(d/r)) #Notice how the float method is used to convert dimensionless scaled quatities to numbers\n",
    "                  #The float method will not work with quantities having some dimension\n",
    "\n",
    "AvgDensity_of_sun = 1*u.Msun/(4*np.pi/3*(1*u.Rsun)**3)\n",
    "print(AvgDensity_of_sun.to(u.kg/u.m**3)) #Side note: note how this quite close to density of water"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last remark to make about astropy is that it can also do 'implicit' conversions, like converting wavelength of light to frequency (in vacuum), by specifying the equivalencies parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((1000 * u.nm).to(u.Hz, equivalencies=u.spectral()))  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filters in Astronomy\n",
    "Astronomers often use various filters for taking images. These filters selectively allow light of certain wavelengths to pass through them while rejecting others. Several filter systems are in use, and the filters used for a certain observation are usually specified along with the observation data. Here, we will take up a specific filter system, SDSS (Sloan Digital Sky Survey). In general the analysis of any other filter system will also be similar. In this exercise, we will model a star as a perfect blackbody and try to find the intensity detected in each filter.\n",
    "\n",
    "Recall that Planck's formula for the black body spectrum is given by:\n",
    "$$B(\\lambda,T) = \\frac{2\\pi hc^2}{\\lambda^5}\\left(\\frac{1}{exp\\left(\\frac{hc}{\\lambda k_BT}\\right)-1}\\right)$$\n",
    "$h$ is the planck's constant, $c$ is the speed of light and $k_B$ is Boltzmann constant. $\\lambda$ and $T$ represent the wavelength at which $B$ is measured and temperature of blackbody, respectively. \n",
    "\n",
    "$B(\\lambda,T)d\\lambda$ gives the  energy emitted per unit time per unit area of blackbody in the form of light having wavelengths within a wavelength interval of (small) width $d\\lambda$ centred at said $\\lambda$. We will first write a function that returns B with its proper unit (which is $Wm^{-3}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy import constants as const #The constants package is a prebuilt package from astropy containing\n",
    "                                       #all usual physical constants. They are prebuilt with their units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "h = const.h\n",
    "print(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The complete list of constants available in astropy is available at https://docs.astropy.org/en/stable/constants/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "c = const.c\n",
    "kB = const.k_B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def B(wl,T): #lambda is a reserved keyword in python, so wavelength is assigned variable wl\n",
    "    exponential = 1/(np.exp(h*c/(wl*kB*T))-1)\n",
    "    prefactor = 2*np.pi*h*c*c/wl**5\n",
    "    return prefactor*exponential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.style.use('default')\n",
    "wl0 = np.linspace(100,10000,1000000)*u.nm\n",
    "T0 = 5000 *u.K\n",
    "plt.xscale('log')\n",
    "plt.xlabel('wavelength in nm')\n",
    "plt.ylabel('$B_\\lambda$')\n",
    "plt.plot(wl0.value,B(wl0,T0).value)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to find out the filter response as a function of wavelength. For most standard filters, these are readily available from appropriate sources, as measured experimentally by some astronomers. For our purpose, we will be using the SDSS filters, and in particular ,the data from Jim Gunn,2001 (available from the SDSS website as a fits file, has been already included in the folder with this notebook)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SDSS_filter = fits.open('filter_curves.fits')\n",
    "SDSS_filter.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the primary HDU is empty, since there is no image and only filled HDUs are the binary table HDUs corresponding to the five filters U/G/R/I/Z (ultraviolet/ green/ red/ infrared/ longer than infrared). We can see that each HDU contains 5 columns, and we can find out what they represent by looking at the headers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SDSS_filter[1].header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the first column is the wavelength in Angstroms and the further columns give the transparency of the filter at said wavelength, under various conditions. 'respt' is the quantum efficiency after going through 1.3 airmasses (i.e atmospheric scattering) while 'resbig' is the quantum efficiency including some infrared scattering effects, mainly applicable for extended objects (>80 pixels).'resnoa' is quantum efficiency with no atmosphere and last column represents assumed atmospheric transparency as a function of atmosphere. For our purposes, we will be using mainly the 'respt' column. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = SDSS_filter[1].data\n",
    "G = SDSS_filter[2].data\n",
    "R = SDSS_filter[3].data\n",
    "I = SDSS_filter[4].data\n",
    "Z = SDSS_filter[5].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.plot(U['wavelength'],U['respt'],color='blue') \n",
    "plt.plot(G['wavelength'],G['respt'],color='green')\n",
    "plt.plot(R['wavelength'],R['respt'],color='red')\n",
    "plt.plot(I['wavelength'],I['respt'],color='black')\n",
    "plt.plot(Z['wavelength'],Z['respt'],color='magenta')\n",
    "plt.ylabel('quantum efficiency')\n",
    "plt.xlabel('wavelength in angstrom')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The intensity of light observed in a certain filter is simply the integral of the product of $B(\\lambda,T)$ and quantum efficiency profile of said filter, corrected for distance, integrated over all wavelengths. \n",
    "$$I_{net} = \\frac{r^2}{d^2}\\int_0^\\infty B(\\lambda,T) q(\\lambda) d\\lambda$$\n",
    "where $q(\\lambda)$ is the quantum efficiency, i.e fraction of photons that are let through where r is radius of the star, d is our distance from it and T is its temperature.\n",
    "\n",
    "We will now try to calculate the intensity of the sirius, as calculated in G filter, as an example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelengths = G['wavelength']*u.AA #Angstroms\n",
    "qG = G['respt'] \n",
    "Tsirius = 9940 * u.K\n",
    "r_sirius = 1.711*u.Rsun\n",
    "d_sirius = 2.64 * u.pc "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will estimate the integral by [trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integral = 0*u.W/u.m**2\n",
    "i = 1\n",
    "while i < len(wavelengths):\n",
    "    integral += B((wavelengths[i-1]+wavelengths[i])/2 , Tsirius)*(qG[i-1]+qG[i])/2*(wavelengths[i]-wavelengths[i-1])\n",
    "    i = i + 1\n",
    "intensity = integral*(r_sirius/d_sirius)**2\n",
    "print(intensity.to(u.W/u.m**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your assignment...\n",
    "...should you accept it, is as follows:\n",
    "\n",
    "Arcturus is one of the prominent stars from the Bootes constellation. An astronomer measures the intensity of Arcturus in the G and R filters of the SDSS system and obtains the values of intensity as in the below cell. Estimate the temperature of Arcturus to within $\\pm 50K$. You may assume that Arcturus behaves like a perfect blackbody. The information provided below is sufficient for solving this question and no additional information (such as radius of Arcturus) is required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "intensity_G = 1.3887e-09 * (u.W / u.m**2)\n",
    "intensity_R = 2.5553e-09 * (u.W / u.m**2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}