problem_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework Problem \\#3 - Finding an exoplanet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hand-in format:** IPython Notebook or python program. Submit via email.\n",
    "\n",
    "As a reminder: please make sure your code is clean, documentated, and understandable. Make sure it runs without errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this problem set, we will analyze a dataset that contains the [radial velocity](https://en.wikipedia.org/wiki/Radial_velocity) of a star versus time. The radial velocity is the velocity with which the star is moving towards or away from us (positive velocities indicate it is moving away from us).\n",
    "\n",
    "If a star is not close to any other objects, there is no reason why the radial velocity should change over time. However, in the case where a second object is orbiting a star (such as a planet or another star), the star and the object will both orbit the center of mass of the system. Therefore, the star will show periodic variations in its velocity over time. \n",
    "\n",
    "These changes in velocity then cause a shift in spectral lines via the [Doppler effect](https://en.wikipedia.org/wiki/Doppler_effect), which we can measure with specographs on telescopes.\n",
    "\n",
    "> The smaller the second object, the less the star will be affected. For example, the Earth causes the Sun to change its velocity with an amplitude of 0.1 m/s over 1 year.\n",
    "\n",
    "Watch the following video to see an example of a large planet orbiting a star, and the effect on the observed spectral lines of the sta:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAEBAAMBAQEAAAAAAAAAAAAAAQIDBAYFB//EAD8QAAIBAgIGBgcGBQQDAAAAAAABAgMRITEEEhRBUZEFFVJhcbETIjIzcoGyNEJTc6HBIySCktEGFmJ0NUPh/8QAGQEBAQEBAQEAAAAAAAAAAAAAAAECAwQF/8QALREBAAICAQMCBQIHAQAAAAAAAAERAhIDBAVBFTETITJCUTNxUmGBscHR8BT/2gAMAwEAAhEDEQA/APz8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAHX1fV7UOb/wOr6vahzZdZS4cgOrYKvahzZer6tr60Ob/AMDWS4cgOrYKvahzY2Cr2oc2KkuHKDqegVVm4fqNhq9qHNipLhyg6dhq9qHNjYavahzYqS4cwOnYavahzY2Gr2oc2KkuHMDp2Kp2ocxsVTtQ5ipLhzA6diqdqHMbFU7UOYqS4cwOjY6najzGx1O1HmKkuHODo2Op2o8xsdTjHmKkuHODo2OpxjzGx1OMeYqS4c4OjY6nGPMbHU4x5ipLhzg6NjqcY8xsdTjHmKkuHODo2OpxjzGx1O1HmKkuHODo2OpxjzGx1OMeYqS4c4OjY6nGPMbHU7UeYqS4c4OjY6nGPMbHU4x5ipLhzg6NjqcY8xsdTjHmKkuHODo2OpxjzGx1OMeYqS4c4OjY6nGPMbHU7UeYqS4c4OjY6nGI2Op2o8xUlw5wdL0KqvvQ5jYqnahzFSXDmB07FU7UOY2Kp2ocxUlw5gdOxVO1DmNiqdqHNipLhzA6dhq9qHNjYavahzYqS4cwOnYavahzZdhq9qHNipLhyg6thq9qHNjYKvahzY1kuHKDq2Cr2oc2Xq+r2oc2NZLhyA6tgq9qHNl6vq9qHNjWS4cgOvq+r2oc3/gdX1e1Dm/8DWS4fRADd23ax2cReGQA9XHPIKhCkIBCkAAAACAKEAAAXwsABC3InZp5kAABQWwIAAsAAAGQAAt7q1gIBfAK18XgAe4AAAAAAvhbcAKQGWFmgMblIW/DAoXwAAFBCgAgAighb4WwAoJcAUABF3ZghQKs9z8QQqKKUxKnYCgACANWYAEDIQBmwGsAIV+q8H80QACFIFGQouBBuFg72TfyAEAIorrEd4AAgAAN4gAAQAW5GAAXIqbTGKtIgABgKAAItwQZgV4AgAoWLIUAAwBQQAUAFFBABQAEULEDcBd4DbeYuwBSN3dwEUAryAF7zEpRSk3AAG7u7zDwdsxZu75gTeQpF4kB5KyFsL4AgAAZPAAT5hsNNK4VALgAL4EG7vIoEm3ZZggAreBC3wAgBAAKYgUAAAQoC+AyA3d4C4QAAAAAC3AgAAqBABUwsWN2QAAWAFICq+4ACFKKCFCBSFfcAuUlwgKAAKBezTLhbMIhlfC1l4mJVe19xRQABAwCCAYYgBfCxL2Abu7gL594AvhYKg3b7ggBgZsgUADWWWJADzulZEGbzAC9r94GYAjK72ywIABAAAAFDd8yAKyTtJMm8AAAAALjfENd4EBQBAUgRbZ2BABVgxhdZkK+4Cq18SBd4AFDd228WyAVK7BCgCu2s7ZEBRQAECkKBQEr9yABYMpCoIFvlYgKKAUCAEIG8N3BAKnZprcS7ta+HAAKZkeDs8wAIAwwGCsMiAKEKQgAABa7ssSAPuAEK3fcQAAMgALlYPPAKhUrpvcQoAbu8ItgASBkk7X3ATxLYtu4WAxsLGUrt3buzECAySxxyHFLIDEFaIAYBAiu18MgQqxYEKTeZPDfcBhuIABQAUVAhQBSFasECkAF4YgACgAIouAUCF3EIIAG7gA7XwyBAoAAIAQAW2F8OZBa9yKEAAAB4Y3ANWvinbgQFln39wEIAAAwwt8ytKys7vyCoDfoWiz0zSI0ob82ekf+kXLRnKFRqpbC+8zOUQtPKrEJd5nOlOFaVOUXrxdmlxIqdR2tCWPcVEVioyVKTthnxKqcsPZx/wCSAxsLGxU27Yxx/wCSL6Nu2Ktu9bIDGwsbI07r2o8zNUnbOOfEo0b7mNje6b/4529pBUn/AMc7e0gNTWOVjFo3ulLdHfbAx9FPsSztkBp3Cxk1bB4GLAxZCkAApAigRvfAZMAA8x4ACkLfcUC3IUAFmVq29MgRQEr5DJgW1gBcCoEKnZpoIvgCFKICq292IQQMrthZ+JABCkCq7KKs/kYlGHgBBfCwZAoG8AQgu6xGmnZ5gNYgGN2XzIZRTk1GKbbeSAxG7dY2ejhD3k1dOzjDF88iekUbalOKaecvW88P0AxhCc3aEXJvgrlVPLWnCKe+9/0WInOTSUptpPBbl4GAGaVJW1pSfGyGtTSVoNvfrPD9DAbgr6/+ndIp0tPSmoxusHxPfvTtG0fRHVqTjGKV7n5SnbFGyVarOKU6kpJbmznOETNrEt2n6U9J02tVTajObaV8rnMCo2gjYnd4mKRsUcO8CpKysbIxuYwTRvpxuVGEYNM2KmzohR1lgsTfHR21kWkt8507srp2R9BaNZOTRpqU7XuKLcLiYSVjfU7jRJXIrB1JpYTlnfMxdSW9p43xSZWjFoKrlheUIu7ve1vIl6bxcGnf7rwt8zFkAztTeU2sd63D0Tl7Li8bYSt+jxNYCMpQlB+tFrdijEzjUnCyjJpXva+BkqkX7dOLxxccH/j9ANYRsVOE7ak7Sbfqzwt88vIwlGUX6ytf9QIAAL5AAooACKrAhQL4XDwdh4AAUhWms0ECkBQABBAUgAgZcVg8AqELhbvJuAPwCdmrq/cA1YKgYBBAABnCGtGUpS1YpZve+CE6ralCmtSm3fVTvfx4mVd/w6UYyk4qN8VbF525foaQICkADPABO2QELk7XIUKBBYO4zZBcbdyKiADYk7GyJpRmmB3UKanTlNRlPVzSWS4s3UtV5RSx8Thg2d9LSK8761SUnJWetj5lSX0dFpuTWCfyR9/ReivTUtayVz4mi1pKV5uKurYQS/Y9JoXScI0UpfI386+TnPu+Z0hoMaEWsrZd58DSI4s9B0ppyrSbVrHn9IqRbeAlYcNSxzzZ0VGmc8zDbXjka2zNmDCol3XSzGHzDIAeZCkCBU2iFitZpLMAZwqyiknacFf1ZYr/AOGAA2yppwdSlrOEUta+cW/27zUbdGx0imtWMlJ6rUsnfA121ZNYYMAACi7u8ABAoSwuEBbXeAHAACphJY3dgEEXVbi5JerexC3du4CAEAAuZAIC4WIFCFGF8cUBCFYzyuFQAEECAA3QkqtNUpytJYU22kld4p+ZpqQlTm4TVpJ2aIbqekNRhCrFVaUW2oSdreDA07u8YWyxN3oqVVr0dXVlZuSqeqr8E/8ANjCpRqUlH0kJRUleLazQGsWwuCAV4sCzulvAU3XACzIBSbygZIziYReJnF2A3QOiFXVyONTM4so+lSqu6xOza3GyvkfJpzsZOqW2ad1bSHK+Jx1Kl8zXKqapyFrSzZpkyuRg3cipi8jFmUnb1cPFGBAvhYxKQoXDACABVFzkowi23uWLYAG5aM4pOtOFNOOsru7fyX7h1adO6oQzSWtOzafdw8+8DJR2dOVRL0uKVOUcVdZvngaEJNyk3Jtt5t7yAUAFFA3AIFQCT3AUCzWDzFgLZrNAlypN2wCBU2k1d45kKUTcACCAAAQrzysQKEtcrIACAdgqAAggKQCWG4FsrZ4gYmylWqUdb0c5R1lZ2drowIBuekKXvKFKT1dVNLVt34WVxLZpXcXVp+rgmlK78cLI0h4sK3PR09Z069KajG+er8sbXYlomkRv/ClLVjrNw9ZJcbrwNJN4FlGUHaUXF8GrDxNy0nSW5RVeo9aOq/XeK4eBXpdRuWtClK8dV3pR/wAZ95BoKbpaRGTbejUbuOrhrK3fg8yuro8k/wCW1fVstWo8+ONwNFy3Nutozb/g1l6uH8VZ8fZyMnsjylXj6vZT9bmsArUmZweJVHRvxq3s/hLPh7WXebIR0bD+NU9m7/hrPhn+oQ1sCORlbRsL16i9W7/h5Phn+pLaM1G9aovVvK1PJ8M/1AmtgY3MktG9XWrVVhjamnZ8PaENmvHXqVrW9a0Fg+7HEDW3cxNyloq1LwrSw9a00se7BkjU0ZOOtQnKy9b+Lm+WAGkxublXpRUP5anJxWOtKXrd+DC0pxUNWlRTimrumnfxvcDQbaejV6tvR0aktbK0W7mS0zSIqKhUdPUy1PVf6GqdSc0lOcpJZJu9ijbHRJtx15U6akm7ymsPlmFT0dajqV27+0qcLtc7GgAbo1KMNVxoa7S9b0krpvuStbmybTV1VGM9RKOr6iUbp8bZ5bzVvxyG4IFGBXvSyAhSFACzWYG6xRQG7gIoQbbzCAoAAqBChAApRAAQCbirHBEAEKQKO6eJCh5ICEKAIMsgAqbiFBAtd52IAAaa+ZC2bV9yIADeIAEKsSAKuWRAAK8wQEGV8LbiEKBUZ5GAuBlclyXAGVyEAFIAAABQA8CACgt213BEAAAoCtbMAAUAFd4BOzTLfG5QAARQFjvAFu7W3AMAC7syFCBQrbwUQFatvuGQYgrd8bJdxABCkAADDDAKgLvIBAikCl8LEAvgQGL2yAAm4NjuFm0BAV2wsiAAFdYoAQABQAACptJ9+YRCCluYlAt8CDC3eALjbuBLgAUgXeAK+RCu2BRACp2dwhi8EgL8C52QEAthcAZPg7YcBhYgAqzxyCdn/khQBUQqWBUC4YWIFmBQABQLgCkKAgVK+RClAhQ7EEGF+ADtuVgIxuA34hUBQi0qWGqZopqMbWmvVGqbAdceGJXVq1bMNWV3gZS9uPzMa3uZ+B1jpoqZSYLXwuVU2xD2/wClHTTO+PRYzF255ZU0LR5vJMy2Wbvhln3H0KSRsstTSfh/Yzn0eGLhPPMPnR6PrP7jM+ra6T9R24noqUVZG6tFbO8N68zll0+MPLPW5XVPM9T6S3hSkJ9D6TCLlOEoxW9ntIRXA0dLJdW1f6fNHLLixiHHHuOc5RjTyr6E0pQu6TS4tCHQmkzvqQcrOzsew0j7Ny80aujvZr/mvyRj4cXR6hyaTlTyr6C0q6Xo3d5BdA6U20qbusz2M/tNL4ZfsSl76v8AEvpQ+HDPqPLV1Dx8ugtKSV6bWNh1Dpf4cj2Gk+6Xxw+pGwzPGnqPLV1DxfUOl/hyIug9KbaUHhgz2hppe9r/ABr6UccrhY7jyzfyeS6j0r8NkXQmktXUG0z2Jpoe4gefLlyxWO4clXTya6F0lq6g2mXqTSewz1ND3EPAzPJn1meLX/v5LeT6j0n8NlXQWlfhyPWxNsTzz3HkjxDth1fJk8a+gtKSu6bSMn0BpasnTabwXeewr+4l8vMzq++0f8x/TIkdy5PxD1YcmWTxkugNLjFydOSSV2zL/bum/hS5HtdMX8lX/LfkdKRznuvLEXUO+Ny8F/tzTfwZch/tzTvwZ8j9ASLYz6vy/iHWMH551Bpd7eile9vmV/6d05f+ifI9rBfzMP8AsT+lnfY3n3Xkxr5Q1GD85fQOmLOjPkYPobS1nSlyP0SpFcDlqxXA1h3XOfDUccPAz6Mrwi5SpySWbaMXoFdZwfI9b0jFbJWw+4/I56sVbI9mHW5ZR7NfBh5d6JUWcWYuhNbmfdqpHJUS4Hqx55ySeGHzHTlvuY6p2TNEsztGVuc4RDVZgyZibYmADItsL7isoXMGULa17pWV8UERBAFF3MhSMCMbikIICh4vKwVCohUahYVFCB0xbgAB6cFYS9uPzMK3uZ+BnL24/Mxre5n4Hpj6ckkh7a8EdVM5oe3/AEo6aZ6uP2ebkddI2/c0n4f2NVI2/c0n4f2Mc3s8eXu+pSyRurfZ34rzRppZI3Vvs78V5o8vJ7Pnz9UO6GRz9L/+Nq/0+aOiBz9L/wDjav8AT5o82ftLzcf6mP7t+kfZuXmjT0d7Nf8ANfkjdpH2bl5o09HezX/NfkjH3QsfpS3z+00vhl+xKXvq/wAS+lFn9ppfDL9iUvfV/iX0ovljx/T/ACaT7pfHD6kbTVpPul8cPqRtMZe6fbDFmml72v8AGvpRuZppe9r/ABr6UebNcfaWZpoe4gbjTQ9xA8XI1H0saHuIeBmYUPcQ8DM+bzN+WUTbE0xN0Dw5PXxFf3Evl5myr77R/wAx/TI11/cS+XmbKvvtH/Mf0yMx/v8As+jxs9M+xV/y5eR1I5dM+xV/y5eR1I4ZfTH/AH4ezBkjIxRkcneHDD7TD/sT+lnecEPtMP8AsT+lnedeXx+zUMJnLVOqZy1cjXG3D5fSP2Ot8D8jnq5HT0j9jrfA/I5quR9Pi9nSHFVOOpmdlU46h9HjSXNUNMjdM0yPZi4ZMGYmTMWsDpDjkhSFdr4I0wWd7PMAqKgUXxbeNwA3AACEWfEpErkEKA87gTJ4lJmw8N5bW2WRbmFxc1GVNW2EMLvcS51jmpdll7cfBmFb3U/AtwzrHUxUxSWQfrf0o6abRy4jWlxOuPXREVTnljb6tKSNmstTSfh/Y+P6SfFl9LUV7tkz62MvDhPBfl6ulONlibq01s7x3rzPHrSKy+/LmXaq7w9JLmcsupxnw889FN3b3tOa4mjpZp9G1f6fNHitu0lZVZ8yvT9JkrOrJrg2csuWJcce3ZRlGVvd6R9m5eaNXR3s1/zX5I8V1jpVremnbxC6S0qN7VZK/BmPiRdnp2ek4295P7TS+GX7Epe+r/EvpR4XrPS739NO/iOs9Lv72ePePi/yZ9Nzqre60n3a+OH1I2HgX0npbzqyfzHWmmfjT5mJ5Jnwem51Vves00ve1/jX0o8R1nplr+ml/cOs9LzVWSv3nLK5I7bnF/N7pmmh7iB4rrPS/wAWXMLpHS8lVlzOGXFlksduzqrezoe4h4GZ4nrHSksKsuY6x0v8WXM8ufR55eWvT87u3t4m6J4LrLS/xZcy9Z6X+LPmeee25z5dsOjyx8vd13/Al8vM2VX/ABtH/Mf0yPAdZ6W8HVnzL1ppmfppYd5I7bn+Xpx4ssX6Fpj/AJKv+XLyOlNH5q+ldMas602n3l6101f++f8AcYntXJMVbvFw/TEypn5n1tpufp5/3Mdb6b+PU/uZj0jk/idIze/g/wCZh/2J/Szu1kt5+Y9a6Xg/SzvnmV9K6b+PU/uZvLtOeVfNqM36RUqR4nLVqR4n5/1npbzrT/uI9P0l51Z8zWPaso+5qOWHrukakdkrY/cfkc9WpG2Z5d6XXkrOpJrxJtNV/flzPZh0U4x7tfGh9yrNcTkqSR8301R5yfMOcrJ63yPTjwapPNDqm0aZM1az4kuztGNOc8ls2YkuU3EOcyJq+KugQpWVSTi7yxWS4gtrYEKKNywCABq2Y34ZAAHlaxCsgEBSEC2DxIMwAC8QQKMWwvgABAABAUgVCghAAAAhRgBC2xBAAAAPPBWBclnmQCFA3AAAAAADwBQBCi2AAZZFvfHeQAACgBdtWvgAAKAEAXdYFAJXYAAoAAFd7JW7yBFwK4tZkRSgAUAAVJvICZg+d1hV7MOTJ1hV7MOTM7Q1rL6Qja+OXmfO6wq9mHJk2+r2YcmNoNZfRdtysQ+ft9Xsw5MbdV7MOTJtBrLvBwbdV7MOTC06qvuw5DaDWXeD5+3VezDkxt1Xsw5MbQay772eVxe2RwbbU7MOTG3VOzDkxtC6y7yNHDttTsw5DbanCPIbQay7gcO21OzDkNtqdmHIbQay7rkOLbanCHIm2VOzHkNoKl2g4tsqcI8htlThHkNoKl3Wu8EQ4tsqcI8htlTsx5DaCpdoOJaZUWSiNsqcI8htBUu2wtvOLbKnCPIbXU4R5DaCpdgOPa6nCI2upwjyG0FOwHFtdThHkXbKnCPIbQVLsLvzwOLbKl76sOQ2ypwiNoKl2kONaZUX3YcmNrqcI8htBUu29rWD4nFtdThHkNrqcI8htBTtBxbZU4R5DbKnCPIbQVLtBxbZU4RG2VOEeQ2gqXcDi22p2YcibbU4R5DaCpdwOHbKnCPIbZU7MeQ2gqXek27LFsHDttTsx5DbalrasOTG0Gsu4HDttTsw5DbanZhyG0GsvoYNN7+BDg22p2YcmNuqdmHJjaE1l9DMHBt1Xsw5E26r2YcmNoNZfQKfP26r2YcmNuq9mHJjaDWX0Cnztuq9mHJjb6vZhyY2g1l9Ep83b6vZhyZesKvZhyY2g1l9EHzusKvZhyY6wrJ+zDky7Qay5AAcnQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAf/9k=\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/-BuwWtMygxU\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7f9a28b983a0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo(\"-BuwWtMygxU\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking for small changes for the radial velocities of stars is one method to find planets outside our own Solar System (\"exoplanets\").\n",
    "\n",
    "> The other common technique is the [transit method](https://de.wikipedia.org/wiki/Transitmethode), which we will not consider.\n",
    "\n",
    "## Our goal\n",
    "\n",
    "In this problem set, we want to find out whether a particular star has a companion object, and if so, we want to estimate the probability that the companion object is a planet.\n",
    "\n",
    "The data file required for this exercise is [data/p3_rv.txt](data/p3_rv.txt), a CSV file (whitespace separated) with three columns (and three header lines prefixed by `#`):\n",
    "\n",
    "* column 1: time [days]\n",
    "* column 2: radial velocity [m/s]\n",
    "* column 3: radial velocity uncertainty [m/s]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1 - Visualizing the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start we want to see whether the star does indeed show periodic variations, and if so, we want to measure the period and amplitude.\n",
    "\n",
    "Read in the data and make a plot of the radial velocity in the file. Include axis labels, units, and error bars.\n",
    "\n",
    "Is the time spacing in the data equal? Do you see any clear periodic variations \"by eye\"?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Part 2 - The periodogram\n",
    "\n",
    "To automatically determine the period of a time series, we can use a [periodogram algorithm](https://en.wikipedia.org/wiki/Periodogram). Adopt the following method:\n",
    "\n",
    "1. take an input `period`\n",
    "\n",
    "2. compute the phase of the radial velocity curve: `phase = (time % period) / period`\n",
    "\n",
    "3. sort both the phase and velocity data, by the phase.\n",
    "\n",
    "4. find the \"string length\" $L$ of a line that would join all the sorted points, from left to right. If $\\phi$ is the sorted phase, $v$ is the velocity, and $N$ is the number of points, then\n",
    "\n",
    "  $$L = \\sum_i^{N-1}{\\sqrt{(\\phi_{i+1} - \\phi_{i})^2 + (v_{i+1} - v_{i})^2}}$$\n",
    "\n",
    "---\n",
    "\n",
    "Implement this algorithm, and find the string length for 10,000 periods logarithmically spaced between 1 and 100 days.\n",
    "\n",
    "Make a log-log plot of string length versus period. What do you see? (There should be several clear minima)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Part 3 - Finding the best-fit period\n",
    "\n",
    "Using the above calculations, find the best-fit period of the first clear minimum (optional question: why the first? what do the other minima correspond to?).\n",
    "\n",
    "Then, plot radial velocity (with error bars) versus phase for the best-fit period you derived.\n",
    "\n",
    "> If you have a good best-fit value, you should see something that looks similar to one period of a sinusoidal function.\n",
    ">\n",
    "> This means that there are indeed periodic (in fact. sinusoidal) variations in the radial velocity curve, indicating that a second object is likely orbiting the star!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 - Fitting the radial velocity amplitude\n",
    "\n",
    "Fit the radial velocity versus phase (using your best-fit period), with the functional form: $f(x) = a \\sin{(2\\pi x + b)}$ where $a$ is the amplitude, and $b$ is a phase shift.\n",
    "\n",
    "> Make sure $a$ is positive.\n",
    "\n",
    "Plot the radial velocity data, with error bars, and overplot the best-fit function.\n",
    "\n",
    "Then, derive the uncertainty on $a$, using the covariance matrix returned by the fitting function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now estimate the mass of the planet that is orbiting the star.\n",
    "\n",
    "Given that the radial velocity curve is roughly sinusoidal, this suggests the planet is on a roughly circular orbit.\n",
    "\n",
    "For two bodies (in our case a star and a planet) orbiting a common center of mass in circular orbits, Kepler's 3rd law of motion states that:\n",
    "\n",
    "$$\\frac{4\\pi^2}{T^2} = \\frac{G M_\\star (1 + q)}{R^3}$$\n",
    "\n",
    "where $T$ is the orbital period, $M_\\star$ is the mass of the star, $q$ is the ratio of the object mass to the mass of the star, and $R$ is the distance between the object and the star.\n",
    "\n",
    "In addition, for two objects orbiting a common center of mass, the amplitudes of the velocities $K$ will be related to the mass of the two objects - with the more massive one moving less:\n",
    "\n",
    "$$\\frac{K_\\rm{star}}{K_{\\rm planet}} = q$$\n",
    "\n",
    "Finally, for a circular orbit:\n",
    "\n",
    "$$K_{\\rm planet} T = 2\\pi R$$\n",
    "\n",
    "Substituting this back into Kepler's law and simplifying gives:\n",
    "\n",
    "$$(1 + q)~q^3 = \\frac{T K_{\\rm star}^3}{2\\pi G M_\\star}$$\n",
    "\n",
    "Finally, there is one final important thing to consider, which is that the velocity amplitude we measure for the star is not $K_{\\rm star}$ but\n",
    "\n",
    "$$K_{\\rm obs,star} = K_{\\rm star}~\\sin{i}$$\n",
    "\n",
    "where $i$ is the viewing angle of the system ($i=90^\\circ$ means that we are observing the system *edge-on*, and $i=0^\\circ$ corresponds to a *face-on* orientation where we would not see the star move towards or away from us). $K_{\\rm obs,star}$ is the amplitude of the velocity - the ``a`` parameter you measured above.\n",
    "\n",
    "Therefore, assuming $q << 1$, the final equation is:\n",
    "\n",
    "$$q \\approx \\left(\\frac{T}{2\\pi G M_\\star}\\right)^{1/3} \\frac{K_{\\rm obs,star}}{\\sin{i}}$$\n",
    "\n",
    "---\n",
    "\n",
    "Compute the value of $q$ assuming the values of $T$ and $K_{\\rm obs,star}$ you found above, and adopting $\\sin{i}=1$. Further, assume that $M_\\star$ is the mass of the Sun ($M_\\star=1.989\\times 10^{30} \\,\\rm{kg}$).\n",
    "\n",
    "Using $q$, derive the mass of the object in units of the mass of Jupiter ($M_{\\rm Jupiter}=1.898\\times 10^{27} \\,\\rm{kg}$). Is the object likely to be a planet?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 6 - Quantifying the uncertainty on the planet mass\n",
    "\n",
    "In reality, your derived $K_{\\rm obs,star}$ value has uncertainties, as does the mass of the star, and the viewing angle.\n",
    "\n",
    "Carry out a Monte-Carlo error propagation simulation, to find the likely distribution of masses for the object assuming that:\n",
    "\n",
    "* the error on the radial velocity is the one you derived above, and is normal (i.e. represents a Gaussian)\n",
    "\n",
    "* the mass of the star is sampled from a uniform distribution, between 0.6 and 1.4 times the mass of the Sun.\n",
    "\n",
    "* the viewing angle can be anywhere between 0 and 90 degrees. However, one cannot simply sample $i$ randomly between 0 and 90 because some viewing angles are more likely than others - to do it properly you need to sample $\\cos{i}$ uniformly between 0 and 1, then derive $i$ and so $\\sin{i}$ from this.\n",
    "\n",
    "From the Monte-Carlo simulation, plot a histogram of the probability that the object has a certain mass, and show only the range from 0 to 13 times the mass of Jupiter.\n",
    "\n",
    "What is the (numerical) probability that the object is less massive than 13 times the mass of Jupiter?\n",
    "> This is usually considered to be the upper limit for the mass of a planet.\n",
    "\n",
    "What degree of confidence do we have that the object is a planet, for 1, 3, and 5 sigma [confidence bounds](https://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule)?1/2/3/4/5-sigma confidence terminology?\n",
    "\n",
    "Based on this, what can you conclude about the object?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your solution here\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}