Generate mock data for 15 parameter in LCDM cosmology.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import bilby\n",
    "from scipy import interpolate\n",
    "from matplotlib import pyplot as plt\n",
    "import random\n",
    "from scipy import integrate\n",
    "outdir='outdir'\n",
    "#######################################################################\n",
    "zmax=9\n",
    "deltaz=0.1\n",
    "c_l=3*10**5\n",
    "def EZ(x,omm):\n",
    "    return 1/np.sqrt((omm*(1+x)**3+(1-omm)))\n",
    "def trapzez(H0,omm):\n",
    "    x=np.arange(0,zmax,deltaz)\n",
    "    s=np.zeros(len(x)+1)\n",
    "    for i in range(len(x)):\n",
    "        s[i]=np.trapz(EZ(x,omm)[0:i+1],x[0:i+1]) \n",
    "    return (s[:-1])*(1+x)*(c_l/H0)\n",
    "def getz(DL,H0,omm):\n",
    "    x=np.arange(0,zmax,deltaz)\n",
    "    f=interpolate.interp1d(trapzez(H0,omm),x)\n",
    "    gez=f(DL)\n",
    "    return gez\n",
    "def PDL(DL,H0,omm):\n",
    "    z=getz(DL,H0,omm)\n",
    "    return (DL*DL)/((1+z)**3)\n",
    "def normPDL(DL,H0,omm):\n",
    "    ss=np.trapz(PDL(np.arange(0,31400,200),H0,omm),np.arange(0,31400,200))\n",
    "    return PDL(DL,H0,omm)/ss   \n",
    "######################################################\n",
    "#mmin is the min value for m1\n",
    "alpha=1.5\n",
    "mmax=40\n",
    "lamda=0.1\n",
    "mpp=35\n",
    "deltapp=1.0\n",
    "beta=-0.2\n",
    "mmin=8\n",
    "deltam=2\n",
    "def fm(m):\n",
    "    return deltam/(m)-(deltam/((m-deltam)))\n",
    "def Sm(m):\n",
    "    return 1/(np.exp(fm(m-2.901))+1)\n",
    "def ppow(m):\n",
    "    return (m**(-1*alpha))*Sm(m)\n",
    "def ppp(m):\n",
    "    return np.exp(-1*(m-mpp)**2/(2*deltapp*deltapp))*Sm(m)\n",
    "def Pm1(m):\n",
    "    return ((1-lamda)*ppow(m)+lamda*ppp(m))\n",
    "def Pm2(m2,m):\n",
    "    return ((m2/m)**beta)*Sm(m2)\n",
    "def normPm1(m):\n",
    "    ss=np.trapz(Pm1(np.arange(mmin,40.5,0.5)),np.arange(mmin,40.5,0.5))\n",
    "    return (Pm1(m)/ss)\n",
    "def normPm2(m2,m):\n",
    "    ss=np.trapz(Pm2(np.arange(mmin,40.5,0.5),40),np.arange(mmin,40.5,0.5))\n",
    "    return Pm2(m2,m)/ss\n",
    "###################################################################\n",
    "#Generate prior data (DL m1 m2) with LCDM cosmology H0=70,omm=0.3\n",
    "choiceN=10000#########events numbers\n",
    "deltadl=30\n",
    "sampdl = random.choices(np.arange(1,30000,deltadl),PDL(np.arange(1,30000,deltadl),70,0.3), k=choiceN)\n",
    "sampz=getz(sampdl,70,0.3)\n",
    "sampm1 = random.choices(np.arange(mmin,40.5,0.5),Pm1(np.arange(mmin,40.5,0.5)), k=choiceN)\n",
    "sampm1z=(np.array(sampz)+1)*np.array(sampm1)\n",
    "temsampm2=np.zeros((choiceN,1))\n",
    "for i in range(choiceN):\n",
    "    temsampm2[i,:] = random.choices(np.arange(8,(sampm1[i]+0.5),0.5),Pm2(np.arange(8,(sampm1[i]+0.5),0.5),sampm1[i]),k=1)\n",
    "sampm2=temsampm2.flatten()\n",
    "sampm2z=np.array(sampm2)*((np.array(sampz)+1))\n",
    "sampq=sampm2/sampm1\n",
    "sampdl=np.array(sampdl)\n",
    "sampm1z=np.array(sampm1z)\n",
    "sampm2z=np.array(sampm2z)\n",
    "################################################################\n",
    "#Generate data for the other 12 parameter\n",
    "label = 'samplesdata'\n",
    "bilby.core.utils.setup_logger(outdir=outdir, label=label)\n",
    "duration = 4\n",
    "sampling_frequency = 2048\n",
    "geocent_time_random=np.random.uniform(10**9,10**10,choiceN)\n",
    "a_1_prior         =bilby.gw.prior.Uniform(minimum=0, maximum=0.8, name='a_1', latex_label='$a_1$', unit=None, periodic_boundary=False)\n",
    "a_2_prior         =bilby.gw.prior.Uniform(minimum=0, maximum=0.8, name='a_2', latex_label='$a_2$', unit=None, periodic_boundary=False)\n",
    "tilt_1_prior      =bilby.core.prior.Sine(name='tilt_1', latex_label='$\\\\theta_1$', unit=None, minimum=0, maximum=3.141592653589793, periodic_boundary=False)\n",
    "tilt_2_prior      =bilby.core.prior.Sine(name='tilt_1', latex_label='$\\\\theta_1$', unit=None, minimum=0, maximum=3.141592653589793, periodic_boundary=False)\n",
    "phi_12_prior      =bilby.gw.prior.Uniform(minimum=0, maximum=6.283185307179586, name='phi_12', latex_label='$\\\\Delta\\\\phi$', unit=None, periodic_boundary=True)\n",
    "phi_jl_prior      =bilby.gw.prior.Uniform(minimum=0, maximum=6.283185307179586, name='phi_jl', latex_label='$\\\\phi_{JL}$', unit=None, periodic_boundary=True)\n",
    "theta_jn_prior    =bilby.core.prior.Sine(name='theta_jn', latex_label='$\\\\theta_{JN}$', unit=None, minimum=0, maximum=3.141592653589793, periodic_boundary=False)\n",
    "psi_prior         =bilby.gw.prior.Uniform(minimum=0, maximum=3.141592653589793, name='psi', latex_label='$\\\\psi$', unit=None, periodic_boundary=True)\n",
    "phase_prior       =bilby.gw.prior.Uniform(minimum=0, maximum=6.283185307179586, name='phase', latex_label='$\\\\phi$', unit=None, periodic_boundary=True)\n",
    "geocent_time_prior=bilby.gw.prior.Uniform(minimum=geocent_time_random- 1,maximum=geocent_time_random + 1,name='geocent_time', latex_label='$t_c$', unit='$s$')\n",
    "ra_prior          =bilby.gw.prior.Uniform(minimum=0, maximum=6.283185307179586, name='ra', latex_label='$\\\\mathrm{RA}$', unit=None, periodic_boundary=True)\n",
    "dec_prior         =bilby.core.prior.Cosine(name='dec', latex_label='$\\\\mathrm{DEC}$', unit=None, minimum=-1.5707963267948966, maximum=1.5707963267948966, periodic_boundary=False)\n",
    "\n",
    "kk=choiceN#number of event\n",
    "inpa=np.zeros((kk,15))\n",
    "smaples_a_1          = a_1_prior.sample(kk)\n",
    "smaples_a_2          = a_2_prior.sample(kk)\n",
    "smaples_tilt_1       = tilt_1_prior.sample(kk)\n",
    "smaples_tilt_2       = tilt_2_prior.sample(kk)\n",
    "smaples_phi_12       = phi_12_prior.sample(kk)\n",
    "smaples_phi_jl       = phi_jl_prior.sample(kk)\n",
    "smaples_theta_jn     = theta_jn_prior.sample(kk)\n",
    "smaples_psi          = psi_prior.sample(kk)\n",
    "smaples_phase        = phase_prior.sample(kk)\n",
    "smaples_geocent_time = geocent_time_prior.sample(kk)\n",
    "smaples_ra           = ra_prior.sample(kk)\n",
    "smaples_dec          = dec_prior.sample(kk)\n",
    "\n",
    "inpa[:,0] =sampm1z#######samples_mass_1\n",
    "inpa[:,1] =sampm2z#######smaples_mass_2\n",
    "inpa[:,2] =smaples_a_1\n",
    "inpa[:,3] =smaples_a_2\n",
    "inpa[:,4] =smaples_tilt_1\n",
    "inpa[:,5] =smaples_tilt_2\n",
    "inpa[:,6] =smaples_phi_12\n",
    "inpa[:,7] =smaples_phi_jl\n",
    "inpa[:,8] =sampdl########smaples_luminosity_distance\n",
    "inpa[:,9] =smaples_theta_jn\n",
    "inpa[:,10]=smaples_psi\n",
    "inpa[:,11]=smaples_phase\n",
    "inpa[:,12]=smaples_geocent_time\n",
    "inpa[:,13]=smaples_ra\n",
    "inpa[:,14]=smaples_dec\n",
    "np.savetxt('LCDMcosmology15p.txt',inpa)"
   ]
  }
 ],
 "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
}