Tutorial_03/Functions_1-solved.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions in Python\n",
    "\n",
    "Till now, we have seen quite a few functions (especially in [Tutorial 2](../Tutorial_02)). The question arises then on what exactly a fuction is and if we can define functions of our own. This tutorial will use what we have learnt of Python and Numpy and help us take it a step further with user defined functions. \n",
    "\n",
    "A function is simply a piece of code that performs a task, for which it might take some input(s) and might return some output(s) (This will of course depend on the task). The input (and the output) could be a number, an array, a string etc.\n",
    "\n",
    "So why do we need such a thing? If there is a function to evaluate a polynomial at a given point, then we can just code it normally. The answer is robustness. \n",
    "\n",
    "Apart from just being correct, codes should ideally be robust, efficient, easy to debug, and of course, self explanatory for someone else reading it. Functions are very important for acheiving all of the above, with utmost simplicity. Functions are also very useful if you want the same task to be performed repeatedly, with no changes, or maybe with minor modifications. For example, if you want to evaluate that polynomial repeatedly, then it is better to define a function to do so, otherwise simple typing mistakes in one of those lines might break your code. Or if you want to now change your polynomial coefficients, and you have 100 repetitions of the code, then you'll have to go line by line to each of those. \n",
    "\n",
    "**In short, functions eventually make your code easier to read and modify, and much easier to debug.**\n",
    "\n",
    "The following references should be useful:\n",
    "https://docs.python.org/3/tutorial/controlflow.html#defining-functions and https://scipy-lectures.org/intro/language/functions.html for exhaustive description on functions.\n",
    "And of course, Google."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Important Note**: This tutorial is in 2 parts:\n",
    "1. Part 1 (this notebook) deals with basics of function definition in Python. The assignment for the day is also in this notebook, and you should be able to finish the assignment with just this notebook\n",
    "2. Part 2 ([Functions_2](./Functions_2.ipynb)) deals with the nuances of function definition, and can be thought of as a supplement to the main notebook. There may be issues that are discussed here which you might face in the future. So do take out some time to read through this part, just enough to recognise and understand the error messages and how you could proceed in debugging. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall from [Tutorial 2](../Tutorial_02) that `numpy` has several functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sin(π/2) = 1.0\n",
      "exp(2) = 7.38905609893065\n",
      "Mean of [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19] is 9.5\n"
     ]
    }
   ],
   "source": [
    "print(f\"sin(π/2) = {np.sin(np.pi/2)}\") #function that returns the sine of a given number\n",
    "\n",
    "print(f\"exp(2) = {np.exp(2)}\") #function that returns the exponential of a given number\n",
    "\n",
    "arr = np.arange(20)\n",
    "print(f\"Mean of {arr} is {np.mean(arr)}\") #function that returns the mean of a given numpy array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets try to define some functions on our own."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hello World\n"
     ]
    }
   ],
   "source": [
    "def my_print(): #function definition\n",
    "    \"\"\"\n",
    "        Comments made like this are called docstrings\n",
    "        They are useful when you want to document. \n",
    "        Run my_print? in a cell and see how the docstring shows up\n",
    "    \"\"\"\n",
    "    print(\"Hello World\") #function task\n",
    "\n",
    "my_print() #calling the function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above function does not accept any input, and does not give any output. It just performs a simple task, i.e. printing \"Hello World\" whenever its called. (Printing something is not an output, its a task that the function performs. More specifically, if the function gives an output, we call it a 'returned' value)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are some more examples of functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25\n"
     ]
    }
   ],
   "source": [
    "def square(x):\n",
    "    \"\"\"\n",
    "        This function has one argument (this is the proper name for input for a function) in Python)\n",
    "        It returns one value. In this case, the task, and the output are both run in the same line\n",
    "    \"\"\"\n",
    "    return x**2 # Function returns output\n",
    "\n",
    "y = 5\n",
    "print(square(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "210\n"
     ]
    }
   ],
   "source": [
    "def my_product(array):\n",
    "    \"\"\"\n",
    "        This function accepts an array as input, and returns the product of its elements\n",
    "    \"\"\"\n",
    "    product = 1\n",
    "    for i in range(array.size):\n",
    "        product = product*array[i]\n",
    "    return product\n",
    "\n",
    "my_array = np.array([2,5,3,7])\n",
    "print(my_product(my_array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also have functions which take multiple arguments as inputs and also return multiple outputs, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Krittika IITB\n"
     ]
    }
   ],
   "source": [
    "def function(param1,param2):\n",
    "    print(param1,param2)\n",
    "    \n",
    "function(\"Krittika\",\"IITB\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the following code, where we return two values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def powers(x):\n",
    "    return x**2,x**3  # Note that we are returning two values\n",
    "\n",
    "squa, cube = powers(5) # This is called unpacking\n",
    "square_and_cube = powers(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25\n"
     ]
    }
   ],
   "source": [
    "print(squa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "125\n"
     ]
    }
   ],
   "source": [
    "print(cube)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(25, 125)\n"
     ]
    }
   ],
   "source": [
    "print(square_and_cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unpacking is not new in particular, and allows you to assign each element of a list, tuple or array to different variables, which can be independently manipulated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above is a very basic introduction to functions. The tutorial is by no means complete, and is intended to give you a start, and teach you the stuff that is most commonly used. You are encouraged to Google things, and explore yourself, to know more about functions - one of the most powerful tools developed in programming."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your assignment ...\n",
    "...should you choose to accept it is the following:\n",
    "\n",
    "You are given a 'galaxies.csv' file which contains data for galaxies observed by the Sloan Digital Sky Survey(SDSS) - Mapping Nearby Galaxies At Apache Point Observatory(MaNGA). The file contains basic properties of each galaxy observed, namely  the following:\n",
    "\n",
    "1. serialno - A Serial Number we have assigned to each galaxy\n",
    "2. objra (in degrees) - The Right Ascension of the galaxy\n",
    "3. objdec (in degrees) - The Declination of the galaxy\n",
    "4. redshift - The observed redshift in the spectra of the galaxy\n",
    "\n",
    "Right Ascension and Declination are a set of [coordinates]( https://skyandtelescope.org/astronomy-resources/right-ascension-declination-celestial-coordinates/) used to describe positions of objects in the sky in astrophysics. It is a spherical coordinate system similar to the Longitude and Latitude used to describe positions on the Earth's Surface. (in fact RA and Dec can be visualized as projections of the Longitude and Latitude on the sky)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Redshift](https://www.esa.int/Science_Exploration/Space_Science/What_is_red_shift) describes how much the spectrum of a galaxy is shifted with respect to when the galaxy is stationary with respect to us. \n",
    "\n",
    "This is due to the Doppler Effect; the apparent wavelength of an object moving away from us is larger than its actual wavelength( wavelength as seen from the object's rest frame). And vice versa for an object moving towards us, and this phenomena is called blueshift."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, for galaxies that are located reasonably far away (like the galaxies observed by SDSS-MaNGA), we can assume that the redshift is primarily caused due to the expansion of space itself (Yes!!! we have entered the realm of cosmology), so we can use cosmological relations to calculate the distances to these objects, given their redshifts.\n",
    "\n",
    "Your task is the following:\n",
    "\n",
    "1. Compute the distance to a galaxy, in suitable units, given the galaxy's serial number.\n",
    "\n",
    "2. Compute the number of galaxies observed by SDSS-MaNGA, in some given distance interval.\n",
    "\n",
    "3. Compute the physical separation between two galaxies, given their serial numbers. Note: The two galaxies might be at different radial distances.\n",
    "\n",
    "The above tasks must be done using functions (which take the input of one (or two, as required) arguments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You can use the following information:\n",
    "1. [Hubble–Lemaître law](https://www.pnas.org/content/112/11/3173): Observations by Edwin Hubble in the early 1900s suggested that far away galaxies all move away from us with velocities (called recession velocity) approximately proportional to their distances. $$v = H_0 \\times d$$\n",
    "Where $v$ is the recession velocity, and $H_0$ is the Hubble constant (for this assignment, you can take it to be 70 (Km/s)/Mpc) and $d$ is the instantaneous distance between the source and the observer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "2. Using the Doppler shift, we can measure the velocities of these galaxies using the formula \n",
    "$$v = z \\times c$$\n",
    "Where $v$ is defined earlier, $c$ is the speed of light in vacuum and $z$ is the redshift. This is an approximate relation, valid for z<<1. Almost all redshifts observed by SDSS-MaNGA are less than 0.15, hence you can use this relation as an approximation. It is encouraged to derive this equation from the standard formula for Doppler effect of light.\n",
    "\n",
    "<We are assuming that the redshifts are small enough that the non linear cosmological effects are negligible.>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. The angular separation between two objects, given their RA and Dec: \n",
    "$$ \\cos(\\theta) = \\sin(\\delta_1)\\sin(\\delta_2) + \\cos(\\delta_1)\\cos(\\delta_2)\\cos(\\alpha_1 - \\alpha_2) $$\n",
    "Where $\\delta_1$ and $\\delta_2$ denotes the declination of the two objects respectively, $\\alpha_1$ and $\\alpha_2$ denotes the right ascension of the two objects respectively, and $\\theta$ denotes the angular separation between them. The above relation can be derived using laws of spherical trigonometry. If you are interested in a derivation, you can check out: http://www.castor2.ca/07_News/headline_062515.html<br><br>\n",
    "You could alternatively try to convert the coordinates of the galaxies to a Cartesian System. In this case, the RA would be analogous to the azimuthal angle, and the delination would be the co-polar (or elevation) angle. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We have seen unpacking before; in this context, it is equivalent to taking a transpose.\n",
    "serialno, objra, objdec, redshift = np.loadtxt('galaxies.csv', delimiter = ',', unpack = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "H0 = 70 #km/s/Mpc\n",
    "c = 3E5 #km/s\n",
    "\n",
    "def distance_from_z(z):\n",
    "    return z*c/H0 #distance in Mpc\n",
    "\n",
    "def get_distance(num): #num denotes the serialno of the galaxy\n",
    "    index = np.where(serialno == num)[0][0] #np.where returns a tuple, whose first element is the array of required indices\n",
    "    return distance_from_z(redshift[index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance of galaxy 1000 is 160.18 Mpc\n"
     ]
    }
   ],
   "source": [
    "serial = 1000\n",
    "# The :2.2f allows you to format the string. \n",
    "print(f\"Distance of galaxy {serial} is {get_distance(serial):5.2f} Mpc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_to_z(distance): #distance in Mpc\n",
    "    return H0*distance/c\n",
    "\n",
    "def count_galaxies(distance1, distance2): #distances in Mpc, distance1 should be less than or equal to distance2\n",
    "    z1 = distance_to_z(distance1)\n",
    "    z2 = distance_to_z(distance2)\n",
    "    indices = np.where(np.logical_and(redshift>=z1, redshift<=z2))\n",
    "    return len(indices[0]) #indices is a tuple, we need the size of its first element (the array of required indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Galaxies observed by MaNGA between 200 Mpc and 210 Mpc is 118\n"
     ]
    }
   ],
   "source": [
    "dist1 = 200\n",
    "dist2 = 210\n",
    "print(f\"Number of Galaxies observed by MaNGA between {dist1} Mpc and {dist2} Mpc is {count_galaxies(200,210)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cosine_law(side1, side2, cosine_value):\n",
    "    return (side1**2 + side2**2 - 2*side1*side2*cosine_value)**0.5\n",
    "\n",
    "def cosine_angular_separation(coordinate1, coordinate2): #angular separation between 2 galaxies, given their ra, dec values\n",
    "    ra1 = coordinate1[0]*np.pi/180\n",
    "    dec1 = coordinate1[1]*np.pi/180\n",
    "    ra2 = coordinate2[0]*np.pi/180\n",
    "    dec2 = coordinate2[1]*np.pi/180\n",
    "    cosine_value = np.sin(dec1)*np.sin(dec2) + np.cos(dec1)*np.cos(dec2)*np.cos(ra1-ra2)\n",
    "    return cosine_value\n",
    "\n",
    "def get_coordinates(num): #get ra and dec of a galaxy, given its serialnum\n",
    "    return objra[np.where(serialno == num)], objdec[np.where(serialno == num)]\n",
    "\n",
    "def get_separation(num1, num2): #get distance between two galaxies, given their serial numbers\n",
    "    coordinate1 = (get_coordinates(num1)[0][0], get_coordinates(num1)[1][0])\n",
    "    coordinate2 = (get_coordinates(num2)[0][0], get_coordinates(num2)[1][0])\n",
    "    cosine_value = cosine_angular_separation(coordinate1, coordinate2)\n",
    "    side1 = get_distance(num1)\n",
    "    side2 = get_distance(num2)\n",
    "    separation = cosine_law(side1, side2, cosine_value)\n",
    "    return separation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The distance between Galaxies 1000 and 2000 is 169.9829610909196 Mpc\n"
     ]
    }
   ],
   "source": [
    "serial1, serial2 = 1000, 2000\n",
    "print(f\"The distance between Galaxies {serial1} and {serial2} is {get_separation(serial1, serial2)} Mpc\") "
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}