example_SNR_findchirp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial will give an example of how to use Riroriro to calculate the signal-to-noise ratio (SNR) of the gravitational waveform from a given binary merger.\n",
    "\n",
    "We import the data of a synthetic gravitational waveform generated using the gravitational wave functions of Riroriro (demonstrated in the tutorial example_GW), and the FINDCHIRP trace data for the corresponding waveform. The example waveform is included in this repository, but the FINDCHIRP data files are not included in Github for file size reasons and can instead be found at this link: https://www.dropbox.com/sh/lt8109c0c2eb10c/AADBnYXMGmWOEqCK-XWFZuJLa\n",
    "\n",
    "You can use any Riroriro gravitational waveform in the format output by waveform_exporter for inputarray (if you choose a system different than the default, remember to also select the corresponding FINDCHIRP trace). We use the data for a binary black hole merger with log(Mc/Msun) = 1.4 and mass ratio q = 0.8. We put this example system at a distance of 1000 Mpc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "inputarray = np.genfromtxt('example_1.4_0.8.dat')\n",
    "findchirp_array = np.genfromtxt('findchirp_traces/mergerwave_1.40000_0.800000.dat')\n",
    "\n",
    "d = 1000.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import riroriro.snrcalculatorfuns as snr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the redshift from the input distance and then adjust the strain amplitude to account for the distance and the frequency to account for the effect of redshift."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = snr.polynomial_redshift(d)\n",
    "inputarray = snr.redshift_distance_adjustment(inputarray,d,z) #note self-reference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We determine the minimum and maximum of the frequency range of the waveform, which we will need at a later step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "freqmax, freqmin = snr.frequency_limits(inputarray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We perform an approximation of a Fourier transform on the strain amplitude, using the FINDCHIRP data for the system with the same parameters as calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "fourieramp = snr.findchirp_fourier(inputarray,findchirp_array,d,z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below plot illustrates the relationship between the transformed amplitude and frequency, similar to figure 18 of the thesis \"Gravitational Waveform Synthesis\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAELCAYAAAAhuwopAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoYklEQVR4nO3de3xdVZn/8c+TW9Nrrr2kbdL0ChQphaa0RVFwKiBQYLxwUQcRBmZ0ZETnpzIjv5eoMy90RH8zBUasgIgoCApSGBAERVBaaNOW0gKVNjRtWuglvTdt0yTP74+9kx5CLuckZ+ec5Hzfr9d5JWeds/d6EnrysNbaez3m7oiIiCQiK9UBiIhI/6PkISIiCVPyEBGRhCl5iIhIwpQ8REQkYUoeIiKSMCUPERFJWE6qA+gLpaWlXllZmeowRET6lerq6p3uPrKj1/pF8jCzi4HzgVHA7e7+dEdtnR1fWVnJ8uXL+yJUEZEBw8xqO3st8mkrM7vbzLab2Zp27eea2TozW29mN3R1Dnf/rbtfA1wJXNpZm4iI9I2+GHncA9wG3NvaYGbZwO3AR4A6YJmZLQaygZvbHX+Vu28Pv78xPC5WR20iIhKhyJOHuz9vZpXtmk8D1rt7DYCZPQBc5O43Axe0P4eZGfBd4El3X9FZm4iI9I1UrXmMAzbHPK8D5nTx/uuA+UCBmU1x9zs6aWtjZtcC1wJUVFQkM3YRkYyXquRhHbR1ur2vuy8EFnbX1u71RcAigKqqKm0dLCKSRKm6z6MOKI95Ph7YmuxOzGyBmS3au3dvsk8tIpLRUpU8lgFTzWyimeUBlwGLk92Juz/m7tcWFBT06Pim5hYeXbWFbfsOJzkyEZH+LfJpKzO7HzgTKDWzOuCb7n6XmX0ReIrgCqu73X1tBH0vABZMmTKlR8ev3bqPLz2wCoBJpUOZO7mEuZNKmDephJHDByUvUBGRfsYyoZJgVVWV9+QmweYW57Wt+1hSs5OlNbt4+a1dHDjSBMCUUcOYN6mEeZNLmDOxmJJhSiYiMrCYWbW7V3X42kBOHjEjj2vefPPNXp+vqbmFNVv3sWRDPUtq6lm+cRcNjc0AHDd6OPPCkcncScUUDsnrdX8iIqmUscmjVU9HHt052tzC6rq9LK2pZ8mGepbX7uLw0RbM4IQxI5g3OZjimj2xmILBuUnvX0QkSkoeESWP9o40NfPK5r3hyGQnKzbtobGphSyDE8cWvCuZDBvUL7YVE5EMlrHJI9nTVok6fLSZlZv2sKSmnqUb6lm5eTdHm53sLOOkcQVt01yzK4sYkqdkIiLpJWOTR6u+Gnl051BjM9W1u4Nprpp6Xtm8h6YWJyfLOLm8sG0BftaEIvJzs1MdrohkOCWPNEke7R080sTy2t1tC/BrtuylucXJy85iZnkhc8NprlMqCpVMRKTPKXmkafJob//hoyzfuJsl4QL82q17aXEYlJPFqRVFwZrJ5BJOHl9IXo6KQIpItDI2eaR6zaO39h46ystv7Wq7muv1d/bhDvm5WVRNKG5bM5kxvoDcbCUTEUmujE0erfrLyKM7exoaWVpzLJms27YfgCF52cyuLA7ufp9cwvvGjiBHyUREeknJY4Akj/bqDxzhpbd2ta2ZrN9+AIDhg3KYPbG4bQH+hLIRZGd1tJGxiEjnukoeuj60HysZNojzTirjvJPKANi+/zBLa4Jk8lJNPX94IyjAOCI/h9MmlrTdZ3L8mOFkKZmISC8M6JFHf1/z6K139h5um+JaUlPPpl0NABQNyWVOazKZXMLUUcMICjOKiByjaasBOm2VqC17DrFkQ31bQtmy5xAAJUPzgj25wpHJ5JFDlUxERMlDyaNjm3c1tI1Klmyo552wbsnI4YOYN6mkbQG+smSIkolIBtKah3SovHgI5cVDuGR2Oe5ObX1DWyJZUlPP4leC4o5jRuS3rZfMm1xCefGQFEcuIqmmkYd0yN3ZsONg275cS2vqqT/YCMC4wsFto5J5k0sYVzg4xdGKSBR6PG1lZvOAzwBnAGXAIWAN8L/Afe7eL4qDK3n0nrvz5vYDwahkQz1L36pnT8NRACqKh7SNSuZOKmFMQX6KoxWRZOhR8jCzJ4GtwKPAcmA7kA9MA84CFgA/dPek1x5Plky/2ipKLS3OG+/sb9vk8aWaevYdDqosTiwd2jYymTupmFHDlUxE+qOeJo9Sd9/ZzYm7fU860Mgjes0tzutvH6uy2L5k79xJxcybVMrcSSrZK9JfJOVqKzMbChx29+ZkBtcXlDz6XlNzC2u37mtbgF+mkr0i/U5PRx5ZwGXAp4HZwBFgELADeAJY5O79Yi5IySP1Ykv2Lq0Jkklryd7jx4xoWzM5TSV7RdJGT5PHn4BnCNY81rh7S9heTLDm8SngEXe/L5Kok0jJI/00NrXwSt2etgX46k27OyzZW1VZxPB8JRORVOhp8sh196PdnLjb96QDJY/0d/hoM6s272lbM1m1aQ+NzS1kZxnvG1fQNjKpmlDEUNV/F+kTvVrzMLO7gFvdfVVM203uflMyg4ySkkf/c6ixmRWbdrdtp7KqXcne1gX4WROKGJynKosiUeht8qgDdhJclntv2LbC3U9NeqQRUfLo/xoam95VZfFVlewViVxvk8cK4EzgF8Am4EvAMnc/JclxJp3u8xi4Wkv2tt5nsmZLULI3LyeLUysKmTeplHmTS5hZrpK9Ij3V2+SxsjVRmNlNwEeAMneflOxAo6KRx8C399BRlr21q21k0nHJ3mJmjC9UyV6ROPV2Y8S2O8jd/SYzWw58JVnBiSRDweBc5k8fzfzpo4GgZG9rlcWlNfV8/6l1QFCyt6ryWJVFlewV6RltjCgZobVkb2stkzfDkr3DBuVw2sTitgX46WNVslekVY9GHma2H+gosxjg7j4iSfGJRK59yd4d+4+0rZcs3dBxyd65k4o5YcwIlewV6UBcI4/YdY/+SCMP6c62fe8u2VtbH5TsLRySy5yJrdNcpUwbrZK9kjmSUQxq4M9tSUYbPSKfi2aO46KZ4wDYGpbsbV2Af2rtNqB9yd5iJo9UMpHMpFt1RTowtnAwH581no/PGg+EJXvDKa4lNfX876tvA0HJ3rmTjlVZVMleyRRdrXl8LOZpYbvnuPvDkUUlkmbaSvZWHSvZ27pmsmRDPY+1K9nbugBfXjxYyUQGpK72tvppF8e5u18VTUjJpzUPiZK7U7PzYNs010s19ew88N6SvXMnFTO+SPXfpf9ISj2PVDGzi4HzgVHA7e7+tJmdQHCneynwrLv/qKtzKHlIX3J31m8/0DYqWVpTz+6wZG958eC2Ka55k0pVslfSWk931f0M8MvWrdg7eH0ywZ3mf+6i47uBC4Dt7v6+mPZzgf8GsoE73f27cfwQRcAt7n51TFsW8JPYto4oeUgqtbQ467btb0skL721i72HgmTSWrJ37qTgLniV7JV00tOrrUqAlWZWDVQTFIHKB6YAHyLYLPGGbvq+B7gNuDcmmGzgdoJtTuqAZWa2mCCR3Nzu+KvcfXv4/Y3hca3nuTDs/7ZuYhBJqaws44SyEZxQNoKrPjCxrWRv66XBj7+ylftf3gTA5JFD20YlKtkr6azLaavwD/2HgfcDZcAh4HXgSXffFFcHZpXA460jDzObB9zk7ueEz/8VwN3bJ47W4w34LvB7d3+mg9f/193P76D9WuBagIqKilm1tbXxhCvS55qaW3gtpv77srd2cTAs2Ttt9LC2aa45E0soGqqSvdJ3enqH+ffc/etmVpjk2h3jgM0xz+uAOV28/zpgPlBgZlPc/Q4zOxP4GEFZ3Cc6OsjdFwGLIJi26n3YItHIyc5ixvhCZowv5B8+NJmjzS28umVv2zTXg8vr+NmSWpXslbTS1ZrHq8CpwEu9qd3Rwcjjk8A57v734fO/A05z9+t62kcXfWtLdun3GptaWF13rMpide1ujjQF9d9PHHssmcyuLFbJXkmqni6Yf59g2mco0BD7EgnsbdXbaatk0IK5DCRHmppZuWlP25rJynYle4N7TIJkopK90hu9refxqLtf1IvOK3l38sgB/gr8DbAFWAZ8yt3X9rSP7ih5yEB2+GgzK2qPVVl8pW4PR5uDkr0zxhe0LcCrZK8kqqcjj18SXO10uBcd309QhbAU2AZ8093vMrPzgP8iuMLqbnf/j5720U3/mraSjNPQ2ER17e62aa7VdUHJ3txsY2Z5IfPCvblOrShSyV7pUk+Tx+nAxQTTVI+7+58iizBiGnlIJjtwpIllG3exNFyAf7WTkr0nlxcwKEfJRI7p7bRVDnAecAawC3jI3dcnPcoIaOQh8l77Docle8ORyWtvHyvZO2tCUdsCvEr2StK2JzGzAuBWYBLwNXd/MTkhRksjD5HO7Wlo5OWY+u9vvLMfCEr2nlpRRFVlEbMri5lZXqgF+AzT6+RhZuXAF4AFBHeb/5xgT6l+cf+EkodI/HYdbOSlcMfgZRt388Y7wcgkO8uYXjaiLZlUTShi1AhtpzKQJSN5zAPOBv7T3Q8lOb7IaNpKpPf2HT7Kyk17WL5xF8s27mLV5j0cPhpseVdRPKQtmcyuLGJS6TCV7R1AerpgfjnwtLvXh8+Hu/v+6MKMjkYeIslztLmFtVv3tSWT6trdbVvQFw7JpWpCEVXhyOR94wp0RVc/1tONEScAD5lZLvAs8KSZvdxfpqpEJBq52VnMLC9kZnkhf3/GJNydjfUNQSLZuJtltbt45vVgP9OcLOP4suHMLC/klPIiZlYUMrFkqEYnA0A8V1sNJ9hb6lzgNIKNEX8HPOXu2yKPMAk08hDpW/UHjrBi0x5Wbd7Nyk17WF23lwNHmgAYkZ/DyeWFnFJeyPvGFdDU4uw8cISGcDNIgGGDcigakkfRkFzGFQ1mbOFgXfmVAkktBmVm04GPAme3bjGSrrTmIZIemlucDTsOsGrTHlZu3sPKTbv567b9tMT55yc7yxhXOJhJI4cyvWwE08eOYHrZCCo1iolU0isJmtnx7v5GryPrIxp5iKSfg0eaWLdtP4NzsykdNoihg7IxjBZ39h9uYndDI7sPNlK35xCb6huo3dXA+u0HeHPbfprCrDM8Pye4nHhCEbMqizilXFuwJFNP1zy68jRQ0fOQRCTTDR0U/OHv7LXOSvQeaWpm/fYDrN26j5Wb9lBdu4sf/H4HAHnZWcyeWMQZU0dyxtRSThgzQiOTiHR1tdXCzo4BPhvvrrrpQCMPkYFtb8NRVmzazYsbdvLCmzvbbnQsHTaIM6aWcsbUUj44bSSlqsyYkJ5eqrsf+BfgSAcv/8DdS5MXYjS05iGSmbbtO8wLb+7k+b/u4M/rd7LrYCNmUDWhiHNOHMM5J46hvHhIqsNMez1NHn8AbuxoCxIze8vdJyY3zOho5CGSuVpanLVb9/HM69t4au07baOSE8pGcM6JoznnxDEcP2Y4QcVridXT5FEMHHb3hg7f0I8oeYhIq9r6gzy9dhu/W/sOKzbtxh0mlAzh/JPKuGjmOI4bMzzVIaaNpF9t1d8oeYhIR7bvP8wzr23nyTVv8+KGeppbnOPHDOfCmWNZMGNsxk9t9XTk8RiwCPidux9t99ok4Epgo7vfndxwk0/JQ0S6s2P/EZ549W0Wv7KV6trdAMyaUMRFM8dy3kllGbnY3tPkMQb4CvBxgjoeO4B8oBLYANzm7o9GEXCyKXmISCI272pg8StbWbxqK+u27Sc7yzjruJFcUlXOWcePypi73ZOxq24lUAYcAv7aX9ZBdLWViPTWG+/s45GVW3h4xRZ27D9C6bBBfHzWOC6pKmfyyGGpDi9SWvPQyENEeqmpuYXn1u3gV8s384c3ttPc4lRNKOKS2eWcf1LZgCyUpeSh5CEiSbR9/2EeWbGFXy3fTM2OgwzNy+ZvTx3HFfMqmTZ64FytpeSh5CEiEXB3qmt3c//Lm3ls9VYam1o4fXIJV8yrZP4Jo8jp52sjyVjzGAxUuPu6ZAfXF5Q8RCRq9QeO8Kvlm7lvSS1b9x5mbEE+nz29kk/NqWB4fm6qw+uRXiWPcNH5FiDP3Sea2Uzg2+5+YdIjjYiSh4j0labmFp59Yzv3/GUjS2rqGT4oh0/PncBV76/sdzXfe5s8qoEPA8+5+ylh22p3n5H0SCOi5CEiqbC6bg8//lMNT655m5ysLD526ji+cOYUKkr6x82Hvd2Svcnd92rfFxGRxMwYX8jtnz6VjTsPcuefa3hoeR2/rq7jktnlXPfhKZQVDE51iD0Wz2rOGjP7FJBtZlPN7FbgPZslpiMzW2Bmi/bu3ZvqUEQkg1WWDuXfLz6JF752Fp+aU8FDyzfzoe8/x3cef42dBzrauDz9xTNtNQT4BnA2QS2Pp4DvuPvh6MNLDk1biUg62byrgYXPvslvVtQxJC+H6+dP5bOnV6bdneu6VFfJQ0TS0IYdB/j3x1/jj+t2MG30MG668EROn5w+pZJ6szFip5lFV1uJiPSeu/Ps69v51uNr2bzrEBeePJZvXXgiRUPzUh1ajxfMbwm/fgwYA9wXPr8c2Ji06EREMpiZMX/6aD4wtZQ7/rSB2/+4nqU19dzyyZP54LSRqQ6vU/GseTzv7h/sri2daeQhIv3F2q17uf6BVby5/QBXnl7JDR89nvzc7JTE0tXII57VmZFh/Y7Wk00E0jcdioj0YyeOLeCx6z7A595fyT0vbuSTdyxh2770uz4pnuTxZeA5M3vOzJ4D/ghcH2VQIiKZLD83m28uOJE7r6iiZscBLrztz6zZkl63HHSbPNz9d8BU4Evh4zh3fyrqwEREMt386aP59edPJycri0/c8SK/f21bqkNq023yMLMrgEuBk8PHpWGbiIhE7ISyETzyT6dz3OjhfOEX1TyTJgkknmmr2TGPM4CbgD67TNfMLjazn5jZo2Z2dkz7UDOrNrML+ioWEZFUGDU8n5///Ryml43gC79YwXPrtqc6pLimra6LeVwDnALEdQGymd1tZtvNbE279nPNbJ2ZrTezG7rp/7dhv1cSjIBafR14MJ44RET6uxH5udx71Rymjh7GtT+vprp2d0rj6cm98A0EayDxuAc4N7bBzLKB24GPAtOBy81supmdZGaPt3uMijn0xvA4zGw+8BqQHuM3EZE+UDAkl/uunkNZQT6fv6+a7Sm8CqvbXXXb3WmeRfAH/6F4Tu7uz5tZZbvm04D17l4Tnv8B4CJ3vxl4zxSUBdv5fhd40t1XhM1nAUPDWA6Z2RPu3tLuuGuBawEqKiriCVdEJO0VDc3jx383i4/9z4t8/hcruP+aueTl9P2eWPFsyX5LzPdNQK271/Wiz3HA5pjndcCcLt5/HTAfKDCzKe5+h7t/A8DMrgR2tk8cAO6+CFgEwU2CvYhXRCStHD9mBN//xMn80y9X8F/P/JWvnXt8n8cQT/I4z92/HttgZt9r35aAjgqDdLWH1kJgYSev3dNlR0EVxAVTpkxJJD4RkbR3/owynls3nh8/X8P5M8o4cWxBn/Yfz1jnIx20fbQXfdYB5THPxwNbe3G+Trn7Y+5+bUFB3/5SRUT6wo3nT6doSB5f/81qmprfMwETqU6Th5l93sxeBY4zs9Uxj7eA1b3ocxkw1cwmmlkecBmwuBfn65SKQYnIQFYwJJebLpzOmi37eHjllj7tu6uRxy+BBQR/2BfEPGa5+2fiObmZ3Q8sIUhAdWZ2tbs3AV8kKCr1OvCgu6/txc/QKY08RGSgO/+kMmaML+C/n3mTI03NfdZvV8nD3X0j8E/A/pgHZlYcz8nd/XJ3L3P3XHcf7+53he1PuPs0d5/s7v/Rux9BRCRzmRn/cvZxbNlziAeXbe7+gCTpbuQBUA0sD79WxzxPe5q2EpFM8MGppZxcXshPX9xIX1WH7TR5uPsF4deJ7j4p/Nr6mNTZcelE01YikgnMjCvmTqBmx0Fe3FDfJ312tWB+alePPolORETicv6MMoqG5PLLlzb1SX9d3efxgy5ec+DDSY4l6XSfh4hkivzcbM6fUcZvqrfQ0NjEkLx4buPrua6mrc7q4pH2iQM0bSUimeW8k8o4dLSZ59btiLyveOp55JvZV8zsYTP7jZldb2b5kUcmIiIJmTOxhOKheX1S8yOecc29BJfo3ho+vxz4OfDJqIISEZHEZWcZp08u4S8bduLuBPvKRiOe5HGcu58c8/yPZvZKVAElk9Y8RCTTnD65lMdXv82GHQeZMmpYZP3Es7fVSjOb2/rEzOYAf4ksoiTSmoeIZJrTJgb3cK/cFG2xqHiSxxzgRTPbaGYbCbYb+ZCZvWpmvdnjSkREkmxi6VCG5GWzduu+SPuJZ9rq3O7fIiIi6SA7y5heNoK1W6PdWSOeGua1wD6gAChpfbh7bfha2tL2JCKSid43roC1W/fR0hLdViXxlKH9DnAlsIFjRZv6xU2C7v4Y8FhVVdU1qY5FRKSvTBs9nIbGZt7ed5hxhYMj6SOeaatLgMnu3hhJBCIiklQTSoYAsKm+IbLkEc+C+RqgMJLeRUQk6SqKg+SxeVdDZH3EM/K4meBy3TXAkdZGd78wsqhERKTHygryyckyancdjKyPeJLHz4DvAa8CfVskt5d0k6CIZKKc7CxGj8jnnb1Hun9zT/uI4z073X1hZBFESAvmIpKpSoblUX8wtcmj2sxuJqhlHjtttSKyqEREpFdKhuax80B01znFkzxOCb/OjWnrF5fqiohkqpJhg1j3zv7Izt9t8nD3syLrXUREIlEyLI+dBxsj2103rlJTZnY+cCLQVsfD3b+d9GhERCQpiofk0djUQkNjM0MHJb+qYDzFoO4ALgWuA4ygjseEpEciIiJJ05owDh5piuT88dwkeLq7XwHsdvdvAfOA8kiiERGRpBg6KBuAg43NkZw/nuRxKPzaYGZjgaPAxEiiSTJtjCgimWpoXupHHo+bWSHwfWAFsBG4P5JokkzFoEQkU0U9bRXP1VbfCb/9jZk9DuS7u/5XXkQkjbUlj8YUJY9Y7n6EmBsFRUQkPQ3NC9Y8DhxJ3ZqHiIj0M4NyguTR2BTNloRKHiIiA1BuTnBj4NHmaJJHp9NWZnZqVwdqbysRkfSVmx2MDZr6OnkAPwi/5gNVwCsENwnOAF4CPhBJRCIi0mutyaOxOZo65p1OW7n7WeG+VrXAqe5e5e6zCDZKXB9JNCIikhR5YfKIatoqnjWP49391dYn7r4GmBlJNCIikhQ52eGaRwoXzF83szvN7Ewz+5CZ/QR4PZJoOmBmF5vZT8zsUTM7O2w708xeMLM7zOzMvopFRKS/yMmKdsE8nuTxOWAt8CXgeuC1sK1bZna3mW0P65/Htp9rZuvMbL2Z3dDVOdz9t+5+DXAlwQaNENQTOUCwHlMXTywiIpnEzMjLzuJoSzRrHvHcYX443Fn3CXdfl+D57wFuA+5tbTCzbOB24CMEf/iXmdliIBu4ud3xV7n79vD7G8PjAF5w9z+Z2Wjgh8CnE4xLRGTAy822yO7z6DZ5mNmFBPta5QETzWwm8G13v7C7Y939eTOrbNd8GrDe3WvC8z8AXOTuNwMXdNC/Ad8Fnmy9PNjdW38bu4FB3cUhIpKJsrKM5lSNPIBvEvzBfw7A3Vd1kBASMQ7YHPO8DpjTxfuvA+YDBWY2xd3vMLOPAecAhQQjm/cws2uBawEqKip6Ea6ISP+UFUEFwVbxJI8md9+bxDKGHZ2o09To7guBhe3aHgYe7qoTd18ELAKoqqqKJvWKiKSxLIMW7+P7PGKsMbNPAdlmNtXMbgVe7EWfdby7mNR4YGsvztcp1fMQkUxmZilNHtcR1C8/QlDHYx/BVVc9tQyYamYTzSwPuAxY3IvzdUr1PEQkk2UZRJQ74rraqgH4RvhIiJndD5wJlJpZHfBNd7/LzL4IPEVwhdXd7r420XPH2f8CYMGUKVOiOL2ISFoLRh7RnDueq62qgH8DKmPf7+4zujvW3S/vpP0J4Im4o+whd38MeKyqquqaqPsSEUk3BnhEQ494Fsx/AXwVeBWI5oJhERFJuiyz1E1bATvcPZI1iahp2kpEMlmUV1vFdZ+Hmd0JPEtMCdrwctm0pmkrEclkKV3zINjH6nggl2PTVk4391mIiEhqmYF3fhtdr8STPE5295Mi6T1imrYSkUwW5ZpHPPd5LDWz6dF0Hy3d5yEimcxSvObxAeCzZvYWwZpHePVX95fqiohI6qTsaqtwR9t/IChF2+9o2kpEMlmUI48up608uLvk/7l7bftHJNEkmaatRCSTBdNE0Zw73jWP2dF0LyIiUcmKcGPEeNY8zgL+0cw2AgfRmoeISL8QYTmPuJLHR6PrXkRE+qNup63C9Y1CYEH4KOwvax6q5yEiEo1uk4eZfYlgc8RR4eM+M7su6sCSQQvmIiLRiGfa6mpgjrsfBDCz7wFLgFujDExERNJXPFdbGdAc87yZjuuQi4hIhohn5PFT4CUzeyR8fjFwV2QRiYhI2ut05GFmEwHc/YcEO+vuAnYDn3P3/+qT6HpJC+YiItHoatrq1wBm9qy7r3D3he7+3+6+so9i6zUtmItIpkvF3lZZZvZNYJqZfeW9AfkPowlJRESSwSJcnu5q5HEZcJggwQzv4CEiIhmq05GHu68Dvmdmq939yT6MSURE0lw8l+pOM7MRFrjLzFaY2dmRRyYiImkrnuRxlbvvA84GRhJcefXdSKMSEZG0Fu9NggDnAT9191fQTYIiIhktnuRRbWZPEySPp8xsONASbVjJofs8RESiEU/yuBq4AZjt7g1AHsHUVdrTfR4iItHodnsSd28xs23AdDOLZzsTERFJE06KKgmGu+heCrzGsQ0SHXg+kohERCQpUl1J8GLgOHc/El0YIiLSn8Sz5lED5EYdiIiI9B/xjDwagFVm9izQNvpw93+OLCoREUlr8SSPxeFDREQEiO9qq5/1RSAiItJ/xHO11VTgZmA6kN/a7u6TIoxLRETSWDwL5j8FfgQ0AWcB9wI/jzKoWGZ2sZn9xMwebd2Q0cyyzOw/zOxWM/tsX8UiIiKBeJLHYHd/FjB3r3X3m4APx3NyM7vbzLab2Zp27eea2TozW29mN3R1Dnf/rbtfA1xJcL8JwEXAOOAoUBdPLCIimSgVlQRbHTazLOBNM/sisAUYFef57wFuIxitAGBm2cDtwEcI/vAvM7PFQDbB9Fisq9x9e/j9jeFxAMcBS9z9x2b2a+DZOOMREZEkiCd5XA8MAf4Z+A7B1FVcU0Xu/ryZVbZrPg1Y7+41AGb2AHCRu98MXND+HGZmBFvAP+nuK8LmOqAx/L65/TEiIhKtLpNHOEq4xN2/ChwgORsijgM2xzyvA+Z08f7rgPlAgZlNcfc7gIeBW83sDDrZJsXMrgWuBaioqEhC2CIi0qrT5GFmOe7eZGazzMzckzZz1tFuK52e290XAgvbtTUQ7PbbKXdfBCwCqKqqimjWT0QkM3U18ngZOBVYCTxqZg8BB1tfdPeHe9hnHVAe83w8sLWH5+qSmS0AFkyZMiWK04uIZKx4rrYqBuoJrrC6AFhAB2sTCVgGTDWziWaWB1xGRHewq56HiEg0uhp5jDKzrwBrCKaVYqeb4poGMrP7gTOBUjOrA77p7neFV209RXCF1d3uvrYnwcfRv0YeIiIR6Cp5ZAPDSHCN4l1vcr+8k/YngCfiOUdvuPtjwGNVVVXXRN2XiEgm6Sp5vO3u3+6zSCKgkYeIZLqorhbqas0jwhpUfUNrHiKSySzCUoJdJY+/iaxXERHp1zpNHu6+qy8DiYKZLTCzRXv37k11KCIiA0o8l+r2W5q2EhGJxoBOHiIiEg0lDxERSdiATh5a8xARicaATh5a8xARicaATh4iIpkuqkqCSh4iIgNUlHd6D+jkoTUPEZFoDOjkoTUPEZFoDOjkISIi0VDyEBGRhCl5iIhIwpQ8REQkYQM6eehqKxHJZGMK8hk5PC+Sc5tHdQdJGqmqqvLly5enOgwRkX7FzKrdvaqj1wb0yENERKKh5CEiIglT8hARkYQpeYiISMKUPEREJGFKHiIikrABnTx0n4eISDQGdPLQrroiItHIiJsEzWwHUJvAIQVAPMOVeN9XCuxMoP+BIt7fT1/qq5iS3U9vz9fT4xM9LpH3x/PeTP3sQHp8fia4+8gOX3F3Pdo9gEVJft/yVP9M6fx7HIgxJbuf3p6vp8cnelwi74/nvZn62UnGf/OoHwN62qoXHkvy+zJVOv5++iqmZPfT2/P19PhEj0vk/en47yOdpPXvJyOmrVLNzJZ7J/vDiEjn9NlJXxp59I1FqQ5ApJ/SZydNaeQhIiIJ08hDREQSpuQhIiIJU/IQEZGEKXmkgJlNMrO7zOzXqY5FpD8xs4vN7Cdm9qiZnZ3qeDKZkkeSmNndZrbdzNa0az/XzNaZ2XozuwHA3Wvc/erURCqSXhL87PzW3a8BrgQuTUG4ElLySJ57gHNjG8wsG7gd+CgwHbjczKb3fWgiae0eEv/s3Bi+Limi5JEk7v48sKtd82nA+nCk0Qg8AFzU58GJpLFEPjsW+B7wpLuv6OtY5Rglj2iNAzbHPK8DxplZiZndAZxiZv+amtBE0lqHnx3gOmA+8Akz+8dUBCaBnFQHMMBZB23u7vWA/uGLdK6zz85CYGFfByPvpZFHtOqA8pjn44GtKYpFpD/RZyfNKXlEaxkw1cwmmlkecBmwOMUxifQH+uykOSWPJDGz+4ElwHFmVmdmV7t7E/BF4CngdeBBd1+byjhF0o0+O/2TNkYUEZGEaeQhIiIJU/IQEZGEKXmIiEjClDxERCRhSh4iIpIwJQ8REUmYkodkFDNrNrNVMY/KVMeULGZ2ipndGX5/pZnd1u7158ysqovjHzCzqVHHKQOD9raSTHPI3Wd29IKZGcG9Ty19G1LS/Bvw7704/kfA14BrkhOODGQaeUhGM7NKM3vdzP4HWAGUm9lXzWyZma02s2/FvPcbYXGiZ8zsfjP7P2F72//Rm1mpmW0Mv882s+/HnOsfwvYzw2N+bWZvmNkvwsSFmc02sxfN7BUze9nMhpvZC2Y2MyaOv5jZjHY/x3Bghru/EsfPfGHMyGudmb0VvvQCMN/M9D+V0i39I5FMM9jMVoXfvwV8GTgO+Jy7fyEsbTqVoJ6EAYvN7IPAQYL9lU4h+NysAKq76etqYK+7zzazQcBfzOzp8LVTgBMJNvv7C/B+M3sZ+BVwqbsvM7MRwCHgToLKedeb2TRgkLuvbtdXFbCmXdulZvaBmOdTANx9MeE+UWb2IPCnsL3FzNYDJ8fxs0mGU/KQTPOuaatwzaPW3ZeGTWeHj5Xh82EEyWQ48Ii7N4THxbNJ39nADDP7RPi8IDxXI/Cyu9eF51oFVAJ7gbfdfRmAu+8LX38I+L9m9lXgKoLKe+2VATvatf3K3b8Y87M+F/uimX2N4PcRW5FvOzAWJQ/phpKHSDCqaGXAze7+49g3mNn1QGcbwTVxbAo4v925rnP3p9qd60zgSExTM8Fn0Trqw90bzOz3BFUoLyEYZbR3qF3fXTKzvwE+CXyw3Uv54blEuqQ1D5F3ewq4ysyGAZjZODMbBTwP/K2ZDQ7XFxbEHLMRmBV+/4l25/q8meWG55pmZkO76PsNYKyZzQ7fPzxm/eFOgiJIy9y9fclWCHaenRLPD2hmE4D/AS5x9/aJYhqg3WulWxp5iMRw96fN7ARgSbiGfQD4jLuvMLNfAauAWoLF5Va3AA+a2d8Bf4hpv5NgOmpFuCC+A7i4i74bzexS4FYzG0wwApgPHHD3ajPbB/y0k2PfMLMCMxvu7vu7+TGvBEqAR8Kfcau7n2dmowmmsd7u5ngRbcku0hNmdhPBH/Vb+qi/scBzwPGdXUpsZl8G9rv7nT3s48vAPne/q8eBSsbQtJVImjOzK4CXgG90cw/Kj3j3Wkqi9gA/68XxkkE08hARkYRp5CEiIglT8hARkYQpeYiISMKUPEREJGFKHiIikjAlDxERSdj/ByPwo1nbJlGsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(inputarray[:,1],fourieramp)\n",
    "plt.xlabel('Frequency (Hz)')\n",
    "plt.ylabel('Transformed amplitude (1/√f(Hz))')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compare the transformed amplitude to the noise spectra of three detectors in order to calculate the SNR. These spectra were sourced from https://dcc.ligo.org/LIGO-T1500293/public and are also included in this repository for convenience.\n",
    "\n",
    "You can also use different noise spectra files, as long as they have the frequency values in the first column and ASD noise levels in the second.\n",
    "\n",
    "We combine the SNR values from each individual detectors to obtain both the quadrature SNR (if all three detectors were operating simultaneously) and the mean-of-individual SNR (if a non-specific single detector was operating)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Livingston SNR is 23.3, Hanford SNR is 18.1, Virgo SNR is 8.9.\n",
      "The quadrature SNR is 30.8.\n",
      "The mean-of-individual SNR is 16.8.\n"
     ]
    }
   ],
   "source": [
    "ind_SNR = np.empty((3))\n",
    "for i in range(3): #three detectors\n",
    "    if i == 0:\n",
    "        noisearray = np.genfromtxt('noise_spectra/o3_l1.txt') #LIGO Livingston\n",
    "    elif i == 1:\n",
    "        noisearray = np.genfromtxt('noise_spectra/o3_h1.txt') #LIGO Hanford\n",
    "    elif i == 2:\n",
    "        noisearray = np.genfromtxt('noise_spectra/o3_v1.txt') #Virgo\n",
    "    #the noise spectra are not included in this repository for copyright reasons,\n",
    "    #but the ones used for this example can be found at the following link:\n",
    "    #https://dcc.ligo.org/LIGO-T1500293/public\n",
    "    \n",
    "    noise_freq_amp = snr.amplitude_interpolation(inputarray,fourieramp,noisearray,freqmax,freqmin)\n",
    "    ind_SNR[i] = snr.individual_detector_SNR(noisearray,noise_freq_amp)\n",
    "\n",
    "quad_SNR = np.sqrt(ind_SNR[0]**2 + ind_SNR[1]**2 + ind_SNR[2]**2)\n",
    "mean_SNR = (ind_SNR[0] + ind_SNR[1] + ind_SNR[2])/3\n",
    "\n",
    "print(\"Livingston SNR is %.1f, Hanford SNR is %.1f, Virgo SNR is %.1f.\\n\"\n",
    "     \"The quadrature SNR is %.1f.\\nThe mean-of-individual SNR is %.1f.\"\n",
    "      % (ind_SNR[0],ind_SNR[1],ind_SNR[2],quad_SNR,mean_SNR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 4
}