astrophysics_activity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Useful imports\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import libstempo.toasim as LT\n",
    "import libstempo.plot as LP\n",
    "from matplotlib.pyplot import *\n",
    "\n",
    "from astropy.coordinates import SkyCoord"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the stochastic background, we expect a number of closer, louder sources to stand out to us. These are high-mass, inspiraling black hole binaries with orbital frequencies around $10^{-8}$ Hz. Their orbital frequency evolves too fast for the binary to contribute significantly to the stochastic background, but, relative to a typical human lifespan, it evolves very little. As such, we term these signals \"continuous gravitational waves\" since PTA experiments don't witness these binaries changing in orbital frequency. \n",
    "\n",
    "Some of these potential CW sources are close and bright enough to be observed electromagnetically. Optical and X-ray observations of active galactic nuclei and quasars can indicate if a galaxy has a black hole at the center large enough for us to detect. We can use this supplemental astronomical data to inform our searches, giving the search pipeline a heads up about where a signal might be lurking. The information goes both ways, too. By achieveing upper limits from our latest dataset we've been able to make statements about the kinds of binaries which are definitely not present in our data, which in turn constrains the mass, frequency, etc. of the candidate binaries astronomers are also interested in.\n",
    "\n",
    "In this exercise, we begin by looking at a continuous wave signal injected into some simulated pulsar timing data. First, we inspect the residuals by eye and do some quick estimates to give us an idea of how the injected signal gets jumbled. And then we run a full MCMC to see how much more information we can recover. Lastly, our results will be used to estimate some properties of the binary and it's host galaxy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.0 Make personalized directory "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-2-3e9052183b4f>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-2-3e9052183b4f>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    INIT = #initials here\u001b[0m\n\u001b[0m                         ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "INIT = #initials here\n",
    "DIR = './cw_test'.format(INIT)\n",
    "!mkdir './cw_test'.format(INIT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.0 Simulate Pulsar Data\n",
    "### 1.1 Generate fake `.par` files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this function creates a parameter file for a simulated pulsar\n",
    "# it randomizes sky location, name, pulse frequency and parallax\n",
    "\n",
    "def make_fake_pulsar(DIR):\n",
    "    '''\n",
    "    Makes a fake pulsar par file\n",
    "    '''\n",
    "    output = \"MODE 1\\n\"\n",
    "    \n",
    "    # Sphere Point Picking\n",
    "    u = np.random.uniform()\n",
    "    v = np.random.uniform()\n",
    "    phi = 2*np.pi*u #using standard physics notation\n",
    "    theta = np.arccos(2*v-1) - np.pi/2\n",
    "\n",
    "    c = SkyCoord(phi,theta,frame='icrs',unit='rad')\n",
    "    cstr = c.to_string('hmsdms')\n",
    "    #print cstr\n",
    "    RAJ = cstr.split(\" \")[0].replace(\"h\",\":\").replace(\"m\",\":\")[:-1]\n",
    "    DECJ = cstr.split(\" \")[1].replace(\"d\",\":\").replace(\"m\",\":\")[:-1]\n",
    "    cstr = cstr.replace(\" \",\"\")\n",
    "    name = \"J\"+RAJ[0:2]+RAJ[3:5]+DECJ[0]+DECJ[1:3]+DECJ[4:6]\n",
    "\n",
    "    output += \"PSR      %s\\n\"%name\n",
    "\n",
    "    \n",
    "    output += \"PEPOCH   50000.0\\n\"    \n",
    "    output += \"POSEPOCH   50000.0\\n\"\n",
    "\n",
    "    period = 0.001*np.random.uniform(1,10) #seconds\n",
    "    output += \"F0       %0.10f 1\\n\"%(1.0/period)\n",
    "\n",
    "    output += \"RAJ      %s 1\\n\"%RAJ\n",
    "    output += \"DECJ     %s 1\\n\"%DECJ\n",
    "\n",
    "    dist = np.random.uniform(0.1,5) #kpc\n",
    "    output += \"PX       %0.5f 1\\n\"%(1.0/dist) #mas\n",
    "\n",
    "    filename = \"%s/%s.par\"%(DIR,name)\n",
    "    print(filename)\n",
    "    with open(filename,'w') as FILE:\n",
    "        FILE.write(output)\n",
    "\n",
    "    return filename.encode('ascii','ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "par1 = make_fake_pulsar(DIR)\n",
    "par2 = make_fake_pulsar(DIR)\n",
    "par3 = make_fake_pulsar(DIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2  Generate fake `.tim` files with an injected single GW source\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this function simulates a pulsar's timing file\n",
    "# and injects the TOAs with a CW source of given chirp and frequency [Hz]\n",
    "\n",
    "def observe(par,noise=0.5,mass=5e9,fgw=1e-8):\n",
    "    ''' Noise in microseconds, mass in solar masses'''\n",
    "    # let's set up some fake TOAs\n",
    "    t = np.arange(53000,56650,30.0) #observing dates for 10 years\n",
    "    t += np.random.randn(len(t)) #observe every 30+/-1 days\n",
    "    psr = LT.fakepulsar(parfile=par,\n",
    "                        obstimes=t,\n",
    "                        toaerr=noise)\n",
    "    LT.add_equad(psr,equad=noise*1e-6,seed=42) #this is mildly correlated white noise \n",
    "    \n",
    "    # Let's add a source at the Virgo cluster: 12h27m +12d43'\n",
    "    LT.add_cgw(psr, gwtheta=(12+43.0/60)*np.pi/180, gwphi=(12*15+27.0/60)*np.pi/180, mc=mass, dist=15, fgw=fgw, phase0=0, psi=0, inc=0, pdist=1.0,\n",
    "               pphase=None, psrTerm=True, evolve=True,\n",
    "               phase_approx=False, tref=0)\n",
    "    psr.savetim('%s.tim'%par.split('.')[0])\n",
    "    print('%s.tim'%par.split('.')[0])\n",
    "    return psr\n",
    "\n",
    "#Virgo cluster location in terms of celestial coordinates; useful for reading enterprise output\n",
    "print((12+43.0/60)*np.pi/180,(12*15+27.0/60)*np.pi/180)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psr1 = observe(par1.decode('utf-8'))\n",
    "psr2 = observe(par2.decode('utf-8'))\n",
    "psr3 = observe(par3.decode('utf-8'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the timing residuals with this injected signal with the `LP.plotres` function. Tip: if you are unsure how to use a particular imported function, you can write and execute `??<name of function>` to print a description of the function's input and output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.0 Quick Estimates \n",
    "\n",
    "### 2.1 Approximate Gravitational Wave Strain\n",
    "\n",
    "We can estimate the GW strain rather trivially from the timing perturbations. Conveniently, the root-mean-square (RMS) residuals are provided at the top of the libstempo plots. The strain $h \\sim \\Delta L/L$ is the fractional change in distance towards our pulsars. Therefore, we have\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\Delta L & \\sim & hL \\sim hcT \\\\\n",
    "& \\approx & 9.5~\\mathrm{m}~\\left(\\frac{h}{10^{-16}}\\right)\\left(\\frac{T}{10~\\mathrm{yr}}\\right),\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "where $L\\sim cT$ is the length scale form the light/GW travel time. The timing precision is then simply\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\Delta t & \\sim & \\frac{\\Delta L}{c} \\sim hT \\\\\n",
    "& \\approx & 32~\\mathrm{ns}~\\left(\\frac{h}{10^{-16}}\\right)\\left(\\frac{T}{10~\\mathrm{yr}}\\right).\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Using the RMS as the timing perturbation scale, what is your estimate of the GW strain?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.0 Construct CW Signal Search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# More imports for GW detection codes\n",
    "import sys\n",
    "# Enterprise\n",
    "sys.path.append(\"/home/jovyan/work/shared/enterprise/\")\n",
    "from enterprise.pulsar import Pulsar\n",
    "import enterprise.signals.parameter as parameter\n",
    "from enterprise.signals import utils\n",
    "from enterprise.signals import signal_base\n",
    "from enterprise.signals import selections\n",
    "from enterprise.signals.selections import Selection\n",
    "from enterprise.signals import white_signals\n",
    "from enterprise.signals import gp_signals\n",
    "import corner\n",
    "from PTMCMCSampler.PTMCMCSampler import PTSampler as ptmcmc\n",
    "import enterprise_extensions.model_utils as utils\n",
    "import enterprise_extensions.models as models\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Add pulsars to model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create enterprise Pulsar objects for each of your pulsars. You can make Pulsar objects with\n",
    "#psr = Pulsar(<PARFILE>, <TIMFILE>, ephem=\"DE436\", timing_code=\"tempo2\")\n",
    "\n",
    "\n",
    "# Here you should append your pulsar objects to this list for use by the model later\n",
    "psrs = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Form timing model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create white noise parameter priors\n",
    "efac = parameter.Constant(1)\n",
    "equad = parameter.Uniform(-8.5,-5.0)\n",
    "\n",
    "##### Signals below #####\n",
    "\n",
    "# white noise parameters\n",
    "ef = white_signals.MeasurementNoise(efac=efac)\n",
    "eq = white_signals.EquadNoise(log10_equad=equad)\n",
    "\n",
    "# timing model\n",
    "tm = gp_signals.TimingModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Initialize CW source parameter priors & add to signal model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The strain amplitude is a function of a few of the source parameters,  \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "h_{c} = \\dfrac{2 \\mathcal{M}_c^{5/3}(\\pi f_{\\rm GW})^{2/3}}{d_{L}}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "(Note: This formula is an averaged over sky position.)\n",
    "\n",
    "You'd anticipate that more massive systems are \"louder,\" as are closer systems. This is reflected in the strain's dependence on the binary's effective mass, called the chrip mass ($\\mathcal{M}_{c}$) and inverse dependence on a particular distance to the binary, called the luminosity distance ($d_L$).\n",
    "\n",
    "If we have no prior knowledge about our continuous wave source, the only information we can recover from our search is about the strain (i.e., $h_{c}$). And since that is a mix of chirp mass and luminosity distance, we need some other information in order to get estimates of these two values independently. Fortunately, there is a way to untangle these values: through the inclusion of the so-called pulsar term.\n",
    "\n",
    "The full influence of a gravitational wave passing through our pulsar array is a combination of the measurements of the signal at the Earth and at each pulsar individually. If the wave jiggles Earth, then the times of arrival for all the pulsars in the array will change. We encapsulate this particular portion of the signal effect into the \"Earth-term\". But, every pulsar experiences a jiggle too, and since the pulsar's pulses take time to reach us, we essentially see the signal in the past, before it reached Earth. This is the \"pulsar-term\" in our residuals. In combination, we uncover the time evolution of the signal. \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\t\\Phi(t) = \\Phi_0 + \\Phi_p + \\frac{2\\pi}{32} \\mathcal{M}^{-5/3} \\left[ f(t_{p,0})^{-5/3} - f(t_p)^{-5/3} \\right]\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\t\\Phi_p = \\frac{2 \\pi}{32} \\mathcal{M}^{-5/3} \\left[ f_0^{-5/3} - f(t_{p,0})^{-5/3} \\right]\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Theoretically, we understand this to provide information on the chirp mass of the system. By including both Earth- and pulsar-terms into our search, we're better able to recover the source properties. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# continuous GW parameters\n",
    "# note that we are pre-initializing them with names here so that they will be shared\n",
    "# across all pulsars in the PTA\n",
    "\n",
    "# Our standard CW search looks for a GW with a specific frequency and we hold log10_fgw constant,\n",
    "# as commented below. However, we will do a search in frequency\n",
    "#freq = 8e-09\n",
    "#log10_fgw = parameter.Constant(np.log10(freq))('log10_fgw')\n",
    "\n",
    "\n",
    "cos_gwtheta = parameter.Uniform(-1, 1)('cos_gwtheta') #position of source\n",
    "gwphi = parameter.Uniform(0, 2*np.pi)('gwphi') #position of source\n",
    "log10_mc = parameter.Uniform(7, 10)('log10_mc') #chirp mass of binary\n",
    "log10_fgw = parameter.Uniform(-10,-7)('log10_fgw') #gw frequency \n",
    "phase0 = parameter.Uniform(0, 2*np.pi)('phase0') #gw phase\n",
    "psi = parameter.Uniform(0, np.pi)('psi') #gw polarization \n",
    "cos_inc = parameter.Uniform(-1, 1)('cos_inc') #inclination of binary with respect to Earth \n",
    "\n",
    "log10_h = parameter.LinearExp(-18, -11)('log10_h') #gw strain (linear exponential for an upper limit calculation)\n",
    "\n",
    "\n",
    "# define CGW waveform and signal\n",
    "cw_wf = models.cw_delay(cos_gwtheta=cos_gwtheta, gwphi=gwphi, log10_mc=log10_mc, \n",
    "                     log10_h=log10_h, log10_fgw=log10_fgw, phase0=phase0, \n",
    "                     psi=psi, cos_inc=cos_inc)\n",
    "cw = models.CWSignal(cw_wf, ecc=False, psrTerm=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Define the full model (includes noise + GW signals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# full signal\n",
    "s = ef + eq + tm + cw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 With model and pulsars, create `enterprise` PTA object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize PTA\n",
    "model = [s(psr) for psr in psrs]\n",
    "pta = signal_base.PTA(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.6 Prepare the MCMC sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare sampler initial condition\n",
    "x0 = np.hstack(p.sample() for p in pta.params)\n",
    "ndim = len(x0)\n",
    "\n",
    "# initial jump covariance matrix\n",
    "cov = np.diag(np.ones(ndim) * 0.1**2)\n",
    "\n",
    "# parameter groupings\n",
    "groups = utils.get_parameter_groups(pta)\n",
    "\n",
    "# define where you want to put the chains from the MCMC\n",
    "chaindir = 'chains/'\n",
    "\n",
    "# set up jump groups by red noise groups (need better way of doing this)\n",
    "sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, groups=groups, \n",
    "                 outDir=chaindir, resume=False)\n",
    "\n",
    "# write parameter file for convenience\n",
    "filename = chaindir + '/params.txt'\n",
    "np.savetxt(filename,list(map(str, pta.param_names)), fmt='%s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add prior draws to proposal cycle -- this helps prevent the sampler's walkers \n",
    "# from getting trapped in local minimum in parameter space\n",
    "jp = utils.JumpProposal(pta)\n",
    "sampler.addProposalToCycle(jp.draw_from_prior, 5)\n",
    "sampler.addProposalToCycle(jp.draw_from_cw_log_uniform_distribution, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.7 Sample!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100000\n",
    "sampler.sample(x0, N, SCAMweight=30, AMweight=15, DEweight=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.0 Visualize Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the MCMC chains and parameter names (if you need them)\n",
    "\n",
    "chain = np.loadtxt(chaindir + '/chain_1.txt')\n",
    "params = list(np.loadtxt(chaindir + '/params.txt', dtype='str'))\n",
    "\n",
    "# Value of burn-in to apply to the chains.\n",
    "burn = int(0.25*chain.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convenience function to help plot the marginalized parameter distributions\n",
    "def plot_param(name):\n",
    "    '''Given one of the CGW function names above, plot the distribution'''\n",
    "    hist(chain[burn:,params.index(name)], 50, normed=True, lw=2, color='C0', histtype='step')\n",
    "    if \"log10\" in name:\n",
    "        xlabel(r\"$\\log_{10} \\mathrm{%s}$\"%(name.split('_')[-1]))\n",
    "    else:\n",
    "        xlabel(name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "corner.corner(chain[burn:,:-5], labels=pars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.0 Post-Processing Analysis \n",
    "\n",
    "### 5.1 Upper and Lower Limits\n",
    "\n",
    "The 11-year dataset was evaluated to find the 95% upper strain limit. This means with 95% confidence 95% of our strain samples lie below a value $h_{c,95}$. This can be transalted to a luminosity distance, and we can limit the mass of the SMBHB at the center of the Virgo cluster. \n",
    "\n",
    "\n",
    "From your analysis find the 95% upper limit of the strain amplitude and from that determine the 95% lower limit on the luminosity distance to the Virgo cluster.\n",
    "\n",
    "For reference, the chirp mass for two equal mass black holes $m_1 = m_2 = M/2$  is\n",
    "\n",
    "$$\\mathcal{M}_c = \\frac{(m_1 m_2)^{3/5}}{(m_1+m_2)^{1/5}} = \\frac{(M/2)^{6/5}}{M^{1/5}} = 2^{-6/5} M \\approx 0.435 M.$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.constants import M_sun, G, c\n",
    "from astropy import units as u\n",
    "\n",
    "def d_luminosity(h, M, f):\n",
    "    \"\"\"\n",
    "    Eq. 17 from 5-yr paper, luminosity distance given \n",
    "    gravitational wave strain, f in Hz, M in solar masses, \n",
    "    and h (not log10(h)).\n",
    "    \n",
    "    :param h: characteristic strain\n",
    "    :param M: binary chirp mass, units of solar masses\n",
    "    :param f: search frequency, units of Hz\n",
    "    \n",
    "    :returns: 95% lower limit on luminosity distance in Mpc\n",
    "    \"\"\"\n",
    "    d = (2 * (M * M_sun * G)**(5/3) * (np.pi * f * u.s**(-1))**(2/3) / (h * c**4)).to(u.Mpc)\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is your 95% strain upper limit and 95% sky-averaged lower luminosty distance limit?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Orbital Separation\n",
    "\n",
    "As the two black holes have a large separation, we can assume the system is in the Newtonian limit and use the Keplerian orbital frequency to relate $\\Omega^2 = GM/a^3$, where $M$ is the total mass of the system and $a$ is the semi-major axis. The gravitational wave frequency is twice the orbital frequency, and so\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "f & = & 2\\left(\\frac{\\Omega}{2\\pi}\\right) \\\\\n",
    "& = &\\frac{1}{\\pi} \\left(\\frac{GM}{a^3}\\right)^{1/2} \\\\\n",
    "& \\approx & 200~\\mathrm{nHz}~\\left(\\frac{M}{10^8~M_\\odot}\\right)^{1/2} \\left(\\frac{a}{\\mathrm{mpc}}\\right)^{-3/2}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "What is the separation of the system?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3 Decay Timescale\n",
    "\n",
    "For a more general derivation, let's assume we have two masses $m_1$ and $m_2$, with total mass $M = m_1 + m_2$ and reduced mass $\\mu = m_1 m_2 /M$, in a circular orbit and separated by a semi-major axis $a$. From Peters (1964; Phys. Rev. 136, B1224), one can write down the change in the semi-major axis as\n",
    "\n",
    "$$\\frac{da}{dt} = -\\frac{64}{5} \\frac{G^3 \\mu M^2}{c^5 a^3}.$$\n",
    "\n",
    "The decay timescale is then \n",
    "\n",
    "$$t_{\\rm gw} = \\int_a^0 \\frac{dt}{da'} da' = \\int_0^a \\frac{5}{64} \\frac{c^5 a'^3}{G^3 \\mu M^2} da' = \\frac{5}{256} \\frac{c^5 a^4}{G^3 \\mu M^2}.$$\n",
    "\n",
    "We can write this expression in terms of a mass ratio $q \\equiv m_1/m_2 \\le 1$. It will be useful to substitute $M = m_1+m_2 = m_2 (q+1)$. Therefore, our mass terms can be written out as\n",
    "\n",
    "$$\\mu M^2 = m_1 m_2 M = m_1 m_2^2 (q+1) = m_2^3 q (q+1) = M^3 \\frac{q}{(q+1)^2}.$$\n",
    "\n",
    "Therefore, our decay timescale becomes\n",
    "\n",
    "$$t_{\\rm gw} = \\frac{5}{256}\\frac{c^5 a^4}{G^3 M^3} \\frac{(q+1)^2}{q}.$$\n",
    "\n",
    "Defining the dimensionless quantity $q_r = q/(q+1)^2$, we can write the scaling relation as\n",
    "\n",
    "$$t_{\\rm gw} = 4.88\\times 10^6~\\mathrm{yrs}~\\left(\\frac{a}{10^3~R_s}\\right)^4 \\left(\\frac{M}{10^8~M_\\odot}\\right)^{-3} q_r^{-1},$$\n",
    "\n",
    "where the Schwarzschild Radius $R_s = 2GM/c^2$ and we use the same fiducial mass of $10^8~M_\\odot$.\n",
    "\n",
    "Given your parameters above, how long will it take the system to merge?\n",
    "What kind of equal-mass binaries exhibit significant frequency evolution over a the course of a typical PTA experiment?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4 Host Galaxy Parameters\n",
    "\n",
    "McConnell and Ma (2013) found a variety of scaling relations between various properties of a galaxy and its black hole mass. Presumably, with a single black hole after merger, the galaxy should still fall close to these relations, and so the relations should hold before the merger as well. The relations they find are\n",
    "\n",
    "The Black Hole Mass-Stellar Velocity Dispersion ($M_{\\mathrm{BH}}$-$\\sigma$) relation:\n",
    "\n",
    "$$\\log_{10}(M_{\\mathrm{BH}}) = 8.32 + 5.64 \\log_{10}\\left(\\frac{\\sigma}{200~\\mathrm{km~s^{-1}}}\\right)$$\n",
    "\n",
    "The Black Hole Mass-Luminosity (V-band) ($M_{\\mathrm{BH}}$-$L_V$) relation:\n",
    "\n",
    "$$\\log_{10}(M_{\\mathrm{BH}}) = 9.23 + 1.11 \\log_{10}\\left(\\frac{L_V}{10^{11}~L_\\odot}\\right)$$\n",
    "\n",
    "The Black Hole Mass-Bulge Stellar Mass ($M_{\\mathrm{BH}}$-$M_{\\mathrm{bulge}}$) relation:\n",
    "\n",
    "$$\\log_{10}(M_{\\mathrm{BH}}) = 8.46 + 1.05 \\log_{10}\\left(\\frac{M_{\\mathrm{bulge}}}{10^{11}~M_\\odot}\\right)$$\n",
    "\n",
    "Given the more informative chirp mass estimate above, calculate values for the different characterisitics of the the host galaxy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## And beyond:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you are an expert, try re-running the code in the following ways and see how your answers change:\n",
    "\n",
    "* Reduce the simulated mass of the system in `observe()` to make the observed strain even weaker. How well can you do?\n",
    "* Add red noise into the `observe()` function using `LT.add_rednoise(psr,A,gamma)`. $A$ will be in GW strain amplitude units and you can vary $\\gamma$ in steepness, typically somewhere between 1 and 5. How does this affect the results? What new correlations do you see amongst the parameters? \n",
    "\n",
    "Remember: try adding more pulsars into your array if you need some additional signal boost!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some space to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}