guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 15:21:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

hod code cleaned out

Browse source
This commit is contained in:
P.M. Sutter 2014-05-22 21:03:56 -05:00
parent 1ec0d26b20
commit bc9c3b5328
15 changed files with 1757 additions and 55 deletions
Download patch file Download diff file Expand all files Collapse all files

29
c_tools/hod/CMakeLists.txt
View file

@ -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

73
c_tools/hod/growthfactor.c Normal file
View file

@ -0,0 +1,73 @@
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#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));
}

300
c_tools/hod/halo_bias.c Normal file
View file

@ -0,0 +1,300 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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;i<nfiles;++i)
fscanf(fp,"%lf %f",&mass[i],&x1);
for(i=0;i<nfiles;++i)
mass[i] = log(mass[i]);
fclose(fp);
for(i=0;i<nfiles;++i)
{
sprintf(fname,"/home/tinker/cosmo/SDSS/zspace_analysis/SO_calibration/halofiles/bias_nl.200.m%d",i);
fp = openfile(fname);
n[i] = filesize(fp);
for(j=1;j<=n[i];++j)
{
fscanf(fp,"%lf %lf %f %f %f",&rad[i][j],&bias[i][j],&x1,&x2,&x3);
rad[i][j] = log(rad[i][j]);
}
fclose(fp);
}
flag = 0;
}
r = log(r);
if(m>mass[nfiles-1])return -1;
for(i=1;i<nfiles;++i)
if(mass[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;
}

155
c_tools/hod/halo_concentration.c Normal file
View file

@ -0,0 +1,155 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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));
}

86
c_tools/hod/halo_mass_conversion.c Normal file
View file

@ -0,0 +1,86 @@
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#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)));
}

338
c_tools/hod/halo_mass_function.c Normal file
View file

@ -0,0 +1,338 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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(m<x[1])return(0);*/
if(m>x[n])
return(halo_mass_function(exp(m))*normhi);
if(m<x[1])
return(halo_mass_function(exp(m))*normlo);
splint(x,y,y2,n,m,&a);
/*if(m<x[1])printf("LESS %e %e %e\n",pow(10.0,m),pow(10.0,a),pow(10.0,y[1]));*/
return(exp(a));
}
void initialize_mass_function(double *a1, double *a2, double *a3, double *a4)
{
int n = 9, i;
double *x, *y, *z, at, ztemp;
x = dvector(1,n);
y = dvector(1,n);
z = dvector(1,n);
// initialize the overdensities
for(i=1;i<=9;i+=2)
x[i] = log(200*pow(2.0,(i-1.0)/2.0));
for(i=2;i<=9;i+=2)
x[i] = log(300*pow(2.0,(i-2.0)/2.0));
//first parameter
y[1] = 1.858659e-01 ;
y[2] = 1.995973e-01 ;
y[3] = 2.115659e-01 ;
y[4] = 2.184113e-01 ;
y[5] = 2.480968e-01 ;
y[6] = 2.546053e-01 ;
y[7] = 2.600000e-01 ;
y[8] = 2.600000e-01 ;
y[9] = 2.600000e-01 ;
spline(x,y,n,1.0E+30,1.0E+30,z);
splint(x,y,z,n,log(DELTA_HALO),a1);
if(DELTA_HALO>=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);
}

95
c_tools/hod/jeans.c Normal file
View file

@ -0,0 +1,95 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#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;i<nr;++i)
{
rx = exp(i*dlogr)*rmin/rs;
s1 = qromo(func_jeans2,rx,cmax,midpnt);
sig = cvir/mu1*rx*(1+rx)*(1+rx)*s1*2;
printf("%f %f\n",rx*rs/rvir,sqrt(sig));
}
exit(0);
}

5
c_tools/hod/main.c
View file

@ -98,11 +98,6 @@ int main(int argc, char **argv)
* arg==999 goes to the test program, superceding tasks.
* arg<0 supercedes the MCMC random number in the batfile.
*/
if(argc>2)
{
if(atoi(argv[2])==999)
test(argc,argv);
}
tasks(argc,argv);
}

39
c_tools/hod/mstar.c Normal file
View file

@ -0,0 +1,39 @@
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#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;
}

204
c_tools/hod/nonlinear_power_spectrum.c Normal file
View file

@ -0,0 +1,204 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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(xklog<kk[1])
return(linear_power_spectrum(xk));
/* If xk larger than highest k, return extrapolation.
*/
if(xklog>kk[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;
}

39
c_tools/hod/tasks.c
View file

@ -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");
}

140
c_tools/hod/tf_eisenstein_hu.c Normal file
View file

@ -0,0 +1,140 @@
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#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;
}

70
c_tools/hod/transfnc.c Normal file
View file

@ -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 <stdio.h>
#include <math.h>
#include <stdlib.h>
#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;
}

59
c_tools/hod/transfunc_file.c Normal file
View file

@ -0,0 +1,59 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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));
}

180
c_tools/hod/xi_matter.c Normal file
View file

@ -0,0 +1,180 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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;i<n;++i)
{
flag=0;
if(y[i]>0)flag2=1;
if(y[i]<0 && !flag2)continue;
if(x[i]<r)continue;
if(fabs(y[i]/y[i-1])>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);
}
}