From bc9c3b53280e3b0429a6390d90fc383be4af48da Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 22 May 2014 21:03:56 -0500 Subject: [PATCH] hod code cleaned out --- c_tools/hod/CMakeLists.txt | 29 ++- c_tools/hod/growthfactor.c | 73 ++++++ c_tools/hod/halo_bias.c | 300 ++++++++++++++++++++++ c_tools/hod/halo_concentration.c | 155 ++++++++++++ c_tools/hod/halo_mass_conversion.c | 86 +++++++ c_tools/hod/halo_mass_function.c | 338 +++++++++++++++++++++++++ c_tools/hod/jeans.c | 95 +++++++ c_tools/hod/main.c | 5 - c_tools/hod/mstar.c | 39 +++ c_tools/hod/nonlinear_power_spectrum.c | 204 +++++++++++++++ c_tools/hod/tasks.c | 39 --- c_tools/hod/tf_eisenstein_hu.c | 140 ++++++++++ c_tools/hod/transfnc.c | 70 +++++ c_tools/hod/transfunc_file.c | 59 +++++ c_tools/hod/xi_matter.c | 180 +++++++++++++ 15 files changed, 1757 insertions(+), 55 deletions(-) create mode 100644 c_tools/hod/growthfactor.c create mode 100644 c_tools/hod/halo_bias.c create mode 100644 c_tools/hod/halo_concentration.c create mode 100644 c_tools/hod/halo_mass_conversion.c create mode 100644 c_tools/hod/halo_mass_function.c create mode 100644 c_tools/hod/jeans.c create mode 100644 c_tools/hod/mstar.c create mode 100644 c_tools/hod/nonlinear_power_spectrum.c create mode 100644 c_tools/hod/tf_eisenstein_hu.c create mode 100644 c_tools/hod/transfnc.c create mode 100644 c_tools/hod/transfunc_file.c create mode 100644 c_tools/hod/xi_matter.c diff --git a/c_tools/hod/CMakeLists.txt b/c_tools/hod/CMakeLists.txt index ef85f62..9931729 100644 --- a/c_tools/hod/CMakeLists.txt +++ b/c_tools/hod/CMakeLists.txt @@ -1,17 +1,24 @@ include_directories(${CMAKE_CURRENT_BINARY_DIR}) -SET(hod_SRCS header.c main.c utility.c sigmac.c transfnc.c transfunc_file.c - nonlinear_power_spectrum.c least_squares.c hod_functions.c - xi_matter.c one_halo_rspace.c two_halo_rspace.c - input_params.c dFdx.c mstar.c halo_concentration.c growthfactor.c - halo_mass_conversion.c nfw_transform.c - pair_density.c tasks.c wp_minimization.c - kaiser_distortions.c - tf_eisenstein_hu.c jeans.c - populate_simulation.c test.c - dark_matter_statistics.c +SET(hod_SRCS header.c main.c utility.c sigmac.c + hod_functions.c + halo_bias.c + transfnc.c + nonlinear_power_spectrum.c + halo_mass_function.c + tf_eisenstein_hu.c + transfunc_file.c + halo_concentration.c + mstar.c + xi_matter.c + halo_mass_conversion.c + least_squares.c + jeans.c + growthfactor.c + input_params.c dFdx.c + tasks.c + populate_simulation.c meshlink2.c nbrsfind2.c i3tensor_2.c - mcmc.c halo_mass_function.c halo_bias.c nrutil.c qromo.c midpnt.c midinf.c polint.c splint.c spline.c zbrent.c qtrap.c trapzd.c cisi.c complex.c amoeba.c amotry.c gaussj.c powell.c linmin.c f1dim.c mnbrak.c brent.c gasdev.c diff --git a/c_tools/hod/growthfactor.c b/c_tools/hod/growthfactor.c new file mode 100644 index 0000000..2c09945 --- /dev/null +++ b/c_tools/hod/growthfactor.c @@ -0,0 +1,73 @@ +#include +#include +#include +#include "header.h" + +/* This calculates the linear growthfactor at redshift z. + * Currently, the default is a flat LCDM universe, but can easily + * be changed. + */ +double x03_g1, + xi3_g1; +double func_D0(double); +double func_Di(double); + +double growthfactor(double z) +{ + int i; + double zp1,x03,xi3,bi,b0,yy,sqx,sqx1,dbdx,x0,xi,lambda_i,htemp,omega_L,omega_i,hubble_i, + astart,fac,redshift; + + redshift=z; + astart=1.0/(z+1); + zp1=redshift+1; + omega_L=1-OMEGA_M; + + hubble_i = sqrt(OMEGA_M/(astart*astart*astart) + + (1.0-OMEGA_M-omega_L)/(astart*astart) +omega_L); + + if(omega_L>0) + { + lambda_i=omega_L/(omega_L+(1.0-OMEGA_M-omega_L)*zp1*zp1+OMEGA_M*pow(zp1,3.0)); + omega_i=OMEGA_M*pow(zp1,3.0)*lambda_i/omega_L; + } + else + { + lambda_i=0; + omega_i=OMEGA_M*zp1/(1.0+OMEGA_M*redshift); + } + + fac=astart; + + if((OMEGA_M < 0.99) && (omega_L > 0.001)) + { + x03_g1=x03=1.0/OMEGA_M-1.0; + xi3_g1=xi3=1.0/omega_i-1.0; + b0 = qromo(func_D0,0.0,1.0,midpnt); + bi = qromo(func_Di,0.0,1.0,midpnt); + b0=b0*sqrt(x03+1.0); + bi=bi*sqrt(xi3+1.0)*astart; + fac=bi/b0; + } + + + if((OMEGA_M < 0.99) && (omega_L == 0)) + { + x0=1.0/OMEGA_M-1.0; + xi=x0*astart; + b0 = 1.+3./x0+3.*sqrt(1+x0)*log(sqrt(1.+x0)-sqrt(x0))/pow(x0,1.5); + bi = 1.+3./xi+3.*sqrt(1+xi)*log(sqrt(1.+xi)-sqrt(xi))/pow(xi,1.5); + fac = bi/b0; + } + return(fac); +} + +double func_Di(double y) +{ + return(pow(1.0+xi3_g1*pow(y,1.2),-1.5)); +} + +double func_D0(double y) +{ + return(pow(1.0+x03_g1*pow(y,1.2),-1.5)); +} diff --git a/c_tools/hod/halo_bias.c b/c_tools/hod/halo_bias.c new file mode 100644 index 0000000..10d2ace --- /dev/null +++ b/c_tools/hod/halo_bias.c @@ -0,0 +1,300 @@ +#include +#include +#include + +#include "header.h" + + +/* This is the bias factor for halos. Mass (in h^{-1} M_sol) is the input. + * This is not the scale-dependent bias factor, which can be calculated later. + * + * There are many variants of the bias in this function. + * The one it is set up to use is from the Tinker etal (in prep, currently) + * bias paper calibrated from the SO halo catalogs from the Tinker et al mass functino + * paper. The parameters of the bias functions are calibrated such that they will + * work with any chosen halo overdensity. + * + * Other variants: + * + * The parameterization of b(M) is taken from Sheth Mo & Tormen (2001 MNRAS 232 1) + * but the actual parameters have been replaced by the values of Appendix A in + * Tinker, Weinberg, Zheng, & Zehavi. 2005 Apj (M/L paper) + * + * See also Seljak & Warren (2004). + * See also Sheth & Tormen (1999) + * See also Mo & White (1996) + */ + +double bias_from_file(double m, double r); + +double bias(double mass) +{ + double rm,sig,k,neff,b,logk,khi,klo,phi,plo,nu,psp,x; + static int flag=0; + static double pnorm, prev_delta=-1; + + // variables for the SO(DELTA) bias functions + static double bias_A, bias_a, bias_B, bias_b; + + /* Original SMT parameters */ + double a=0.707,bb=0.5,c=0.6; + + pnorm=SIGMA_8/sigmac(8.0); + rm=pow(3.0*mass/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + sig=pnorm*sigmac(rm); + + /* Tinker et al parameters */ + a=0.707; + bb=0.35; + c=0.8; + + /* Fitting to Mike Warren's simulations. */ + bb=0.28; + + /* Use the globel parameters. */ + a=BIAS_A; bb=BIAS_B; c=BIAS_C; + + /* First normalize the power spectrum + */ + pnorm=SIGMA_8/sigmac(8.0); + rm=pow(3.0*mass/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + sig=pnorm*sigmac(rm); + + + + /* This is from Tinker etal in prep for SO halos + */ + if(DELTA_HALO != prev_delta) + { + bias_A = pow(fabs(log10(DELTA_HALO)-2.69),2)*0.16 + 0.785; + bias_a = (log10(DELTA_HALO) - 2.28)*0.45; + bias_B = 0.4*(log10(DELTA_HALO)-2.3)+0.7*log10(DELTA_HALO)*exp(-pow(4/log10(DELTA_HALO),4.0))+1.23; + bias_b = 2.4; + fprintf(stderr,"BIAS PARAMS: %f %f %f %f\n",bias_A, bias_a, bias_B, bias_b); + prev_delta = DELTA_HALO; + } + + a = pow(sig,-bias_a); + b = 1 - bias_A*a/(a+1) + bias_B*pow(sig,-bias_b); + + return(b); + + /* This is the Seljak & Warren (2004) bias. + * There's an alpha_s in the correction term, which here I set to zero. (RUNNING.) + * See App. A of Tinker et al (M/L) for a comparison of this bias with above bias. + */ + x=mass/MSTAR; + b=(0.53+0.39*pow(x,0.45)+0.13/(40*x+1) + 5.0e-4*pow(x,1.5)); + b=b+log10(x)*(0.4*(OMEGA_M-0.3+SPECTRAL_INDX-1)+0.3*(SIGMA_8-0.9+HUBBLE-0.7)+0.8*0.0); + return(b); + + /* This is the Sheth Mo Tormen bias. + * (possible that parameter values have been changed to better fit simulations, + * ie from Tinker etal 2005 ML paper). + */ + nu=DELTA_CRIT/sig; + b=1+1.0/(sqrt(a)*DELTA_CRIT)*(sqrt(a)*a*nu*nu + sqrt(a)*bb*pow(a*nu*nu,1-c) - + (pow(a*nu*nu,c)/(pow(a*nu*nu,c)+bb*(1-c)*(1-c/2.)))); + return(b); + + + + /* This is Sheth & Tormen + */ + nu = DELTA_CRIT/sig; + nu = nu*nu; + return(1 + (0.707*nu - 1)/DELTA_CRIT + 2*0.3/DELTA_CRIT/(1+pow(0.707*nu,0.3))); + + /* This is the old Mo & White (1996) formula + */ + return(1+DELTA_CRIT/sig/sig-1/DELTA_CRIT); + + +} + +/* Just like the halo mass function, we'll set this up such that + * you just need to interpolate. + * + * Now this has been changed to calculate the spatial-scale dependence + * of the bias factor. If you don't want the scale dep. b, then just + * input r<0. + * + * If the global flag LINEAR_PSP==0, uses the scale dependence calculated for + * for halo bias relative to the non-linear matter \xi_m(r): + * + * f^2(r) = (1.0+xi*1.17)^1.49/(1.0+xi*0.69)^2.09 --> b(r) = b0*f(r) + * + * For LINEAR_PSP==1, use scale dependence determined for the linear P(k): + * + * f(r) = 1 + exp[-(r/A)^0.7] --> where A is a parameter that we're gonna have to + * determine in more detail, but now A \approx 1 + */ +double bias_interp(double m, double r) +{ + static int flag=0,prev_cosmo=0; + static double *x,*y,*y2, pnorm; + int i,n=100; + double dm,max=16.3,min=9,a,b,m1,m2,dm1,xi,power,rm,sig,b1,b2,mass,rvir,a1; + + if(!flag || RESET_COSMOLOGY!=prev_cosmo) + { + if(!ThisTask && OUTPUT) + fprintf(stdout,"RESET: resetting bias for %f %f\n",OMEGA_M,SIGMA_8); + if(!flag) + { + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + } + flag=1; + dm=(double)(max-min)/n; + for(i=1;i<=n;++i) + { + x[i]=pow(10.0,min+i*dm); + y[i]=log(bias(x[i])); + x[i] = log(x[i]); + //printf("BIAS %f %f\n",x[i]/2.302, y[i]/2.302); + continue; + + // no longer need to do this part, since we're taking into account + // halo overdensity in the bias formula. + if(DELTA_HALO!=200) + { + x[i]=log(halo_mass_conversion2(x[i],halo_c200(x[i]),200.0,DELTA_HALO)); + } + else + { + x[i]=log(x[i]); + } + } + spline(x,y,n,2.0E+30,2.0E+30,y2); + prev_cosmo=RESET_COSMOLOGY; + pnorm=SIGMA_8/sigmac(8.0); + + } + + + m=log(m); + splint(x,y,y2,n,m,&a); + a = exp(a); + + // if we're using systematic errors in an MCMC, adjust parameter a1 (amplitude) + if(USE_ERRORS) + a *= M2N.bias_amp; + + // if no scale-dependence required, return the large-scale value + if(r<0) return a; + + + /* + * SCALE DEPENDENT BIAS!!! + *---------------------------------------------------------------------- + */ + + /* first, calculate the standard Zheng-like, mass-independent scale bias + * based on the non-linear correlation function + * (re-calibrated using the SO catalogs) + */ + xi = xi_interp(r); + b = pow(1.0+xi*0.92,2.08)*pow(1.0+xi*0.74,-2.37); + if(b<0.6)b=0.6; + + /* Now the mass-dependent term. + * Where the mass-dependence comes in the form of the large-scale bias itself + */ + /* + sig = sigmac_radius_interp(r); + if(a<0) + b *= pow(1 + 0.028*pow(a,1.53)*pow(sig,1.57),2.59)/ + pow(1 + 0.253*pow(a,1.24)*pow(sig,1.71),0.397); + */ + + // if we're using systematic errors in an MCMC, adjust parameter a1 (amplitude) + if(USE_ERRORS) { + b = (b-1)*M2N.scalebias_amp + 1; + if(b<0)b=0; + } + return b*a; + +} + + +/* This is the integrand which qromo or qtrap would call + * to calculate the large-scale galaxy bias. + * The integral is a number-weighted average of the halo + * bias function, integrated over the halo mass function. + */ +double func_galaxy_bias(double m) +{ + m=exp(m); + return(dndM_interp(m)*N_avg(m)*bias_interp(m,-1.)*m); +} + +/* + * This is to get the bias from a series of files. + * + */ +double bias_from_file(double m, double r) +{ + FILE *fp; + static int nfiles=10, *n, flag =1; + static double *mass, **rad, **bias; + float x1,x2,x3; + char fname[1000]; + int i,j,i1,i2; + double b; + + if(flag) + { + n = ivector(0,nfiles-1); + mass = dvector(0,nfiles-1); + rad = dmatrix(0,nfiles-1,1,30); + bias = dmatrix(0,nfiles-1,1,30); + + fp = openfile("/home/tinker/cosmo/SDSS/zspace_analysis/SO_calibration/halofiles/sigma.dat"); + + for(i=0;imass[nfiles-1])return -1; + + for(i=1;i=m)break; + i2 = i; + i1 = i-1; + + + for(j=2;j<=n[i1];++j) + if(rad[i1][j]>r)break; + x1 = (bias[i1][j] - bias[i1][j-1])/(rad[i1][j] - rad[i1][j-1])*(r-rad[i1][j-1]) + bias[i1][j-1]; + + for(j=2;j<=n[i2];++j) + if(rad[i2][j]>r)break; + x2 = (bias[i2][j] - bias[i2][j-1])/(rad[i2][j] - rad[i2][j-1])*(r-rad[i2][j-1]) + bias[i2][j-1]; + + b = (x2 - x1)/(mass[i2] - mass[i1])*(m - mass[i1]) + x1; + + // if(r>log(4)) + // printf("%e %e %f %f %f %f %f\n",m,r,b,x2,x1,mass[i2],mass[i1]); + return b; + +} diff --git a/c_tools/hod/halo_concentration.c b/c_tools/hod/halo_concentration.c new file mode 100644 index 0000000..71f0f9b --- /dev/null +++ b/c_tools/hod/halo_concentration.c @@ -0,0 +1,155 @@ +#include +#include +#include + +#include "header.h" + +double collapse_redshift(double z); +double cvir_pnorm_g1, + cvir_sigma_g1; + +/* This calculates and tabulates the halo concentrations + * as a function of halo mass. Uses the "Bullock model", + * described in a little more detail below. + */ + +double halo_concentration(double m) +{ + static int flag=1,n,prev_cosmo=0; + static double *x,*y,*y2; + int i; + float x1,x2,cfac; + double a,dm,x3,x4; + FILE *fp; + char fname[1000]; + + if(flag || RESET_COSMOLOGY!=prev_cosmo) + { + MSTAR = mstar(); + if(OUTPUT) + fprintf(stdout,"Calc cvir with DELTA_HALO= %f\n",DELTA_HALO); + prev_cosmo=RESET_COSMOLOGY; + n=50; + if(flag) + { + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + } + cvir_pnorm_g1=SIGMA_8/sigmac(8.0); + + flag=0; + + dm=(log(HOD.M_max)-log(1.0E8))/(n-1); + for(i=1;i<=n;++i) + { + x[i]=exp((i-1)*dm+log(1.0E8)); + y[i]=cvir_model(x[i]); + x[i]=log(halo_mass_conversion(x[i],&y[i],DELTA_HALO)); + y[i]=log(y[i]); + //y[i] = log(10); + } + spline(x,y,n,1.0E+30,1.0E+30,y2); + } + cfac = 0.13*log10(m/1.0e12) + 0.125; + cfac = 1; + if(m>80*MSTAR) m=80*MSTAR; + m=log(m); + splint(x,y,y2,n,m,&a); + return(exp(a)*cfac); + +} + + +/* This is the cvir model, which should reproduce the results + * of cvir3.f, which can be obtained off James Bullock's web site. + * + * Basic idea; For a halo of mass M, find the collapse redshift + * of that halo be finding the redshift at which rms mass variance + * sigma(Mf) = DELTA_CRIT (using linear growth factor), where Mf is some + * fraction of the redshift zero halo mass. + * + * cvir = k*(1+z_coll(M*f)) + * + * Model params: + * - k = 3.12 - fitting parameter (taken from cvir3.f) + * - f = 0.001 - fraction of halo mass that collapsed at z_coll + */ + +double cvir_model(double mass) +{ + double cvir,z_coll,zbrent(),rad; + double k=3.12; + double f=0.001; + + rad=pow(3.0*f*mass/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + cvir_sigma_g1=cvir_pnorm_g1*sigmac(rad); + + if(collapse_redshift(0.0)<0)return(k); + z_coll=zbrent(collapse_redshift,0.0,200.0,1.0E-3); + cvir=k*(1+z_coll); + return(cvir); +} + +/* This is the input function to find where sig(M)*D(z)-DELTA_CRIT = 0 + * which is the redshift at which a halo of mass M collapsed. + */ +double collapse_redshift(double z) +{ + double D; + D=growthfactor(z); + return(cvir_sigma_g1*D-DELTA_CRIT); +} + +/* Some quantities are specific to an overdensity of 200 (i.e., the Jenkins mass + * function and the halo bias relation of Tinker et al. ) + * Specfically, these are actually for FOF 0.2 halos, for which 200 is the + * current best approximation. (Although Warren et al quotes 250, which is the most recent + * value.) + * + * Therefore, the halo concentrations for Delta=200 halos need to be easily + * accessible for halo mass conversion of these quantities. The values are + * tabulates here, as opposed to the above routine which tabulates the + * concentrations for a user-specified overdensity. + */ + +double halo_c200(double m) +{ + static int flag=1,n,prev_cosmo=0; + static double *x,*y,*y2; + int i; + float x1,x2; + double a,dm,x3,x4; + FILE *fp; + char fname[1000]; + + if(flag || RESET_COSMOLOGY!=prev_cosmo) + { + if(OUTPUT) + fprintf(stdout,"Calc cvir with DELTA_HALO= %f\n",200.0); + prev_cosmo=RESET_COSMOLOGY; + n=50; + if(flag) + { + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + } + cvir_pnorm_g1=SIGMA_8/sigmac(8.0); + + flag=0; + + dm=(log(HOD.M_max)-log(1.0E8))/(n-1); + for(i=1;i<=n;++i) + { + x[i]=exp((i-1)*dm+log(1.0E8)); + y[i]=cvir_model(x[i]); + x[i]=log(halo_mass_conversion(x[i],&y[i],200.0)); + y[i]=log(y[i]); + } + } + m=log(m); + splint(x,y,y2,n,m,&a); + return(exp(a)); + +} diff --git a/c_tools/hod/halo_mass_conversion.c b/c_tools/hod/halo_mass_conversion.c new file mode 100644 index 0000000..34cce65 --- /dev/null +++ b/c_tools/hod/halo_mass_conversion.c @@ -0,0 +1,86 @@ +#include +#include +#include +#include "header.h" + + +/* This is a routine which takes a virial mass (as defined by the virial + * overdensity relation (see fitting formula of Bryan & Normax) + * and converts if to a mass of some given overdensity. This is all + * taken from the appendix of Hu & Kravtsov. 2003ApJ...584..702H + */ + +/* Constants from conversion formula in Hu & Kravtsov (2003) + */ +#define a1 0.5116 +#define a2 -0.4283 +#define a3 -3.13E-3 +#define a4 -3.52E-5 + +double HK_func(double x); + +double halo_mass_conversion(double mvir, double *cvir1, double delta_halo) +{ + double x,y,z,vx,vy,vz,x1,x2,x3,x4,omega_m,delta=180,mstar,delta_vir,delta_fof,p,f, + mp,v,mass,lambda,error,tol=1.0E-3,cprev,m_delta,delta_fac,cvir,rdelta,rvir; + + cvir=*cvir1; + omega_m=OMEGA_M; + lambda=1-OMEGA_M; + mstar=MSTAR; + delta=delta_halo; + + x=omega_m-1; + delta_vir=(18*PI*PI+82*x-39*x*x)/(1+x); + + /* Now convert from virial mass to m_delta. + */ + f=delta/delta_vir*HK_func(1.0/cvir); + p=a2+a3*log(f)+a4*log(f)*log(f); + x1=1.0/sqrt(a1*pow(f,2*p)+0.5625)+2*f; + m_delta=mvir*(delta/delta_vir*pow(1.0/(x1*cvir),3.0)); + + /* Now convert the cvir to c_delta. + */ + rvir=pow(3*mvir/(4*delta_vir*PI*RHO_CRIT*OMEGA_M),1.0/3.0); + rdelta=pow(3*m_delta/(4*delta*PI*RHO_CRIT*OMEGA_M),1.0/3.0); + *cvir1=cvir*rdelta/rvir; + + return(m_delta); +} + +/* This is a slight modification to the above routine-- instead of converting from + * the virial mass, it converts from a specified Delta to the other Delta. + * (so here the "_vir" quantities are any arbitrary input overdensity.) + */ + +double halo_mass_conversion2(double mvir, double cvir1, double delta_vir, double delta_halo) +{ + double x,y,z,vx,vy,vz,x1,x2,x3,x4,omega_m,delta=180,mstar,delta_fof,p,f, + mp,v,mass,lambda,error,tol=1.0E-3,cprev,m_delta,delta_fac,cvir,rdelta,rvir; + + cvir=cvir1; + omega_m=OMEGA_M; + lambda=1-OMEGA_M; + mstar=MSTAR; + delta=delta_halo; + + /* Now convert from virial mass to m_delta. + */ + f=delta/delta_vir*HK_func(1.0/cvir); + p=a2+a3*log(f)+a4*log(f)*log(f); + x1=1.0/sqrt(a1*pow(f,2*p)+0.5625)+2*f; + m_delta=mvir*(delta/delta_vir*pow(1.0/(x1*cvir),3.0)); + + /* Now convert the cvir to c_delta. + */ + rvir=pow(3*mvir/(4*delta_vir*PI*RHO_CRIT*OMEGA_M),1.0/3.0); + rdelta=pow(3*m_delta/(4*delta*PI*RHO_CRIT*OMEGA_M),1.0/3.0); + + return(m_delta); +} + +double HK_func(double x) +{ + return(x*x*x*(log(1+1.0/x)-1.0/(1+x))); +} diff --git a/c_tools/hod/halo_mass_function.c b/c_tools/hod/halo_mass_function.c new file mode 100644 index 0000000..139bc9b --- /dev/null +++ b/c_tools/hod/halo_mass_function.c @@ -0,0 +1,338 @@ +#include +#include +#include + +#include "header.h" + +/* This is a calculation of the differential halo mass function of + * Jenkins et al 2001 (MNRAS 321, 372). + * + * Also included is the mass function from my own analysis, which is + * a variant of the Sheth-Tormen mass function (leaving all 3 parameters + * as free parameters), with a 2-parameter high-mass cutoff of the form + * suggested by Reed et al. + * + */ +void initialize_mass_function(double *a1, double *a2, double *a3, double *a4); + +double halo_mass_function(double mass) +{ + double sig,logm,a,slo,shi,rm,rlo,rhi,mlo,mhi,dsdM,n,nuprime,nufnu,p,A; + static int flag=0,SO180=0,SO324=0,WARREN=0,ST=0,JENKINS=0; + static double pnorm, prev_delta; + double btemp = -1; + + static double + a1 = 0.325277, + a2 = 0.492785, + a3 = 0.310289, + a4 = 1.317104, + a5 = 2.425681; + + /* Jenkins et al. SO 180 best-fit + */ + if(SO180) + { + JENKINS_A = 0.301; + JENKINS_B = 0.64; + JENKINS_C = 3.82; + } + if(SO324) + { + JENKINS_A = 0.316; + JENKINS_B = 0.67; + JENKINS_C = 3.82; + } + if(JENKINS) //their .2 fof function + { + JENKINS_A = 0.315; + JENKINS_B = 0.61; + JENKINS_C = 3.8; + } + + + /* First normalize the power spectrum + */ + pnorm=SIGMA_8/sigmac(8.0); + rm=pow(3.0*mass/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + sig=pnorm*sigmac(rm); + logm=log10(mass); + + mlo=0.99*mass; + mhi=1.01*mass; + rlo=pow(3.0*mlo/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + rhi=pow(3.0*mhi/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + + slo=pnorm*sigmac(rlo); + shi=pnorm*sigmac(rhi); + dsdM=(shi-slo)/(mhi-mlo); + + if(SO324)goto JENKINS_FUNCTION; + if(SO180)goto JENKINS_FUNCTION; + if(WARREN)goto WARREN_FUNCTION; + if(ST)goto ST_FUNCTION; + if(JENKINS)goto JENKINS_FUNCTION; + + /* Tinker et al. (in prep) for SO 200 + */ + if(DELTA_HALO != prev_delta) + { + initialize_mass_function(&a1,&a2,&a3,&a4); + prev_delta = DELTA_HALO; + + // if we're using systematic errors in an MCMC, adjust parameter a1 (amplitude) + if(USE_ERRORS) + a1 *= M2N.mf_amp; + + fprintf(stderr,"MF PARAMS for DELTA=%f %f %f %f %f\n",DELTA_HALO,a1,a2,a3,a4); + } + + n = -a1*(pow(sig/a3,-a2)+1)*exp(-a4/sig/sig)*OMEGA_M*RHO_CRIT/mass/sig*dsdM; + return(n); + + /* Jenkins et al. FOF .2 best-fit (unless SO180==1) + */ + JENKINS_FUNCTION: + a=-JENKINS_A*OMEGA_M*RHO_CRIT/mass/sig; + n=a*dsdM*exp(-pow(fabs(JENKINS_B-log(sig)),JENKINS_C)); + return(n); + + /* Warren et al. (calibrated only on concordance cosmology, FOF.2) + */ + WARREN_FUNCTION: + n = -0.7234*(pow(sig,-1.625)+0.2538)*exp(-1.198/sig/sig)*OMEGA_M*RHO_CRIT/mass/sig*dsdM; + return(n); + + + + /* Need to find the derivative dlog(sig)/dlog(M) + */ + mlo=0.99*logm; + mhi=1.01*logm; + rlo=pow(3.0*pow(10.0,mlo)/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + rhi=pow(3.0*pow(10.0,mhi)/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + slo=log10(pnorm*sigmac(rlo)); + shi=log10(pnorm*sigmac(rhi)); + dsdM=(shi-slo)/(mhi-mlo); + + ST_FUNCTION: + + /* This is a bunch of Sheth-Tormen stuff. + * NB! because I'm skipping the above derivative (dlogs/dlogM), i'm using the lower + */ + nuprime=0.841*DELTA_CRIT/sig; + nufnu=0.644*(1+1.0/pow(nuprime,0.6))*(sqrt(nuprime*nuprime/2/PI))*exp(-nuprime*nuprime/2); + //n=RHO_CRIT*OMEGA_M/mass*mass*nufnu*fabs(dsdM); + n=RHO_CRIT*OMEGA_M/mass*nufnu*fabs(dsdM)/sig; + return(n); + + + + +} + + +/* It may be a bit costly to run the above function every time you need + * dn/dM, so here we put the values into an array and then interpolate. + * + * The currentrange of masses calculated is 10^9 to 10^16.7. The tabulation is + * done in log(M), so the spline interpolation will perform a power-law fit + * to masses outside this range. + */ +double dndM_interp(double m) +{ + static int flag=0,prev_cosmo=0; + static double *x,*y,*y2; + int i,n=2000; + double dm,max=16.7,min=8,a,m1,m2,dm1; + + if(!flag || RESET_COSMOLOGY!=prev_cosmo) + { + if(!ThisTask && OUTPUT) + fprintf(stdout,"RESET: resetting mass function for %f %f\n",OMEGA_M,SIGMA_8); + fflush(stdout); + + if(!flag) + { + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + } + flag=1; + dm=(double)(max-min)/n; + for(i=1;i<=n;++i) + { + x[i]=pow(10.0,min+i*dm); + y[i]=log(halo_mass_function(x[i])); + x[i]=log(x[i]); + continue; + + /* + m1 = halo_mass_conversion2(x[i]*1.001,halo_concentration(x[i]*1.001),DELTA_HALO,200); + m2 = halo_mass_conversion2(x[i]*0.999,halo_concentration(x[i]*0.999),DELTA_HALO,200.0); + dm1 = (x[i]*1.001 - x[i]*0.999)/(m1-m2); + y[i] = y[i] + log(dm1); + x[i]=log(halo_mass_conversion2(x[i],halo_concentration(x[i]),DELTA_HALO,200.0)); + printf("%e %e %e %e\n",exp(x[i]),exp(y[i]),0.1*exp(y[i]),0.1); + continue; + */ + if(DELTA_HALO!=200) + { + m1 = halo_mass_conversion2(x[i]*1.001,halo_c200(x[i]*1.001),200.0,DELTA_HALO); + m2 = halo_mass_conversion2(x[i]*0.999,halo_c200(x[i]*0.999),200.0,DELTA_HALO); + dm1 = (x[i]*1.001 - x[i]*0.999)/(m1-m2); + y[i] = y[i] + log(dm1); + x[i]=log(halo_mass_conversion2(x[i],halo_c200(x[i]),200.0,DELTA_HALO)); + } + else + x[i]=log(x[i]); + + } + spline(x,y,n,2.0E+30,2.0E+30,y2); + prev_cosmo=RESET_COSMOLOGY; + } + //exit(0); + m=log(m); + splint(x,y,y2,n,m,&a); + return(exp(a)); + +} + + +/* Reads mass function in from a file and spline interpolates. + */ +double nbody_massfunction(double m) +{ + static int flag=0,n; + static double *x,*y,*y2,log10_2,normhi,normlo; + float x1,x2,x3; + int i; + double a,dx; + char aa[1000]; + FILE *fp; + + if(!flag) + { + log10_2=log10(2.0); + flag++; + if(!(fp=fopen(Files.MassFuncFile,"r"))) + endrun("ERROR opening MassFuncFile"); + i=0; + n = filesize(fp); + fprintf(stderr,"Read %d lines from [%s]\n",n,Files.MassFuncFile); + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + for(i=1;i<=n;++i) + { + fscanf(fp,"%f %f",&x1,&x2); + x[i]=log(x1); + y[i]=log(x2); + fgets(aa,1000,fp); + } + spline(x,y,n,2.0E+30,2.0E+30,y2); + fclose(fp); + fprintf(stderr,"Minimum halo mass in N-body dndM= %e\n",exp(x[1])); + + normhi = exp(y[n])/halo_mass_function(exp(x[n])); + normlo = exp(y[1])/halo_mass_function(exp(x[1])); + } + + m=log(m); + // if(m>x[n])return(0); + /*if(mx[n]) + return(halo_mass_function(exp(m))*normhi); + if(m=1600) *a1 = 0.26; + + //second parameter + y[1] = 1.466904e+00 ; + y[2] = 1.521782e+00 ; + y[3] = 1.559186e+00 ; + y[4] = 1.614585e+00 ; + y[5] = 1.869936e+00 ; + y[6] = 2.128056e+00 ; + y[7] = 2.301275e+00 ; + y[8] = 2.529241e+00 ; + y[9] = 2.661983e+00 ; + + spline(x,y,n,1.0E+30,1.0E+30,z); + splint(x,y,z,n,log(DELTA_HALO),a2); + + //third parameter + y[1] = 2.571104e+00 ; + y[2] = 2.254217e+00 ; + y[3] = 2.048674e+00 ; + y[4] = 1.869559e+00 ; + y[5] = 1.588649e+00 ; + y[6] = 1.507134e+00 ; + y[7] = 1.464374e+00 ; + y[8] = 1.436827e+00 ; + y[9] = 1.405210e+00 ; + + spline(x,y,n,1.0E+30,1.0E+30,z); + splint(x,y,z,n,log(DELTA_HALO),a3); + + + //fourth parameter + y[1] = 1.193958e+00; + y[2] = 1.270316e+00; + y[3] = 1.335191e+00; + y[4] = 1.446266e+00; + y[5] = 1.581345e+00; + y[6] = 1.795050e+00; + y[7] = 1.965613e+00; + y[8] = 2.237466e+00; + y[9] = 2.439729e+00; + + spline(x,y,n,1.0E+30,1.0E+30,z); + splint(x,y,z,n,log(DELTA_HALO),a4); + + // now adjust for redshift + if(!(REDSHIFT>0))return; + + ztemp = REDSHIFT; + if(REDSHIFT>3) ztemp = 3.0; + *a1 *= pow(1+ztemp,-0.14); + *a2 *= pow(1+ztemp,-0.14); + at = -pow(0.75/log10(DELTA_HALO/75),1.2); + at = pow(10.0,at); + *a3 *= pow(1+ztemp,-at); + +} diff --git a/c_tools/hod/jeans.c b/c_tools/hod/jeans.c new file mode 100644 index 0000000..1f4831b --- /dev/null +++ b/c_tools/hod/jeans.c @@ -0,0 +1,95 @@ +#include +#include +#include +#include +#include + +#include "header.h" + +double *mu_x,*mu_y,*mu_z; +int mu_n; + +/* Equation (10) of vdB et al, assuming alpha = 1 (standard NFW) + */ +double func_jeans1(double x) +{ + return(x/(1+x)/(1+x)); +} + +/* Equation (11) of vdB et al. + */ +double func_jeans2(double x) +{ + double y; + splint(mu_x,mu_y,mu_z,mu_n,x,&y); + return(y/(x*x*x*(1+x)*(1+x))); +} +double jeans_dispersion(double mass, double rx, double v[]) +{ + int i,j,k,n,nc,nr; + double rmin,rmax,dlogr; + static double cmin,cmax,dlogc; + double rvir,rs,s1,sig,mu1,vcirc,cvir; + + static long IDUM2 = -4555; + + static double fac=-1; + static int prev_cosmo=-1; + + if(fac<0) + { + mu_n = nc = 1000; + cmin = 0.1; + cmax = 100.0; + dlogc = (log(cmax) - log(cmin))/(nc-1); + + mu_x = dvector(1,nc); + mu_y = dvector(1,nc); + mu_z = dvector(1,nc); + + for(i=1;i<=nc;++i) + { + mu_x[i] = exp((i-1)*dlogc)*cmin; + mu_y[i] = qromo(func_jeans1,0,mu_x[i],midpnt); + } + spline(mu_x,mu_y,nc,1.0E+30,1.0E+30,mu_z); + fac=sqrt(4.499E-48)*3.09E19; + } + + cvir = halo_concentration(mass); + rvir = pow(3*mass/(4*PI*DELTA_HALO*OMEGA_M*RHO_CRIT),THIRD); + rs = rvir/cvir; + rx = rx/rs; + vcirc = fac*fac*mass/rvir; + splint(mu_x,mu_y,mu_z,mu_n,cvir,&mu1); + + s1 = qromo(func_jeans2,rx,cmax,midpnt); + sig = sqrt(cvir*vcirc/mu1*rx*(1+rx)*(1+rx)*s1); + + if(isnan(sig)) + sig = sqrt(vcirc/2); + + for(i=0;i<3;++i) + v[i]=gasdev(&IDUM2)*sig*VBIAS; + + return sig; + + nr = 100; + rmin = rvir*0.01; + rmax = rvir*0.99; + dlogr = (log(rmax) - log(rmin))/(nr-1); + rs = rvir/cvir; + + vcirc = fac*fac*mass/rvir; + splint(mu_x,mu_y,mu_z,mu_n,cvir,&mu1); + + for(i=0;i2) - { - if(atoi(argv[2])==999) - test(argc,argv); - } tasks(argc,argv); } diff --git a/c_tools/hod/mstar.c b/c_tools/hod/mstar.c new file mode 100644 index 0000000..5aced0f --- /dev/null +++ b/c_tools/hod/mstar.c @@ -0,0 +1,39 @@ +#include +#include +#include +#include "header.h" + +/* This function calculates M_star (non-linear mass scale) for the given + * cosmological paramters. + */ + +double sigma_Mmdelta_c(double lnM); +double pnorm1; + +double mstar() +{ + double sig,lnMmin,lnMmax,M_star; + + sig=sigmac(8.0); + pnorm1 = SIGMA_8/sig; + + lnMmin=log(1e7); + lnMmax=log(1e18); + M_star=zbrent(sigma_Mmdelta_c,lnMmin,lnMmax,1e-5); + M_star=exp(M_star); + if(!ThisTask) + fprintf(stderr,"M_star = %e h^{-1}M_sol\n",M_star); + return(M_star); +} + +/*** solve for M_* ***/ +double sigma_Mmdelta_c(double lnM) +{ + double sig,M,rm; + + M=exp(lnM); + rm=pow(3.0*M/(4.0*PI*OMEGA_M*RHO_CRIT),1.0/3.0); + sig=pnorm1*sigmac(rm); + + return sig-DELTA_CRIT; +} diff --git a/c_tools/hod/nonlinear_power_spectrum.c b/c_tools/hod/nonlinear_power_spectrum.c new file mode 100644 index 0000000..e08e141 --- /dev/null +++ b/c_tools/hod/nonlinear_power_spectrum.c @@ -0,0 +1,204 @@ +#include +#include +#include + +#include "header.h" + +/* This is my implementation of the Halo Model fitting function + * described in Appendix C of Smith, Peacock, et al. 2003 MNRAS 341, 1311 + * + * Comparison to Smith et al. program halofit.f output is successful. + */ + +/* Internal functions + */ +double func_nl(double r); +void tabulate_pk_nl(double *kk, double *pknl, int nk); + +/* This returns the linear power + * Delta = 4pi*k^3 P(k)/(2pi)^3 + */ +double linear_power_spectrum(double xk) +{ + static double *kk,*pknl,*y2,pnorm=-1,ahi,bhi; + static int flag=1,nk=1000,prev_cosmology=0; + double a,psp,x1[4],y1[4]; + int i; + + if(pnorm<0 || prev_cosmology!=RESET_COSMOLOGY) + { + pnorm=SIGMA_8/sigmac(8.0); + pnorm*=pnorm; + prev_cosmology=RESET_COSMOLOGY; + } + if(ITRANS>0) + psp=pow(xk,SPECTRAL_INDX)*pow(transfnc(xk),2.); + else + psp=pow(xk,SPECTRAL_INDX); + psp=psp*pnorm*xk*xk*xk/(2*PI*PI); + return(psp); +} + +double nonlinear_power_spectrum(double xk) +{ + static double *kk,*pknl,*y2,pnorm=-1,ahi,bhi; + static int flag=1,nk=1000,prev_cosmology=0; + double a,psp,x1[4],y1[4],xklog; + int i; + + if(flag || RESET_COSMOLOGY!=prev_cosmology) + { + prev_cosmology=RESET_COSMOLOGY; + pnorm=SIGMA_8/sigmac(8.0); + flag=0; + kk=dvector(1,nk); + pknl=dvector(1,nk); + y2=dvector(1,nk); + tabulate_pk_nl(kk,pknl,nk); + spline(kk,pknl,nk,1.0E+30,1.0E+30,y2); + + /* This takes the last four points in the power spectrum at high k + * and fits a power law to them for extrapolation. + */ + for(i=0;i<4;++i) + { + x1[i]=(kk[nk-3+i]); + y1[i]=(pknl[nk-3+i]); + } + least_squares(x1,y1,4,&ahi,&bhi); + + } + + xklog=log(xk); + + /* If xk is less than the smallest k value tabulates, return linear power. + */ + if(xklogkk[nk]) + return((exp((ahi+bhi*xklog)))); + + splint(kk,pknl,y2,nk,xklog,&a); + return(exp(a)); + +} +void tabulate_pk_nl(double *kk, double *pknl, int nk) +{ + double r_nl,pnorm,DeltaQ,DeltaLin,DeltaH,psp,neff,s1,s2,n1,n2,ncurve, + lnk_lo,lnk_hi,dlnk,xk1,y,xk,fy; + double an,bn,cn,f1b,f2b,f3b,alpha,beta,gamma,nu,mu; + int i,j,n=10000; + + /* First normalize the linear power spectrum. + */ + pnorm=SIGMA_8/sigmac(8.0); + + /* Calculate the non-linear scale. + */ + r_nl = exp(zbrent(func_nl,log(0.01),log(10.0),1.0E-4)); + if(OUTPUT) + fprintf(stdout,"R_NL= %f\n",r_nl); + + /* Calculate the effective spectral index at the non-linear scale. + */ + s1=pnorm*sigmac(-r_nl*0.999); + s2=pnorm*sigmac(-r_nl*1.001); + neff=-(3+2*(s2-s1)/0.002); + if(OUTPUT) + fprintf(stderr,"neff= %f\n",neff); + + /* Spectral curvature. + */ + lnk_hi=10.0; + lnk_lo=-10.0; + dlnk=(lnk_hi-lnk_lo)/n; + s1=0; + for(i=1;i<=n;++i) + { + xk1=exp(lnk_lo+(i-0.5)*dlnk); + y=xk1*r_nl; + DeltaLin=linear_power_spectrum(xk1); + s1+=DeltaLin*y*y*(1-y*y)*exp(-y*y)*dlnk; + } + ncurve=(3+neff)*(3+neff)+4*s1; + if(OUTPUT) + fprintf(stderr,"ncurve= %f\n",ncurve); + + /* Coefficients of the model. + */ + an=pow(10.0,1.4861 + 1.8369*neff + 1.6762*neff*neff + 0.7940*neff*neff*neff + + 0.1670*pow(neff,4.0) - 0.6202*ncurve); + bn=pow(10.0,0.9463 + 0.9466*neff + 0.3084*neff*neff - 0.9400*ncurve); + + cn= pow(10.0,-0.2807 + 0.6669*neff + 0.3214*neff*neff - 0.0793*ncurve); + + gamma = 0.8649 + 0.2989*neff + 0.1631*ncurve; + alpha = 1.3884 + 0.3700*neff - 0.1452*neff*neff; + beta = 0.8291 + 0.9854*neff + 0.3401*neff*neff; + mu = pow(10.0,-3.5442 + 0.1908*neff); + nu = pow(10.0, 0.9589 + 1.2857*neff); + + /* Testing the parameter dependence. + * TESTING TESTING + */ + /* + alpha *= 0.3; + beta *= 0.3; + */ + + /* Omega-dependent functions (FLAT LAMBDA COSMOLOGY) + */ + f1b = pow(OMEGA_M,-0.0307); + f2b = pow(OMEGA_M,-0.0585); + f3b = pow(OMEGA_M,+0.0743); + + + /* Tabulate the power spectrum + */ + lnk_lo=log(0.001); + lnk_hi=log(1000.0); + dlnk=(lnk_hi-lnk_lo)/(nk-1); + + for(i=1;i<=nk;++i) + { + xk=exp(lnk_lo+(i-1)*dlnk); + y=xk*r_nl; + fy=y/4.0 + y*y/8.0; + + /* TEST */ + /* fy*=1.2; */ + + DeltaLin=linear_power_spectrum(xk); + + DeltaQ = DeltaLin*pow(1+DeltaLin,beta)/(1+alpha*DeltaLin)*exp(-fy); + + DeltaH = an*pow(y,3*f1b)/(1+bn*pow(y,f2b)+pow(cn*f3b*y,3-gamma)); + + DeltaH*= 1.0/(1+mu/y+nu/(y*y)); + + kk[i]=log(xk); + pknl[i]=log(DeltaQ + DeltaH); + } +} + +/* This is a function to find the non-linear scale. Similar to M_star, + * but here we're using a Guassian smoothing window rather than top hat. + * (And we're finding the scale at which the variance = 1, not 1.686). + */ +double func_nl(double r) +{ + static int prev_cosmology=0; + double sig; + static double pnorm=-1; + + if(pnorm<0 || RESET_COSMOLOGY!=prev_cosmology) + pnorm=SIGMA_8/sigmac(8.0); + prev_cosmology=RESET_COSMOLOGY; + + r=exp(r); + sig=pnorm*sigmac(-r); + return sig-1.0; +} diff --git a/c_tools/hod/tasks.c b/c_tools/hod/tasks.c index ffc8318..7f0ac4c 100644 --- a/c_tools/hod/tasks.c +++ b/c_tools/hod/tasks.c @@ -37,18 +37,10 @@ void tasks(int argc, char **argv) - /* This is for chi^2 minimization of data for the projected correlation - * function. - */ - if(Task.wp_minimize && !HOD.color) - wp_minimization(argv[1]); - /* This is for Monte-Carlo Markov Chain exploration of the posterior * distribution of the parameters, either real-space or redshift-space, * depending on what MCMC is set to. */ - if(Task.MCMC) - mcmc_minimization(); /* This is to output the shape of the mean occupation function and the * scatter about the mean. @@ -85,26 +77,6 @@ void tasks(int argc, char **argv) * 5 - projected correlation function (1st column is now r_p) * 6 - projected correlation function without z-space correction */ - if(Task.real_space_xi || Task.All) - { - fprintf(stderr,"\n\nCALCULATING REAL-SPACE CORRELATION FUNCTION.\n"); - fprintf(stderr, "--------------------------------------------\n\n"); - - sprintf(fname,"%s.r_space",Task.root_filename); - fp=fopen(fname,"w"); - dr=(log(70.0)-log(0.01))/49.0; - for(i=0;i<50;++i) - { - r=exp(i*dr+log(0.01)); - x1=one_halo_real_space(r); - x2=two_halo_real_space(r); - x3=projected_xi(r); - x4 = projected_xi_rspace(r); - fprintf(fp,"%f %e %e %e %e %e\n",r,x1,x2,x1+x2,x3,x4); - fflush(fp); - } - fclose(fp); - } /* This takes a halofile from a simulation and populates the halos * with galaxies according the HOD specified in the batch file. @@ -123,9 +95,6 @@ void tasks(int argc, char **argv) * 2 - linear P(k) [Mpc/h]^3 * 3 - non-linear P(k) [Mpc/h]^3 */ - if(Task.matter_pk) - output_matter_power_spectrum(); - /* Output the linear and non-linear dark matter power spectrum. * Non-linear xi(r) is Fourier transform of Smith et al (above) * @@ -134,8 +103,6 @@ void tasks(int argc, char **argv) * 2 - linear xi(r) * 3 - non-linear xi(r) */ - if(Task.matter_xi) - output_matter_correlation_function(); /* Output the matter variance as a function of scale. * @@ -145,8 +112,6 @@ void tasks(int argc, char **argv) * 3 - sigma(r) [non-linear, using Smith] * 4 - mass [M_sol/h] mass = (4/3)*PI*r^3*rho_bar */ - if(Task.sigma_r) - output_matter_variance(); /* Output the halo concentrations using Bullock et al (2001) model * @@ -154,8 +119,6 @@ void tasks(int argc, char **argv) * 1 - mass [Msol/h] --> mass specified by DELTA_HALO (input file). * 2 - halo concentration. (for the specified halo definition) */ - if(Task.cvir) - output_halo_concentrations(); /* Output the halo mass function using the results of Tinker et al (2008) * @@ -163,8 +126,6 @@ void tasks(int argc, char **argv) * 1 - mass [Msol/h] * 2 - dn/dM [Mph/c]^-3 (the differential mass function). */ - if(Task.dndM) - output_halo_mass_function(); endrun("finished with tasks"); } diff --git a/c_tools/hod/tf_eisenstein_hu.c b/c_tools/hod/tf_eisenstein_hu.c new file mode 100644 index 0000000..d657322 --- /dev/null +++ b/c_tools/hod/tf_eisenstein_hu.c @@ -0,0 +1,140 @@ +#include +#include +#include +#include "header.h" + +// Transfer function of Eisenstein & Hu 1998 +// (Equation numbers refer to this paper) + +double calc_tf_eh(double k); + +double tf_eisenstein_hu(double k) +{ + static int flag=0,prev_cosmo=0; + static double *x,*y,*y2; + int n=1000,i; + double a,rlo=1.0E-4,rhi=1.0E+4,dlogr,klo; + double xi_linear_int(); + + if(!flag || RESET_COSMOLOGY!=prev_cosmo) + { + if(!flag) + { + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + } + flag=1; + + dlogr = (log(rhi)-log(rlo))/(n-1); + for(i=1;i<=n;++i) + { + x[i] = exp((i-1)*dlogr)*rlo; + y[i] = log(calc_tf_eh(x[i])); + //printf("TK %e %e %e %e %e\n",x[i],exp(y[i]),exp(y[i]),exp(y[i]),exp(y[i])); + x[i] = log(x[i]); + } + spline(x,y,n,2.0E+30,2.0E+30,y2); + prev_cosmo=RESET_COSMOLOGY; + } + + splint(x,y,y2,n,log(k),&a); + return(exp(a)); +} + +double calc_tf_eh(double k) +{ + double rk,e,thet,thetsq,thetpf,b1,b2,zd,ze,rd,re,rke,s,rks,q,y,g; + double ab,a1,a2,ac,bc,f,c1,c2,tc,bb,bn,ss,tb,tk_eh; + double h,hsq,om_mhsq,om_b,om_m; + + // set up cosmology + h = HUBBLE; + om_m = OMEGA_M; + om_b = OMEGA_B; + + // convert k to Mpc^-1 rather than hMpc^-1 + rk=k*h; + hsq=h*h; + om_mhsq=om_m*hsq; + + // constants + e=exp(1.); + thet=2.728/2.7; + thetsq=thet*thet; + thetpf=thetsq*thetsq; + + // Equation 4 - redshift of drag epoch + b1=0.313*pow(om_mhsq,-0.419)*(1.+0.607*pow(om_mhsq,0.674)); + b2=0.238*pow(om_mhsq,0.223); + zd=1291.*(1.+b1*pow(om_b*hsq,b2))*pow(om_mhsq,0.251) + /(1.+0.659*pow(om_mhsq,0.828)); + + // Equation 2 - redshift of matter-radiation equality + ze=2.50e4*om_mhsq/thetpf; + + // value of R=(ratio of baryon-photon momentum density) at drag epoch + rd=31500.*om_b*hsq/(thetpf*zd); + + // value of R=(ratio of baryon-photon momentum density) at epoch of matter-radiation equality + re=31500.*om_b*hsq/(thetpf*ze); + + // Equation 3 - scale of ptcle horizon at matter-radiation equality + rke=7.46e-2*om_mhsq/(thetsq); + + // Equation 6 - sound horizon at drag epoch + s=(2./3./rke)*sqrt(6./re)*log((sqrt(1.+rd)+sqrt(rd+re))/(1.+sqrt(re))); + + // Equation 7 - silk damping scale + rks=1.6*pow(om_b*hsq,0.52)*pow(om_mhsq,0.73)*(1.+pow(10.4*om_mhsq,-0.95)); + + // Equation 10 - define q + q=rk/13.41/rke; + + // Equations 11 - CDM transfer function fits + a1=pow(46.9*om_mhsq,0.670)*(1.+pow(32.1*om_mhsq,-0.532)); + a2=pow(12.0*om_mhsq,0.424)*(1.+pow(45.0*om_mhsq,-0.582)); + ac=pow(a1,(-om_b/om_m))*pow(a2,pow(-(om_b/om_m),3.)); + + // Equations 12 - CDM transfer function fits + b1=0.944/(1.+pow(458.*om_mhsq,-0.708)); + b2=pow(0.395*om_mhsq,-0.0266); + bc=1./(1.+b1*(pow(1.-om_b/om_m,b2)-1.)); + + // Equation 18 + f=1./(1.+pow(rk*s/5.4,4.)); + + // Equation 20 + c1=14.2 + 386./(1.+69.9*pow(q,1.08)); + c2=14.2/ac + 386./(1.+69.9*pow(q,1.08)); + + // Equation 17 - CDM transfer function + tc=f*log(e+1.8*bc*q)/(log(e+1.8*bc*q)+c1*q*q) + + (1.-f)*log(e+1.8*bc*q)/(log(e+1.8*bc*q)+c2*q*q); + + // Equation 15 + y=(1.+ze)/(1.+zd); + g=y*(-6.*sqrt(1.+y)+(2.+3.*y)*log((sqrt(1.+y)+1.)/(sqrt(1.+y)-1.))); + + // Equation 14 + ab=g*2.07*rke*s/pow(1.+rd,0.75); + + // Equation 23 + bn=8.41*pow(om_mhsq,0.435); + + // Equation 22 + ss=s/pow(1.+pow(bn/rk/s,3.),1./3.); + + // Equation 24 + bb=0.5+(om_b/om_m) + (3.-2.*om_b/om_m)*sqrt(pow(17.2*om_mhsq,2.)+1.); + + // Equations 19 & 21 + tb=log(e+1.8*q)/(log(e+1.8*q)+c1*q*q)/(1+pow(rk*s/5.2,2.)); + tb=(tb+ab*exp(-pow(rk/rks,1.4))/(1.+pow(bb/rk/s,3.)))*sin(rk*ss)/rk/ss; + + // Equation 8 + tk_eh=(om_b/om_m)*tb+(1.-om_b/om_m)*tc; + + return tk_eh; +} + diff --git a/c_tools/hod/transfnc.c b/c_tools/hod/transfnc.c new file mode 100644 index 0000000..02d6557 --- /dev/null +++ b/c_tools/hod/transfnc.c @@ -0,0 +1,70 @@ +/* PROGRAM TRANSFERFUNCTION + + --- transfnc(xk) + --- compute the transfer function T(k) of the power spectrum + --- T(k) defined as P(k) = k^{xindx}*T^{2}(k) + + * itrans=type of transfer function + 0 -> no change (returns 1) + 4 -> Efstathiou, Bond & White transfer function with Gamma as + specified (eqn. 7) + 5 -> Eisnstein & Hu + 11 -> read in TF from file (usually CMBFAST) + + NOTE: xk is in h/Mpc and is defined as k=2pi/lambda (not k=1/lambda) +*/ + +#include +#include +#include +#include "header.h" + +/* EBW CDM parameters */ + +#define aa 6.4 +#define bb 3.0 +#define cc 1.7 +#define xnu 1.13 +#define twopi 6.283185 +#define xnuinv -0.884956 + +double transfnc(double xk) +{ + double transf; + double q,t1,t2; + + if(xk==0.) + { + transf=1.; + return (double)transf; + } + + switch(ITRANS) + { + case 0: + transf=1.; + break; + + case 4: + q = xk/GAMMA; + t1=aa*q+pow((bb*q),1.5)+(cc*q)*(cc*q); + t2=pow(t1,xnu); + transf=pow((1.+t2),xnuinv); + break; + + case 5: + transf = tf_eisenstein_hu(xk); + break; + + case 11: + transf = transfunc_file(xk); + break; + + default: + fprintf(stderr,"transfnc> Unrecognized transfer function %d \n",ITRANS); + exit(-1); + break; + + } + return (double)transf; +} diff --git a/c_tools/hod/transfunc_file.c b/c_tools/hod/transfunc_file.c new file mode 100644 index 0000000..b5ccf32 --- /dev/null +++ b/c_tools/hod/transfunc_file.c @@ -0,0 +1,59 @@ +#include +#include +#include +#include "header.h" + +/* This routine reads in a transfer function from a file. + * - col 1 = k [h/Mpc] + * - col 2 = T(k) + * other columns not used. + * + * The TF is considered to be un-normalized, and will normalize all values + * by the entry in the first line. + * + * The TF is stored in arrays an spline interpolation is used at all k values. + * The interpolation is in log(k)/log(TF) to preserve the power-law dependence + * of T(k) on k outside the k-range of the file. + */ + +double transfunc_file(double xk) +{ + static double *x,*y,*y2; + static int flag=1,n; + int i; + double t,x0; + FILE *fp; + char a[1000]; + float x1,x2; + + if(flag) + { + flag=0; + + fp=openfile(Files.TF_file); + n=filesize(fp); + + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + for(i=1;i<=n;++i) + { + fscanf(fp,"%f %f",&x1,&x2); + x[i]=x1; + y[i]=x2; + if(i==1)x0=y[i]; + fgets(a,1000,fp); + y[i]/=x0; + x[i]=log(x[i]); + y[i]=log(y[i]); + } + fclose(fp); + spline(x,y,n,1.0E+30,1.0E+30,y2); + } + xk=log(xk); + splint(x,y,y2,n,xk,&t); + return(exp(t)); + +} + + diff --git a/c_tools/hod/xi_matter.c b/c_tools/hod/xi_matter.c new file mode 100644 index 0000000..2bec7b4 --- /dev/null +++ b/c_tools/hod/xi_matter.c @@ -0,0 +1,180 @@ +#include +#include +#include + +#include "header.h" + +/* This calculates and tabulates both the linear and non-linear + * matter correlation function. + * + * If the parameter BOX_SIZE is set, then the lower limits of the + * Fourier transform in 1/BOX_SIZE, else it starts at k=0. + * + * The integral (Smith etal 2003 MNRAS.341.1311S Eq [4]) to transform: + * + * Int_0^\infty Delta(k) sin(rk)/rk/k dk + * + * This integral sometimes gives qromo problems when calculating xi(r) at large r. + * Therefore, we cut off the integral at k = 10^3, which has negligible effect on the + * correlation function at the relevent scales. + */ + +double r_g4; +double xi_int(double xk); + +/* Calculates and tabulates the non-linear matter correlation function. + * Since this is only really needed for the scale-dependence of the bias, + * which is basically 1 at scales r>~8 Mpc/h, I won't calculate this much + * past that value. + */ +double xi_interp(double r) +{ + static int flag=0,prev_cosmo=0; + static double *x,*y,*y2; + int n=30,i,j; + double a,xi_int(),rhi=95,rlo=0.1,dlogr,klo,s1,s2,tolerance=1.0e-6; + + if(!flag || RESET_COSMOLOGY!=prev_cosmo) + { + if(!flag) + { + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + } + flag=1; + dlogr = (log(rhi)-log(rlo))/(n-1); + + for(i=1;i<=n;++i) + { + klo = 0; + if(BOX_SIZE>0)klo = 1/BOX_SIZE; + r_g4 = x[i] = exp((i-1)*dlogr)*rlo; + + j=1; + s1 = qromo(xi_int,klo,j/r_g4,midpnt); + s2 = s1; + klo = j/r_g4; + while(mabs(s1)>tolerance*mabs(s2)) { + j+=16; + s1 = qromo(xi_int,klo,j/r_g4,midpnt); + s2 += s1; + klo = j/r_g4; + } + y[i]=s2; + } + check_for_smoothness(x,y,n,0); + spline(x,y,n,2.0E+30,2.0E+30,y2); + prev_cosmo=RESET_COSMOLOGY; + } + + splint(x,y,y2,n,r,&a); + return(a); +} + + +/* This is the integrand of Smith et al Eq. [4] + */ +double xi_int(double xk) +{ + double xk1,xk2,psp; + + if(xk==0)return(0); + + /* power spectrum at xk + */ + psp=nonlinear_power_spectrum(xk); + + /* Integrand of Fourier transform + */ + xk1=r_g4*xk; + psp*=sin(xk1)/xk1/xk; + return psp; +} + + +double xi_linear_interp(double r) +{ + static int flag=0,prev_cosmo=0; + static double *x,*y,*y2; + int n=100,i; + double a,rlo=0.1,rhi=150,dlogr,klo; + double xi_linear_int(); + + if(!flag || RESET_COSMOLOGY!=prev_cosmo) + { + if(!flag) + { + x=dvector(1,n); + y=dvector(1,n); + y2=dvector(1,n); + } + flag=1; + + dlogr = (log(rhi)-log(rlo))/(n-1); + klo = 0; + if(BOX_SIZE>0)klo = 1/BOX_SIZE; + for(i=1;i<=n;++i) + { + r_g4 = x[i] = exp((i-1)*dlogr)*rlo; + y[i] = qromo(xi_linear_int,klo,1.0/r_g4,midpnt)+ + qromo(xi_linear_int,1.0/r_g4,1.0E+3,midpnt); + } + check_for_smoothness(x,y,n,0); + spline(x,y,n,2.0E+30,2.0E+30,y2); + prev_cosmo=RESET_COSMOLOGY; + } + + splint(x,y,y2,n,r,&a); + return(a); +} + +double xi_linear_int(xk) +double xk; +{ + double psp; + double xk1,xk2; + + /* power spectrum at xk + */ + psp=linear_power_spectrum(xk); + + /* Now take Fourier transform + */ + xk1=r_g4*xk; + psp*=sin(xk1)/xk1/xk; + + return psp; +} + +/* Sometimes one or two of the correlation function values will + * totally crap out, [a side-effect of qromo] so here we check for that and if + * we find it then we interpolate the values in log-space. + */ +void check_for_smoothness(double *x, double *y, int n, double r) +{ + int i,flag,flag2=0; + double m,b,new; + + for(i=2;i0)flag2=1; + if(y[i]<0 && !flag2)continue; + if(x[i]3.0 && fabs(y[i]/y[i+1])>3.0)flag=1; + if(fabs(y[i]/y[i-1])<0.2 && fabs(y[i]/y[i+1])<0.2)flag=1; + if(y[i]<0 && (y[i-1]>=0 && y[i+1]>=0))flag=1; + if(y[i+1]<0)flag=0; + if(!flag)continue; + + m=(log(y[i+1])-log(y[i-1]))/(log(x[i+1])-log(x[i-1])); + b=log(y[i+1])-m*log(x[i+1]); + new=m*log(x[i])+b; + + fprintf(stderr,"SMOOTHING: %e %e %e r= %f new=%e\n",y[i-1],y[i],y[i+1],x[i],exp(new)); + y[i]=exp(new); + } + + +}