- "source": "#############################################\n#\n# extract position space correlation functions\n#\n#############################################\n\n\n\ndef J0(r,nu):\n return -1.*np.sin(pi*nu/2.)*r**(-3.-1.*nu)*special.gamma(2+nu)/(2.*pi**2.)\nNmax = 256\nbk = -1.1001\nkmax = 100.\nk0 = 1.e-4\nrtab = np.zeros(Nmax)\nrmin = 0.01\nrmax = 1000.\nDelta = log(kmax/k0) / (Nmax - 1)\nDelta_r = log(rmax/rmin) / (Nmax - 1)\n\nPdiscrin0 = np.zeros(Nmax);\nPdiscrin1 = np.zeros(Nmax);\nPdiscrin2 = np.zeros(Nmax);\njsNm = np.arange(-Nmax//2,Nmax//2+1,1)\netam = bk + 2*1j*pi*(jsNm)/(1.*Nmax)/Delta\nkbins3 = np.zeros(Nmax);\n\nfor i in range(Nmax):\n kbins3[i] = k0 * exp(Delta * i)\n kinloop1 = kbins3[i] * h\n rtab[i] = rmin * exp(Delta_r * i)\n theory0 = (M1.pk_lin(kinloop1,z_pk))*h**3.\n theory1 = (M.pk(kinloop1,z_pk)[0]+M.pk(kinloop1,z_pk)[14]+2*cs*M.pk(kinloop1,z_pk)[10]/h**2.)*h**3.\n theory2 = (M1.pk(kinloop1,z_pk)[0]+M1.pk(kinloop1,z_pk)[14]+2*cs*M1.pk(kinloop1,z_pk)[10]/h**2.)*h**3.\n Pdiscrin0[i] = theory0 * exp( -1.*(kinloop1/4.)**4. -1.*bk*i*Delta)\n Pdiscrin1[i] = theory1 * exp( -1.*(kinloop1/4.)**4. -1.*bk*i*Delta)\n Pdiscrin2[i] = theory2 * exp( -1.*(kinloop1/4.)**4. -1.*bk*i*Delta)\n\ncm0 = np.fft.fft(Pdiscrin0)/ Nmax\ncm1 = np.fft.fft(Pdiscrin1)/ Nmax\ncm2 = np.fft.fft(Pdiscrin2)/ Nmax\ncmsym0 = np.zeros(Nmax+1,dtype=np.complex128)\ncmsym1 = np.zeros(Nmax+1,dtype=np.complex128)\ncmsym2 = np.zeros(Nmax+1,dtype=np.complex128)\n\nfor i in range(Nmax+1):\n if (i+2 - Nmax//2) < 1:\n cmsym0[i] = k0**(-etam[i])*np.conjugate(cm0[-i + Nmax//2])\n cmsym1[i] = k0**(-etam[i])*np.conjugate(cm1[-i + Nmax//2])\n cmsym2[i] = k0**(-etam[i])*np.conjugate(cm2[-i + Nmax//2])\n else:\n cmsym0[i] = k0**(-etam[i])* cm0[i - Nmax//2]\n cmsym1[i] = k0**(-etam[i])* cm1[i - Nmax//2]\n cmsym2[i] = k0**(-etam[i])* cm2[i - Nmax//2]\n \ncmsym0[-1] = cmsym0[-1] / 2\ncmsym0[0] = cmsym0[0] / 2\ncmsym1[-1] = cmsym1[-1] / 2\ncmsym1[0] = cmsym1[0] / 2\ncmsym2[-1] = cmsym2[-1] / 2\ncmsym2[0] = cmsym2[0] / 2\n\nxi0 = np.zeros(Nmax)\nxi1 = np.zeros(Nmax)\nxi2 = np.zeros(Nmax)\n\nfor i in range(Nmax):\n for j in range(Nmax + 1):\n xi0[i] = xi0[i] + np.real(cmsym0[j]*J0(rtab[i],etam[j]))\n xi1[i] = xi1[i] + np.real(cmsym1[j]*J0(rtab[i],etam[j]))\n xi2[i] = xi2[i] + np.real(cmsym2[j]*J0(rtab[i],etam[j]))\n\nxi0inter = interpolate.InterpolatedUnivariateSpline(rtab,xi0)\nxi1inter = interpolate.InterpolatedUnivariateSpline(rtab,xi1)\nxi2inter = interpolate.InterpolatedUnivariateSpline(rtab,xi2)\n\nrvec = np.logspace(0.,np.log10(3),1000) # array of rvec in Mpc/h\nrvec = 10**rvec\n# xi_lin = np.zeros(len(rvec))\n# xi_loop = np.zeros(len(rvec))\n# xi_loopir = np.zeros(len(rvec))\n# for i in range(len(rvec)):\n# xi_lin[i]=xi0inter(rvec[i])*rvec[i]\n# xi_loop[i]=xi1inter(rvec[i])*rvec[i]\n# xi_loopir[i]=xi2inter(rvec[i])*rvec[i]\n \nxi_lin = []\nxi_loop =[]\nxi_loopir =[]\nfor r in rvec:\n xi_lin.append(xi0inter(r)*r)\n xi_loop.append(xi1inter(r)*r)\n xi_loopir.append(xi2inter(r)*r)\n\nfig_xi, ax_xi = plt.subplots()\nax_xi.plot(rvec,xi_lin,color='purple',linestyle='-.',label='linear')\nax_xi.plot(rvec,xi_loop,color='b',linestyle='--',label='1-loop, no IR resummation')\nax_xi.plot(rvec,xi_loopir,color='red',linestyle='-',label='1-loop, IR resummation')\nax_xi.set_xlim([60.,140.])\nax_xi.set_ylim([-0.03,0.15])\nax_xi.set_xlabel(r'$r \\,\\,\\,\\, [h^{-1}\\mathrm{Mpc}]$')\nax_xi.set_ylabel(r'$r\\cdot \\xi(r)\\,\\,\\,\\, [h^{-1}\\mathrm{Mpc}]$')\nax_xi.legend(fontsize='16',ncol=1,loc='upper left')\nfig_xi.tight_layout()\nfig_xi.savefig('real_Xi.pdf')"
0 commit comments