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

first pass at cleaning HOD

Browse source
This commit is contained in:
P.M. Sutter 2014-05-22 21:45:03 -04:00
parent 134d9b1df6
commit 57715fe0de
24 changed files with 0 additions and 5963 deletions
Download patch file Download diff file Expand all files Collapse all files

161
c_tools/hod/dark_matter_statistics.c
View file

@ -1,161 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "header.h"
void output_matter_power_spectrum()
{
int k;
double dlogk,xk,plin,pnl;
FILE *fp;
char aa[100];
fprintf(stderr,"\n\nCALCULATING MATTER POWER SPECTRA.\n");
fprintf(stderr, "---------------------------------\n\n");
sprintf(aa,"%s.matter_pk",Task.root_filename);
fp = fopen(aa,"w");
for(k=-300;k<=100;++k)
{
xk=pow(10.0,k/100.0);
plin = linear_power_spectrum(xk)/(xk*xk*xk)*TWOPI*PI;
pnl = nonlinear_power_spectrum(xk)/(xk*xk*xk)*TWOPI*PI;
fprintf(fp,"%e %e %e\n",xk,plin,pnl);
}
fclose(fp);
}
void output_matter_correlation_function()
{
int k,nr=50;
double dlogr,r,xlin,xnl,rmin = 0.05,rmax = 80.0;
FILE *fp;
char aa[100];
fprintf(stderr,"\n\nCALCULATING MATTER CORRELATION FUNCTIONIONS.\n");
fprintf(stderr, "--------------------------------------------\n\n");
sprintf(aa,"%s.matter_xi",Task.root_filename);
fp = fopen(aa,"w");
dlogr = (log(rmax) - log(rmin))/(nr-1);
for(k=0;k<nr;++k)
{
r = exp(k*dlogr)*rmin;
xlin = xi_linear_interp(r);
xnl = xi_interp(r);
fprintf(fp,"%e %e %e\n",r,xlin,xnl);
}
fclose(fp);
}
void output_halo_correlation_function(double mass)
{
int k,nr=50;
double dlogr,r,xlin,xnl,rmin = 0.15,rmax = 40.0;
FILE *fp;
char aa[100];
fprintf(stderr,"\n\nCALCULATING HALO CORRELATION FUNCTIONION (M=%.2e)\n",mass);
fprintf(stderr, "-------------------------------------------------\n\n");
sprintf(aa,"%s.halo_xi",Task.root_filename);
fp = fopen(aa,"w");
rmin = pow(3*mass/(4*DELTA_HALO*PI*RHO_CRIT*OMEGA_M),THIRD);
dlogr = (log(rmax) - log(rmin))/(nr-1);
for(k=0;k<nr;++k)
{
r = exp(k*dlogr)*rmin;
xlin = bias_interp(mass,r);
xnl = xi_interp(r);
fprintf(fp,"%e %e\n",r,xnl*xlin*xlin);
}
fclose(fp);
}
void output_matter_variance()
{
int k,nr=50;
double dlogr,r,slin,snl,mass,rmin = 0.05,rmax = 80.0,pnorm,pnorm_nl;
FILE *fp;
char aa[100];
fprintf(stderr,"\n\nCALCULATING MATTER VARIANCE.\n");
fprintf(stderr, "----------------------------\n\n");
sprintf(aa,"%s.sigma_r",Task.root_filename);
fp = fopen(aa,"w");
dlogr = (log(rmax) - log(rmin))/(nr-1);
pnorm = SIGMA_8/sigmac(8.0);
pnorm_nl = pnorm*sigmac(80.0)/nonlinear_sigmac(80.0);
for(k=0;k<nr;++k)
{
r = exp(k*dlogr)*rmin;
mass = 4./3.*PI*r*r*r*RHO_CRIT*OMEGA_M;
slin = sigmac(r)*pnorm;
snl = nonlinear_sigmac(r)*pnorm_nl;
fprintf(fp,"%e %e %e %e\n",r,slin,snl,mass);
}
fclose(fp);
}
void output_halo_concentrations()
{
int k,nr=100;
double x,dlogm,m,mvir,cdelta,mmin=1e8,mmax=1e16,delta_vir;
FILE *fp;
char aa[100];
fprintf(stderr,"\n\nCALCULATING HALO CONCENTRATIONS.\n");
fprintf(stderr, "--------------------------------\n\n");
sprintf(aa,"%s.cvir",Task.root_filename);
fp = fopen(aa,"w");
dlogm = (log(mmax) - log(mmin))/(nr-1);
x=OMEGA_M-1;
delta_vir=(18*PI*PI+82*x-39*x*x)/(1+x);
for(k=0;k<nr;++k)
{
m = exp(k*dlogm)*mmin;
cdelta = halo_concentration(m);
fprintf(fp,"%e %e\n",m,cdelta);
}
fclose(fp);
}
void output_halo_mass_function()
{
int k,nr=100;
double x,dlogm,m,mvir,cdelta,mmin=1e8,mmax=1e16,delta_vir;
FILE *fp;
char aa[100];
fprintf(stderr,"\n\nCALCULATING HALO MASS FUNCTION.\n");
fprintf(stderr, "-------------------------------\n\n");
sprintf(aa,"%s.massfunc",Task.root_filename);
fp = fopen(aa,"w");
dlogm = (log(mmax) - log(mmin))/(nr-1);
for(k=0;k<nr;++k)
{
m = exp(k*dlogm)*mmin;
x = dndM_interp(m);
fprintf(fp,"%e %e\n",m,x);
}
fclose(fp);
}

73
c_tools/hod/growthfactor.c
View file

@ -1,73 +0,0 @@
#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
View file

@ -1,300 +0,0 @@
#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;
}

178
c_tools/hod/halo_bias_error.c
View file

@ -1,178 +0,0 @@
#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.
*
* 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 (submitted)
*/
double bias(double mass)
{
static int prev_cosmo=-1;
static long IDUM=-3333;
double rm,sig,k,neff,b,logk,khi,klo,phi,plo,nu,psp,x,fac,fac1,fac2;
static int flag=0,na=4;
static double pnorm,*xa,*ya,*za;
int i;
/* Original SMT parameters */
double a=0.707,bb=0.5,c=0.6;
/* 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;
if(!flag)
{
xa=dvector(1,na);
ya=dvector(1,na);
za=dvector(1,na);
flag=1;
xa[1] = log(1/2.51);
xa[2] = log(1/1.49);
xa[3] = log(1/0.905);
xa[4] = log(1/0.501);
}
if(prev_cosmo != RESET_COSMOLOGY)
{
prev_cosmo=RESET_COSMOLOGY;
for(i=1;i<=4;++i)
ya[i] = gasdev(&IDUM);
spline(xa,ya,na,1.0E+30,1.0E+30,za);
}
/* 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);
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.))));
splint(xa,ya,za,na,log(1/sig),&fac1);
fac2 = 0.0158*pow(sig,-2.7) + 0.00251*pow(sig,2.2);
fac = 1 + fac1*fac2;
if(fac<0.01)fac=0.01;
b*=fac;
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 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 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;
int i,n=100;
double dm,max=16.3,min=9,a,b,m1,m2,dm1,xi;
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]));
if(DELTA_HALO!=200)
{
m1 = halo_mass_conversion2(x[i]*1.001,halo_c200(x[i]*1.001),DELTA_HALO,200.0);
m2 = halo_mass_conversion2(x[i]*0.999,halo_c200(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_c200(x[i]),DELTA_HALO,200.0));
}
else
x[i]=log(x[i]);
}
spline(x,y,n,2.0E+30,2.0E+30,y2);
prev_cosmo=RESET_COSMOLOGY;
}
m=log(m);
splint(x,y,y2,n,m,&a);
if(r<0 || r>10)return(exp(a));
if(LINEAR_PSP)
{
b=exp(-pow(r,0.7)) + 1;
}
else
{
xi=xi_interp(r);
b=pow(1.0+xi*1.17,1.49*0.5)*pow(1.0+xi*0.69,-2.09*0.5);
}
a=exp(a)*b;
return(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)
{
return(dndM_interp(m)*N_avg(m)*bias_interp(m,-1.));
}

155
c_tools/hod/halo_concentration.c
View file

@ -1,155 +0,0 @@
#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
View file

@ -1,86 +0,0 @@
#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
View file

@ -1,338 +0,0 @@
#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);
}

146
c_tools/hod/halo_mass_function_error.c
View file

@ -1,146 +0,0 @@
#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.
*
*/
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, prev_cosmo=-1,na=4;
static double pnorm,*xa,*ya,*za;
double fac1,fac2,fac;
static long IDUM=-555;
int i;
static double
a1 = 0.325277,
a2 = 0.492785,
a3 = 0.310289,
a4 = 1.317104,
a5 = 2.425681;
if(!flag)
{
xa=dvector(1,na);
ya=dvector(1,na);
za=dvector(1,na);
flag=1;
xa[1] = log(1/2.51);
xa[2] = log(1/1.49);
xa[3] = log(1/0.905);
xa[4] = log(1/0.501);
}
if(prev_cosmo != RESET_COSMOLOGY)
{
prev_cosmo=RESET_COSMOLOGY;
for(i=1;i<=4;++i)
ya[i] = gasdev(&IDUM);
spline(xa,ya,na,1.0E+30,1.0E+30,za);
}
/* 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(BEST_FIT)
{
n = a1*sqrt(2*a2/PI)*(1+pow(sig*sig/(a2*DELTA_CRIT*DELTA_CRIT),a3))
*DELTA_CRIT/sig*exp(-a2*DELTA_CRIT*DELTA_CRIT/(2*sig*sig))
*exp(-a4/sig/pow(fabs(cosh(2*sig)),a5));
splint(xa,ya,za,na,log(1/sig),&fac1);
fac2 = 0.00562*pow(sig,-2.9) + 0.00158*pow(sig,2.2);
fac = 1 + fac1*fac2;
if(fac<0.01)fac=0.01;
n = n*fac;
/* printf("%e %e %e %e\n",sig,n,n/fac,fac); */
n = -n*OMEGA_M*RHO_CRIT/mass/sig*dsdM;
return(n);
}
/* Jenkins et al.
*/
a=-JENKINS_A*OMEGA_M*RHO_CRIT/mass/sig;
n=a*dsdM*exp(-pow(fabs(JENKINS_B-log(sig)),JENKINS_C));
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.
*/
double dndM_interp(double m)
{
static int flag=0,prev_cosmo=0;
static double *x,*y,*y2;
int i,n=1000;
double dm,max=16.7,min=9,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]));
if(DELTA_HALO!=200)
{
m1 = halo_mass_conversion2(x[i]*1.001,halo_c200(x[i]*1.001),DELTA_HALO,200.0);
m2 = halo_mass_conversion2(x[i]*0.999,halo_c200(x[i]*0.999),DELTA_HALO,200.0);
dm1 = (x[i]*1.001 - x[i]*0.999)/(m1-m2);
y[i] = y[i] + log10(dm1);
x[i]=log10(halo_mass_conversion2(x[i],halo_c200(x[i]),DELTA_HALO,200.0));
}
else
x[i]=log(x[i]);
}
spline(x,y,n,2.0E+30,2.0E+30,y2);
prev_cosmo=RESET_COSMOLOGY;
}
m=log(m);
splint(x,y,y2,n,m,&a);
return(exp(a));
}

95
c_tools/hod/jeans.c
View file

@ -1,95 +0,0 @@
#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);
}

432
c_tools/hod/kaiser_distortions.c
View file

@ -1,432 +0,0 @@
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "header.h"
/* Local variables for the coordinate in redshift space.
*/
double rs_g8,rp_g8;
/* Internal functions.
*/
double func_kaiser(double z);
double xi_bar(double r);
double f_xibar(double r);
double xi_prime(double r, double mu);
double xi_t(double r);
double Lpoly(double x, int i);
double f_xi2bar(double r);
double xi_harmonic(double r, int i);
double scale_dependent_dispersion(double r);
double scale_dependent_skewness(double r);
void initial_dispersion_model_values(double *a, double **pp, double *yy);
void dispersion_model_input(void);
double chi2_dispersion_model(double *a);
/* Internal variables.
*/
struct DISPERION_MODEL_STRUCT {
int n;
double **rs,**rp;
double **x,**e;
} XX;
double *rr_g8, *aa_g8;
int nparams_g8, niter_g8=0;
void fit_dispersion_model()
{
int n,niter,i,j;
double *a,**pp,*yy,FTOL=1.0E-3,chi2min,s1,dlogm,m;
FILE *fp;
char aa[1000];
fprintf(stderr,"\n\nCHI2 MINIMIZATION OF xi(SIGMA,PI) WITH DISPERSION MODEL.\n");
fprintf(stderr, "--------------------------------------------------------\n\n");
if(POWELL)
FTOL=1.0E-4;
else
FTOL=1.0E-4;
nparams_g8 = n = 8;
dispersion_model_input();
wp.ncf=n;
a=dvector(1,n);
rr_g8 = dvector(1,n);
aa_g8 = dvector(1,n);
for(i=1;i<=n;++i)
rr_g8[i] = (log(10)-log(0.1))/(n-1.)*(i-1) + log(0.1);
BETA = 0.5;
if(POWELL)
pp=dmatrix(1,n,1,n);
else
pp=dmatrix(1,n+1,1,n);
yy=dvector(1,n+1);
initial_dispersion_model_values(a,pp,yy);
if(POWELL)
{
if(OUTPUT)printf("dispersion_model> starting powell.\n");
powell(a,pp,n,FTOL,&niter,&chi2min,chi2_dispersion_model);
chi2min = chi2_dispersion_model(a);
}
else
{
if(OUTPUT)printf("dispersion_model> starting amoeba.\n");
amoeba(pp,yy,n,FTOL,chi2_dispersion_model,&niter);
for(i=1;i<=n;++i)a[i]=pp[1][i];
chi2min = chi2_dispersion_model(a);
}
}
void initial_dispersion_model_values(double *a, double **pp, double *yy)
{
static int flag=0;
int i,j;
double d[100];
if(!flag) {
for(i=1;i<=nparams_g8;++i)
a[i] = 500;
printf("INITIAL VALUES: ");
for(i=1;i<=wp.ncf;++i)printf("%e ",a[i]);
printf("\n");
}
flag++;
/* Make the starting stepsize 10% of the initial values.
*/
for(i=1;i<=wp.ncf;++i)
d[i]=a[i]*0.1/flag;
if(POWELL)
{
for(i=1;i<=wp.ncf;++i)
{
for(j=1;j<=wp.ncf;++j)
{
pp[i][j]=0;
if(i==j)pp[i][j]+=d[j];
}
}
}
else
{
for(j=1;j<=wp.ncf;++j)
pp[1][j]=a[j];
yy[1]=chi2_dispersion_model(a);
for(i=1;i<=wp.ncf;++i)
{
a[i]+=d[i];
if(i>1)a[i-1]-=d[i-1];
yy[i+1]=chi2_dispersion_model(a);
for(j=1;j<=wp.ncf;++j)
pp[i+1][j]=a[j];
}
a[wp.ncf]-=d[wp.ncf];
}
}
void dispersion_model_input()
{
FILE *fp;
int i,j;
fp = openfile("xi2d.data");
XX.n = (int)sqrt(filesize(fp)*0.99)+1;
XX.rs = dmatrix(1,XX.n,1,XX.n);
XX.rp = dmatrix(1,XX.n,1,XX.n);
XX.x = dmatrix(1,XX.n,1,XX.n);
XX.e = dmatrix(1,XX.n,1,XX.n);
for(i=1;i<=XX.n;++i)
for(j=1;j<=XX.n;++j)
fscanf(fp,"%lf %lf %lf %lf",&XX.rs[i][j],&XX.rp[i][j],&XX.x[i][j],&XX.e[i][j]);
fclose(fp);
}
double chi2_dispersion_model(double *a)
{
int i,j;
double xx,r,chi2=0;
for(i=1;i<=nparams_g8;++i)
{
aa_g8[i] = a[i];
if(a[i]<=0)return(1.0E7);
}
for(i=2;i<=XX.n;++i)
for(j=2;j<=XX.n;++j)
{
xx = kaiser_distortion(XX.rs[i][j],XX.rp[i][j]);
/* printf("%f %f %f %f %f\n",XX.rs[i][j],XX.rp[i][j],XX.x[i][j],xx,xx); */
chi2 += (xx-XX.x[i][j])*(xx-XX.x[i][j])/(XX.e[i][j]*XX.e[i][j]);
}
niter_g8++;
printf("ITER %d %e ",niter_g8,chi2);
for(i=1;i<=nparams_g8;++i)
printf("%e ",a[i]);
printf("\n");
fflush(stdout);
if(niter_g8%100==0)
{
for(i=1;i<=100;++i)
{
r = exp((log(40)-log(0.1))/(100-1.)*(i-1) + log(0.1));
xx = scale_dependent_dispersion(r);
printf("PVD %d %f %f\n",niter_g8,r,xx);
}
fflush(stdout);
}
return(chi2);
}
double kaiser_distortion(double rs, double rp)
{
double s1,rlim=60,rr;
static int prev=-1;
if(prev==-1 || prev!=RESET_COSMOLOGY)
{
prev=RESET_COSMOLOGY;
/*
BETA = pow(OMEGA_M,0.6)/qromo(func_galaxy_bias,HOD.M_low,HOD.M_max,midpnt)*
GALAXY_DENSITY;
*/
printf("kaiser> BETA= %f SIGV= %f\n",BETA,SIGV);
}
rs_g8=rs;
rp_g8=rp;
rr=sqrt(rs*rs+rp*rp);
/* if(rr*1.5>10)rlim=rr*5.5; */
/* s1=qromo(func_kaiser,-rlim,rlim,midpnt); */
s1=qtrap(func_kaiser,-rlim,rlim,1.0E-3);
return(s1-1);
}
double linear_kaiser_distortion(double r, double z)
{
double mu;
mu= sqrt(1 - z/r*z/r);
mu = z/r;
return xi_prime(r,mu);
}
double func_kaiser(double z)
{
double r,mu,v,x,sigv;
r=sqrt(z*z+rs_g8*rs_g8);
mu=z/r;
v=(rp_g8-z)*100.0;
sigv = scale_dependent_dispersion(r);
/*
if(v*z<0)
sigv*=scale_dependent_skewness(r);
*/
/*
x=(1+one_halo_real_space(r)+two_halo_real_space(r))*
exp(-ROOT2*fabs(v)/sigv)/ROOT2/sigv*100.0;
*/
x=(1+xi_prime(r,mu)+one_halo_real_space(r))*exp(-ROOT2*fabs(v)/sigv)/ROOT2/sigv*100.0;
return(x);
}
double scale_dependent_skewness(double r)
{
return 1.5;
}
double scale_dependent_dispersion(double r)
{
static int niter=-1,flag=0;
static double *z,*x,*y;
double a;
int i;
if(niter!=niter_g8)
{
niter = niter_g8;
if(!flag)
{
flag++;
x=dvector(1,nparams_g8+1);
y=dvector(1,nparams_g8+1);
z=dvector(1,nparams_g8+1);
y[1] = 0;
x[1] = log(0.01);
for(i=2;i<=nparams_g8+1;++i)
x[i] = rr_g8[i-1];
}
for(i=2;i<=nparams_g8+1;++i)
y[i] = aa_g8[i-1];
spline(x,y,nparams_g8+1,1.0E+30,1.0E+30,z);
}
r=log(r);
if(r>rr_g8[nparams_g8])return(aa_g8[nparams_g8]);
splint(x,y,z,nparams_g8+1,r,&a);
return(a);
}
/* This function tabulates the spherically averaged correlation function.
* Equation (A7) in Madgwick et al.
*/
double xi_bar(double r)
{
static int flag=0,prev=-1;
static double *x,*y,*y2;
int n=100,i;
static double rmin,rmax,lgr,a;
if(!flag || prev!=RESET_KAISER)
{
if(OUTPUT>1)
printf("Tabulating xi_bar...\n");
prev=RESET_KAISER;
if(!flag)
{
x=malloc(n*sizeof(double));
x--;
y=malloc(n*sizeof(double));
y--;
y2=malloc(n*sizeof(double));
y2--;
}
flag=1;
rmin=-1.9;
rmin=log10(R_MIN_2HALO);
rmax=log10(80.0);
for(i=1;i<=n;++i)
{
lgr=(rmax-rmin)*(i-1.0)/(n-1.0)+rmin;
x[i]=pow(10.0,lgr);
y[i]=qtrap(f_xibar,0.0,x[i],1.0E-4)*3/(x[i]*x[i]*x[i]);
//printf("XIBAR %f %e\n",x[i],y[i]);
}
spline(x,y,n,1.0E30,1.0E30,y2);
}
if(r<x[1])return(-1);
splint(x,y,y2,n,r,&a);
//if(a<xi_t(r))return(xi_t(r));
return(a);
}
double f_xibar(double r)
{
double x;
x=two_halo_real_space(r);
/* x=one_halo_real_space(r)+two_halo_real_space(r); */
return(x*r*r);
}
/* This function tabulates the xi(r)*r^4 dr function.
* Equation (A8)
*/
double xi_2bar(double r)
{
static int flag=0,prev=-1;
static double *x,*y,*y2,rmin;
int n=100,i;
double rmax,lgr,a;
if(!flag || prev!=RESET_KAISER)
{
if(OUTPUT>1)
printf("Tabulating xi_2bar...\n");
prev=RESET_KAISER;
if(!flag)
{
x=malloc(n*sizeof(double));
x--;
y=malloc(n*sizeof(double));
y--;
y2=malloc(n*sizeof(double));
y2--;
}
flag=1;
rmin=-1.9;
rmin = log10(R_MIN_2HALO);
rmax=log10(80.0);
for(i=1;i<=n;++i)
{
lgr=(rmax-rmin)*(i-1.0)/(n-1.0)+rmin;
x[i]=pow(10.0,lgr);
y[i]=qtrap(f_xi2bar,0.0,x[i],1.0E-4)*5/(x[i]*x[i]*x[i]*x[i]*x[i]);
//printf("XI2BAR %f %e\n",x[i],y[i]);
}
spline(x,y,n,1.0E30,1.0E30,y2);
}
if(r<x[1])return(-1);
splint(x,y,y2,n,r,&a);
// if(a<xi_t(r))return(xi_t(r));
return(a);
}
double f_xi2bar(double r)
{
double x;
x=two_halo_real_space(r);
/* x=one_halo_real_space(r)+two_halo_real_space(r); */
return(x*r*r*r*r);
}
/* These are Legendre polynomials
*/
double Lpoly(double x, int i)
{
switch(i) {
case 0: return(1.0);
case 2: return(0.5*(3*x*x-1));
case 4: return((35*x*x*x*x-30*x*x+3)/8.0);
/*case 4: return((384*x*x*x*x+1152*x*x*(x*x-1)+144*(x*x-1)*(x*x-1))/(16.*24.));*/
default: return(0.0);
}
}
double xi_harmonic(double r, int i)
{
switch(i){
case 0: return((1+2.*BETA/3.+BETA*BETA/5.)*xi_t(r));
case 2: return((4.*BETA/3.+4.*BETA*BETA/7.)*(xi_t(r)-xi_bar(r)));
case 4: return(8.*BETA*BETA/35.*(xi_t(r)+2.5*xi_bar(r)-3.5*xi_2bar(r)));
default: return(0);
}
}
double xi_prime(double r, double mu)
{
double x0,x2,x4;
x0=xi_harmonic(r,0)*Lpoly(mu,0);
x2=xi_harmonic(r,2)*Lpoly(mu,2);
x4=xi_harmonic(r,4)*Lpoly(mu,4);
return(x0+x2+x4);
}
double xi_t(double r)
{
return(two_halo_real_space(r));
return(one_halo_real_space(r)+two_halo_real_space(r));
}

469
c_tools/hod/mcmc.c
View file

@ -1,469 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#ifdef PARALLEL
#include <mpi.h>
#endif
#include "header.h"
/* External functions from wp_minimization.c
*/
void wp_input(void);
double mcmc_initialize(double *a, double **cov1, double *avg1, double *start_dev);
/* Internal functions.
*/
double chi2_wp_wrapper(double *a);
int USE_IWEIGHT = 0;
/******************************************************************
*
* HOD.free[] also controls which variables will be held constant/vary
* during MCMC minimization. Since this routine will also so z-space
* minimization if requested, indices>6 are cosmological.
*
* i variable
* --- --------
* [1] -> M_min
* [2] -> M1
* [3] -> alpha
* [4] -> M_cut
* [5] -> sigmaM
* [6] -> CVIR_FAC
* [7] -> MaxCen (or M_cen_max)
* [8] -> M_sat_break
* [9] -> alpha1
*
* [10]-> OMEGA_M
* [11]-> SIGMA_8
* [12]-> VBIAS
* [13]-> VBIAS_C
* [14]-> GAMMA
* [15]-> SPECTRAL_INDX
*
* [0] -> The galaxy_density will be considered data with errors on it,
* and therefore no variable will be fixed by the galaxy density.
*
*/
void mcmc_minimization()
{
double stepfac=1;
double error=1,tolerance=0,**cov1,**tmp,*a,*avg1,chi2,chi2prev,
**evect,*eval,*aprev,*atemp,**tmp1,*opar,x1,fsat,**chain,*start_dev,*eval_prev;
int n,i,j,k,nrot,niter=0,count=0,imax_chain=100000,NSTEP=50,NSTEP_MAX=10000,convergence=0;
long IDUM=-555;
int *pcheck,pcnt,ptot=20,firstflag=1,*iweight,total_weight,icvir;
double t0,tprev,temp,chi2a,chi2b;
FILE *fpmcmc;
char fname[1000];
sprintf(fname,"%s.MCMC",Task.root_filename);
fpmcmc = fopen(fname,"w");
opar=dvector(1,100);
MCMC=Task.MCMC;
pcheck=calloc(ptot,sizeof(int));
wp_input();
Work.imodel=2;
Work.chi2=1;
OUTPUT=0;
srand48(32498793);
/* Find the number of free parameters in the minimization
* for the real-space correlation function.
*/
for(n=0,i=1;i<100;++i)
{
n+=HOD.free[i];
if(OUTPUT)
printf("mcmc_min> free[%i] = %d\n",i,HOD.free[i]);
}
wp.ncf=n;
/* Find out which free parameter is for CVIR_FAC
*/
j=0;
if(HOD.free[6])
for(i=0;i<6;++i)
if(HOD.free[i])j++;
icvir=j+1;
if(HOD.free[0])
{
wp.ngal = GALAXY_DENSITY;
wp.ngal_err = 0.1*wp.ngal;
FIX_PARAM = 0;
}
if(OUTPUT)
printf("mcmc_min> %d free parameters\n",n);
a=dvector(1,n);
start_dev=dvector(1,n);
aprev=dvector(1,n);
atemp=dvector(1,n);
cov1=dmatrix(1,n,1,n);
avg1=dvector(1,n);
tmp=dmatrix(1,n,1,n);
tmp1=dmatrix(1,n,1,1);
evect=dmatrix(1,n,1,n);
eval=dvector(1,n);
eval_prev=dvector(1,n);
chain=dmatrix(1,imax_chain,1,n);
iweight = ivector(1,imax_chain);
for(i=1;i<=imax_chain;++i)
iweight[i] = 0;
IDUM=IDUM_MCMC;
chi2prev=mcmc_initialize(a,cov1,avg1,start_dev);
niter++;
for(i=1;i<=n;++i)
{
aprev[i] = a[i];
chain[1][i] = a[i];
}
pcnt=0;
pcheck[pcnt]=1;
stepfac=1;
while(niter<NSTEP)
{
pcnt++;
if(pcnt==ptot)
{
for(j=i=0;i<ptot;++i)j+=pcheck[i];
stepfac = stepfac*pow(0.9,5-j);
if(!ThisTask)printf("STEPFAC %f %d %d\n",stepfac,j,count);
pcnt=0;
}
/* stepfac=0.7; */
for(i=1;i<=n;++i)
a[i] = (1+gasdev(&IDUM)*start_dev[i]*stepfac)*aprev[i];
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i])OMEGA_M = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])SPECTRAL_INDX = a[++j];
}
if(VBIAS_C<0)continue;
/* Hard-wire CVIR variation
*/
if(HOD.free[6])
CVIR_FAC = a[icvir];
chi2=chi2_wp_wrapper(a);
pcheck[pcnt]=1;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
{
if(USE_IWEIGHT)
iweight[niter+1]++;
pcheck[pcnt]=0;
continue;
}
niter++;
iweight[niter]++;
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask){
fprintf(fpmcmc,"%d %d ",niter,count);
for(i=1;i<=n;++i)
printf("%e ",a[i]);
printf("%e\n",chi2);fflush(fpmcmc);
}
}
stepfac=1.6/sqrt(n);
pcnt=-1;
t0 = second();
NSTEP = niter;
while(niter<imax_chain)
{
stepfac=1.6/sqrt(n);
if(convergence)goto SKIP_MATRIX;
for(j=1;j<=n;++j)
{
avg1[j]=0;
for(k=1;k<=n;++k)
cov1[j][k]=0;
}
total_weight = 0;
for(i=1;i<=niter;++i)
{
for(j=1;j<=n;++j)
{
avg1[j]+=chain[i][j]*iweight[i];
for(k=1;k<=n;++k)
cov1[j][k]+=chain[i][j]*chain[i][k]*iweight[i];
}
total_weight+=iweight[i];
}
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
tmp[i][j] = cov1[i][j]/total_weight - avg1[i]*avg1[j]/(total_weight*total_weight);
jacobi(tmp,n,eval,evect,&nrot);
gaussj(evect,n,tmp1,1);
SKIP_MATRIX:
for(i=1;i<=n;++i)
atemp[i] = gasdev(&IDUM)*sqrt(eval[i])*stepfac;
for(i=1;i<=n;++i)
for(a[i]=0,j=1;j<=n;++j)
a[i] += atemp[j]*evect[j][i];
for(i=1;i<=n;++i)
a[i] += aprev[i];
/* We seem to be having a problem with this.
* So, broadcast the model params from the root processor.
*/
#ifdef PARALLEL
MPI_Bcast(&a[1],n,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD);
#endif
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i])OMEGA_M = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])SPECTRAL_INDX = a[++j];
/* if(HOD.free[++i])SIGV = a[++j]; */
}
if(VBIAS_C<0)continue;
/* Hard-wire CVIR variation
*/
if(HOD.free[6])
CVIR_FAC = a[icvir];
chi2=chi2_wp_wrapper(a);
tprev = t0;
t0 = second();
++count;
pcheck[pcnt]=0;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
{
if(USE_IWEIGHT)
iweight[niter+1]++;
continue;
}
pcheck[pcnt]=1;
niter++;
if(!convergence)NSTEP = niter;
iweight[niter]++;
if(niter%NSTEP_MAX==0 && !convergence && niter>NSTEP_MAX)
{
convergence = 1;
for(i=1;i<=n;++i)
{
x1=fabs(eval[i]-eval_prev[i])/eval_prev[i];
if(x1>0.01)convergence = 0;
printf("CONVERGENCE CHECK %d %d %e %e %e\n",niter/NSTEP_MAX,i,x1,eval[i],eval_prev[i]);
}
for(i=1;i<=n;++i)
eval_prev[i] = eval[i];
convergence = 0;
if(convergence)
printf("CONVERGENCE ACCOMPLISHED %d %d \n",niter,count);
}
if(niter==NSTEP_MAX)
{
for(i=1;i<=n;++i)
eval_prev[i] = eval[i];
}
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask) {
fprintf(fpmcmc,"%d %d ",niter,count);
for(i=1;i<=n;++i)
printf("%e ",a[i]);
printf("%e\n",chi2);fflush(fpmcmc);
}
}
}
double chi2_wp_wrapper(double *a)
{
static int flag=1;
static double *b;
int i,j;
if(flag)
{
b=dvector(1,100);
flag=0;
}
for(j=0,i=1;i<=N_HOD_PARAMS;++i) {
if(HOD.free[i] && i!=5) {
if(a[++j]<=0) { printf("NEG %d %d %e\n",i,j,a[j]); return(1.0E7); } }
if(HOD.free[i] && i==5) {
++j; }
}
i=0;j=0;
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_min */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M1 */
if(HOD.free[++i]){j++;b[j]=a[j];} /* alpha */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_cut */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* sigma_logM */
if(HOD.free[++i]){j++;b[j]=a[j];} /* cvir_fac */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* MaxCen */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_sat_break */
if(HOD.free[++i]){j++;b[j]=a[j];} /* alpha1 */
return(chi2_wp(b));
}
double mcmc_initialize(double *a, double **cov1, double *avg1, double *start_dev)
{
int i,j=0;
double x1,x2,omega_m;
long IDUM = -556;
omega_m = 1;
if(MCMC>1)
omega_m = OMEGA_M;
i=0;j=0;
if(HOD.free[++i]){ a[++j]=log10(HOD.M_min/omega_m);start_dev[j]=0.001; }
if(HOD.free[++i]){ a[++j]=log10(HOD.M1/omega_m);start_dev[j]=0.001; } //.0005
if(HOD.free[++i]){ a[++j]=HOD.alpha;start_dev[j]=0.03; } //.005
if(HOD.free[++i]){ a[++j]=log10(HOD.M_cut/omega_m);start_dev[j]=0.01; } //.001
if(HOD.free[++i]){ a[++j]=log10(HOD.sigma_logM);start_dev[j]=0.01; }
if(HOD.free[++i]){ a[++j]=CVIR_FAC;start_dev[j]=0.02; }
if(HOD.pdfc==7) {
if(HOD.free[++i])a[++j]=log10(HOD.M_cen_max/omega_m); start_dev[j]=0.001; }
else {
if(HOD.free[++i])a[++j]=HOD.MaxCen; start_dev[j]=0.02; }
if(HOD.free[++i]){ a[++j]=log10(HOD.M_sat_break/omega_m);start_dev[j]=0.001; }
if(HOD.free[++i]){ a[++j]=HOD.alpha1;start_dev[j]=0.02; }
if(MCMC>1)
{
if(HOD.free[++i])a[++j]=OMEGA_M;
if(HOD.free[++i])a[++j]=SIGMA_8;
if(HOD.free[++i])a[++j]=VBIAS;
if(HOD.free[++i])a[++j]=VBIAS_C;
if(HOD.free[++i])a[++j]=GAMMA;
if(HOD.free[++i])a[++j]=SPECTRAL_INDX;
}
if(!ThisTask)
{
printf("INITIAL VALUES: ");
for(i=1;i<=wp.ncf;++i)printf("%e ",a[i]);
printf("\n");
}
for(i=1;i<=wp.ncf;++i)
{
avg1[i]=a[i];
for(j=1;j<=wp.ncf;++j)
cov1[i][j]=a[i]*a[j];
}
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i]){ OMEGA_M = a[++j]; start_dev[j] = 0.01; }
if(HOD.free[++i]){ SIGMA_8 = a[++j]; start_dev[j] = 0.01; }
if(HOD.free[++i]){ VBIAS = a[++j]; start_dev[j] = 0.01; }
if(HOD.free[++i]){ VBIAS_C = a[++j]; start_dev[j] = 0.02; }
if(HOD.free[++i]){ GAMMA = a[++j]; start_dev[j] = 0.015; }
if(HOD.free[++i]){ SPECTRAL_INDX = a[++j]; start_dev[j] = 0.02; }
}
x1=chi2_wp_wrapper(a);
x2=0;
if(!ThisTask) {
printf("TRY 0 ");
for(i=1;i<=wp.ncf;++i)
printf("%.4e ",a[i]);
printf("%e\n",x1+x2);fflush(stdout);
printf("INITIAL CHI2: %e %e\n",x1,x2);
fflush(stdout);
}
return(x1+x2);
}

539
c_tools/hod/mcmc_color.c
View file

@ -1,539 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#ifdef PARALLEL
#include <mpi.h>
#endif
#include "header.h"
/* External functions from wp_minimization.c
*/
void wp_input(void);
void wp_color_input(void);
double chi2_wp_color(double *a);
double chi2_wp_color_wrapper(double *a);
/* inteernal functions
*/
double mcmc_color_initialize(double *a, double **cov1, double *avg1, double *start_dev);
void mcmc_color_minimization()
{
double stepfac=1;
double error=1,tolerance=0,**cov1,**tmp,*a,*avg1,chi2,chi2prev,
**evect,*eval,*aprev,*atemp,**tmp1,*opar,x1,fsat,**chain,*start_dev,*eval_prev;
int n,i,j,k,nrot,niter=0,count=0,imax_chain=100000,NSTEP=50,NSTEP_MAX=10000,convergence=0;
long IDUM=-555;
int *pcheck,pcnt,ptot=20,firstflag=1,*iweight,total_weight;
double t0,tprev,temp,chi2a,chi2b;
opar=dvector(1,100);
MCMC=Task.MCMC;
pcheck=calloc(ptot,sizeof(int));
Work.imodel=2;
Work.chi2=1;
OUTPUT=1;
fprintf(stderr,"\n\nMCMC OF W_P(R_P) COLOR DATA..........\n");
fprintf(stderr, "--------------------------------------------\n\n");
HOD.blue_fraction = 0.5857; /* <- Millenium fraction (sloan);; SDSS fraction -> 0.565; */
HOD.blue_fraction = 0.6555; /* <- Millenium fraction (B-V>0.8) */
HOD.blue_fraction = 0.565; /* SDSS -19,-20 */
HOD.blue_fraction = 0.492; /* SDSS -20,-21 */
HOD.blue_fraction = 0.379; /* SDSS -21 */
wp_color.ON = 1;
wp_color_input();
srand48(32498793);
/* Find the number of free parameters in the minimization
* for the real-space correlation function.
*/
for(n=0,i=1;i<12;++i)
{
n+=HOD.free[i];
/* if(i>N_HOD_PARAMS && HOD.free[i])MCMC=3;*/
if(OUTPUT)
printf("mcmc_min> free[%i] = %d\n",i,HOD.free[i]);
}
// add 3 parameters for the color stuff
n+=3;
wp.ncf=n;
if(HOD.free[0])
{
wp.ngal = GALAXY_DENSITY;
wp.ngal_err = 0.1*wp.ngal;
FIX_PARAM = 0;
}
if(OUTPUT)
printf("mcmc_min> %d free parameters\n",n);
a=dvector(1,n);
start_dev=dvector(1,n);
aprev=dvector(1,n);
atemp=dvector(1,n);
cov1=dmatrix(1,n,1,n);
avg1=dvector(1,n);
tmp=dmatrix(1,n,1,n);
tmp1=dmatrix(1,n,1,1);
evect=dmatrix(1,n,1,n);
eval=dvector(1,n);
eval_prev=dvector(1,n);
chain=dmatrix(1,imax_chain,1,n);
iweight = ivector(1,imax_chain);
for(i=1;i<=imax_chain;++i)
iweight[i] = 0;
IDUM=IDUM_MCMC;
//chi2prev=mcmc_initialize(a,cov1,avg1,start_dev);
chi2prev=mcmc_color_initialize(a,cov1,avg1,start_dev);
//initial_color_values(a,pp,yy);
niter++;
for(i=1;i<=n;++i)
{
aprev[i] = a[i];
chain[1][i] = a[i];
}
pcnt=0;
pcheck[pcnt]=1;
stepfac=1;
while(niter<NSTEP)
{
pcnt++;
if(pcnt==ptot)
{
for(j=i=0;i<ptot;++i)j+=pcheck[i];
stepfac = stepfac*pow(0.9,5-j);
if(!ThisTask)printf("STEPFAC %f %d %d\n",stepfac,j,count);
pcnt=0;
}
/* stepfac=0.7; */
for(i=1;i<=n;++i)
a[i] = (1+gasdev(&IDUM)*start_dev[i]*stepfac)*aprev[i];
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i])OMEGA_M = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])SPECTRAL_INDX = a[++j];
/* if(HOD.free[++i])SIGV = a[++j]; */
}
if(VBIAS_C<0)continue;
/* Hard-wire CVIR variation
*/
if(HOD.free[6])
CVIR_FAC = a[3];
chi2=chi2_wp_color_wrapper(a);
if(!ThisTask){
printf("TRY %d ",++count);
for(i=1;i<=n;++i)
printf("%.4e ",a[i]);
printf("%e\n",chi2);fflush(stdout);
}
pcheck[pcnt]=1;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
{
/* This for loop puts the prev element in the chain is
* the current trial point is rejected.
*/
/* For the initialization, don't use this: we need
* separate elements for estimating the covariance matrix.
*/
/*
for(i=1;i<=n;++i)
a[i] = aprev[i];
chi2 = chi2prev;
*/
pcheck[pcnt]=0;
continue;
}
niter++;
iweight[niter]++;
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask){
printf("ACCEPT %d %d ",niter,count);
for(i=1;i<=n;++i)
printf("%e ",a[i]);
printf("%e\n",chi2);fflush(stdout);
printf("HSTATS %d %e %e %e %e\n",niter,HOD.M_min,number_weighted_halo_mass(),
number_weighted_central_mass(),
qromo(func_satellite_density,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY);
printf("FSAT %d %e %e %e %e\n",niter,HOD.M_min,wp.fsat_red,wp.fsat_blue,wp.fsat_all);
}
}
RESTART_POINT:
stepfac=1.6/sqrt(n);
pcnt=-1;
t0 = second();
NSTEP = niter;
while(niter<imax_chain)
{
pcnt++;
if(pcnt==ptot)
{
for(j=i=0;i<ptot;++i)j+=pcheck[i];
//stepfac=1.6/sqrt(n);
//stepfac=0.6/sqrt(n);
// stepfac = stepfac*pow(0.9,6-j);
stepfac = 0.25;
stepfac = 0.5;
stepfac=1.6/sqrt(n);
if(!ThisTask)printf("STEPFAC %f %d %d\n",stepfac,j,count);
pcnt=0;
}
stepfac=1.6/sqrt(n);
//stepfac = 0;
if(convergence)goto SKIP_MATRIX;
// if(niter>NSTEP_MAX && niter%NSTEP_MAX!=0)goto SKIP_MATRIX;
for(j=1;j<=n;++j)
{
avg1[j]=0;
for(k=1;k<=n;++k)
cov1[j][k]=0;
}
total_weight = 0;
for(i=1;i<=niter;++i)
{
for(j=1;j<=n;++j)
{
avg1[j]+=chain[i][j]*iweight[i];
for(k=1;k<=n;++k)
cov1[j][k]+=chain[i][j]*chain[i][k]*iweight[i];
}
total_weight+=iweight[i];
}
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
tmp[i][j] = cov1[i][j]/total_weight - avg1[i]*avg1[j]/(total_weight*total_weight);
jacobi(tmp,n,eval,evect,&nrot);
gaussj(evect,n,tmp1,1);
SKIP_MATRIX:
if(RESTART==4)convergence = 1;
if(RESTART && count==0)stepfac=0;
for(i=1;i<=n;++i)
atemp[i] = gasdev(&IDUM)*sqrt(eval[i])*stepfac;
for(i=1;i<=n;++i)
for(a[i]=0,j=1;j<=n;++j)
a[i] += atemp[j]*evect[j][i];
for(i=1;i<=n;++i)
a[i] += aprev[i];
/* We seem to be having a problem with this.
* So, broadcast the model params from the root processor.
*/
#ifdef PARALLEL
MPI_Bcast(&a[1],n,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD);
#endif
// CHeck that the chi2 for the last point in the restart chain
// is recovered.
/*
if(RESTART && !count)
for(i=1;i<=n;++i)
a[i] = aprev[i];
*/
if(RESTART==5)
{
j = count+1;
if(j>4000)exit(0);
for(i=1;i<=n;++i)
a[i] = chain[j][i];
}
/*
if(!ThisTask)
for(i=1;i<=n;++i)
{
printf("COV %d %d %e ",count,i,sqrt(eval[i]));
for(j=1;j<=n;++j)
printf("%e ",evect[j][i]);
printf("\n");
}
*/
/* Using only variances
*/
//for(i=1;i<=n;++i)
// a[i] = aprev[i] + gasdev(&IDUM)*sqrt(tmp[i][i])*stepfac;
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i])OMEGA_M = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])SPECTRAL_INDX = a[++j];
/* if(HOD.free[++i])SIGV = a[++j]; */
}
if(VBIAS_C<0)continue;
/* Hard-wire CVIR variation
*/
if(HOD.free[6])
CVIR_FAC = a[3];
chi2=chi2_wp_color_wrapper(a);
tprev = t0;
t0 = second();
++count;
if(!ThisTask) {
printf("TRY %d ",count);
for(i=1;i<=n;++i)
printf("%.4e ",a[i]);
if(RESTART==2) {
printf("%e %e %.2f\n",chi2,chi2/(1+exp(-count/100.0)),
timediff(tprev,t0));fflush(stdout); }
else {
printf("%e %.2f\n",chi2,
timediff(tprev,t0));fflush(stdout); }
}
if(0) {
printf("CPU%02d %d ",ThisTask,count);
for(i=1;i<=n;++i)
printf("%.4e ",a[i]);
if(RESTART==2) {
printf("%e %e %.2f\n",chi2,chi2/(1+exp(-count/100.0)),
timediff(tprev,t0));fflush(stdout); }
else {
printf("%e %.2f\n",chi2,
timediff(tprev,t0));fflush(stdout); }
}
pcheck[pcnt]=0;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
{
/*
for(i=1;i<=n;++i)
a[i] = aprev[i];
chi2 = chi2prev;
*/
continue;
}
pcheck[pcnt]=1;
// if(NSTEP<NSTEP_MAX)NSTEP++;
niter++;
if(!convergence)NSTEP = niter;
iweight[niter]++;
if(niter%NSTEP_MAX==0 && !convergence && niter>NSTEP_MAX)
{
convergence = 1;
for(i=1;i<=n;++i)
{
x1=fabs(eval[i]-eval_prev[i])/eval_prev[i];
if(x1>0.01)convergence = 0;
printf("CONVERGENCE CHECK %d %d %e %e %e\n",niter/NSTEP_MAX,i,x1,eval[i],eval_prev[i]);
}
for(i=1;i<=n;++i)
eval_prev[i] = eval[i];
convergence = 0;
if(convergence)
printf("CONVERGENCE ACCOMPLISHED %d %d \n",niter,count);
}
if(niter==NSTEP_MAX)
{
for(i=1;i<=n;++i)
eval_prev[i] = eval[i];
}
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask) {
printf("ACCEPT %d %d ",niter,count);
for(i=1;i<=n;++i)
printf("%e ",a[i]);
printf("%e\n",chi2);fflush(stdout);
printf("FSAT %d %e %e %e %e\n",niter,HOD.M_min,wp.fsat_red,wp.fsat_blue,wp.fsat_all);
if(MCMC==1)
{
printf("HSTATS %d %e %e %e %e\n",niter,HOD.M_min,number_weighted_halo_mass(),
number_weighted_central_mass(),
qromo(func_satellite_density,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY);
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
}
}
}
}
double chi2_wp_color_wrapper(double *a)
{
static int flag=1;
static double *b;
int i,j;
if(flag)
{
b=dvector(1,100);
flag=0;
}
for(j=0,i=1;i<=N_HOD_PARAMS;++i) {
if(HOD.free[i] && i!=5) {
if(a[++j]<=0) { printf("NEG %d %d %e\n",i,j,a[j]); return(1.0E7); } }
if(HOD.free[i] && i==5) {
++j; }
}
i=0;j=0;
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_min */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M1 */
if(HOD.free[++i]){j++;b[j]=a[j];} /* alpha */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_cut */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* sigma_logM */
if(HOD.free[++i]){j++;b[j]=a[j];} /* cvir_fac */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* MaxCen */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_sat_break */
if(HOD.free[++i]){j++;b[j]=a[j];} /* alpha1 */
// now the 3 color params
++j;
b[j] = a[j];
++j;
b[j] = a[j];
++j;
b[j] = a[j];
muh(1);
return(chi2_wp_color(b));
}
double mcmc_color_initialize(double *a, double **cov1, double *avg1, double *start_dev)
{
int i,j=0;
double x1,x2,omega_m;
long IDUM = -556;
i=0;j=0;
if(HOD.free[++i]){ a[++j]=log10(HOD.M_min); start_dev[j]=0.001; }
if(HOD.free[++i]){ a[++j]=log10(HOD.M1); start_dev[j] = 0.001; }
if(HOD.free[++i]){ a[++j]=HOD.alpha; start_dev[j] = 0.03; }
if(HOD.free[++i]){ a[++j]=log10(HOD.M_cut); start_dev[j] = 0.01; }
if(HOD.free[++i]){ a[++j]=log10(HOD.sigma_logM); start_dev[j] = 0.01; }
if(HOD.free[++i]){ a[++j]=CVIR_FAC; start_dev[j] = 0.02; }
if(HOD.pdfc>=7){
if(HOD.free[++i]){ a[++j]=log10(HOD.M_cen_max); start_dev[j] = 0.001; }}
else {
if(HOD.free[++i]) { a[++j]=log10(HOD.MaxCen); start_dev[j] = 0.001; }}
a[++j] = HOD.fblue0_sat;
start_dev[j] = 0.01;
a[++j] = HOD.sigma_fblue_sat;
start_dev[j] = 0.01;
a[++j] = HOD.sigma_fblue_cen;
start_dev[j] = 0.01;
if(!ThisTask)
{
printf("INITIAL VALUES: ");
for(i=1;i<=wp.ncf;++i)printf("%e ",a[i]);
printf("\n");
}
for(i=1;i<=wp.ncf;++i)
{
avg1[i]=a[i];
for(j=1;j<=wp.ncf;++j)
cov1[i][j]=a[i]*a[j];
}
x1=chi2_wp_color_wrapper(a);
x2 = 0;
if(!ThisTask) {
printf("TRY 0 ");
for(i=1;i<=wp.ncf;++i)
printf("%.4e ",a[i]);
printf("%e\n",x1+x2);fflush(stdout);
printf("INITIAL CHI2: %e %e\n",x1,x2);
fflush(stdout);
}
return(x1+x2);
}

632
c_tools/hod/mcmc_exp.c
View file

@ -1,632 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#ifdef PARALLEL
#include <mpi.h>
#endif
#include "header.h"
/* External functions from wp_minimization.c
*/
void wp_input(void);
double mcmc_initialize(double *a, double **cov1, double *avg1, double *start_dev);
/* Internal functions.
*/
double chi2_wp_wrapper(double *a);
void mcmc_restart2(double *start_dev, int np);
int mcmc_restart3(double **chain, int n, double *chi2_prev, int *iweight);
int USE_IWEIGHT = 0;
/******************************************************************
*
* HOD.free[] also controls which variables will be held constant/vary
* during MCMC minimization. Since this routine will also so z-space
* minimization if requested, indices>6 are cosmological.
*
* i variable
* --- --------
* [1] -> M_min
* [2] -> M1
* [3] -> alpha
* [4] -> M_cut
* [5] -> sigmaM
* [6] -> CVIR_FAC
* [7] -> MaxCen (or M_cen_max)
* [8] -> M_sat_break
* [9] -> alpha1
*
* [10]-> OMEGA_M
* [11]-> SIGMA_8
* [12]-> VBIAS
* [13]-> VBIAS_C
* [14]-> GAMMA
* [15]-> SPECTRAL_INDX
*
* [0] -> The galaxy_density will be considered data with errors on it,
* and therefore no variable will be fixed by the galaxy density.
*
*/
void mcmc_minimization()
{
double stepfac=1;
double error=1,tolerance=0,**cov1,**tmp,*a,*avg1,chi2,chi2prev,
**evect,*eval,*aprev,*atemp,**tmp1,*opar,x1,fsat,**chain,*start_dev,*eval_prev;
int n,i,j,k,nrot,niter=0,count=0,imax_chain=100000,NSTEP=150,NSTEP_MAX=1000,convergence=0;
long IDUM=-555;
int *pcheck,pcnt,ptot=20,firstflag=1,*iweight,total_weight;
double t0,tprev,temp;
opar=dvector(1,100);
MCMC=Task.MCMC;
pcheck=calloc(ptot,sizeof(int));
if(MCMC>1 && !COVAR)
wp.esys=0.05;
if(!ThisTask)printf("ESYS %f %d\n",wp.esys,MCMC);
wp_input();
Work.imodel=2;
Work.chi2=1;
/*
OUTPUT=0;
*/
srand48(32498793);
/* Find the number of free parameters in the minimization
* for the real-space correlation function.
*/
for(n=0,i=1;i<100;++i)
{
n+=HOD.free[i];
/* if(i>N_HOD_PARAMS && HOD.free[i])MCMC=3;*/
if(OUTPUT)
printf("mcmc_min> free[%i] = %d\n",i,HOD.free[i]);
}
wp.ncf=n;
if(HOD.free[0])
{
wp.ngal = GALAXY_DENSITY;
wp.ngal_err = 0.1*wp.ngal;
FIX_PARAM = 0;
}
if(OUTPUT)
printf("mcmc_min> %d free parameters\n",n);
a=dvector(1,n);
start_dev=dvector(1,n);
aprev=dvector(1,n);
atemp=dvector(1,n);
cov1=dmatrix(1,n,1,n);
avg1=dvector(1,n);
tmp=dmatrix(1,n,1,n);
tmp1=dmatrix(1,n,1,1);
evect=dmatrix(1,n,1,n);
eval=dvector(1,n);
eval_prev=dvector(1,n);
chain=dmatrix(1,imax_chain,1,n);
iweight = ivector(1,imax_chain);
for(i=1;i<=imax_chain;++i)
iweight[i] = 0;
IDUM=IDUM_MCMC;
if(RESTART)
{
niter = mcmc_restart3(chain,n,&chi2prev,iweight);
if(niter < NSTEP)
{
if(ThisTask==0)
fprintf(stderr,"Not enough points in restart chain: %d<=%d\n",niter,NSTEP);
exit(0);
}
for(i=1;i<=n;++i)
aprev[i] = chain[niter][i];
goto RESTART_POINT;
}
chi2prev=mcmc_initialize(a,cov1,avg1,start_dev);
niter++;
for(i=1;i<=n;++i)
{
aprev[i] = a[i];
chain[1][i] = a[i];
}
pcnt=0;
pcheck[pcnt]=1;
if(RESTART==2)
{
mcmc_restart2(start_dev,n);
}
stepfac=1;
while(niter<NSTEP)
{
pcnt++;
if(pcnt==ptot)
{
for(j=i=0;i<ptot;++i)j+=pcheck[i];
stepfac = stepfac*pow(0.9,5-j);
if(!ThisTask)printf("STEPFAC %f %d %d\n",stepfac,j,count);
pcnt=0;
}
/* stepfac=0.7; */
for(i=1;i<=n;++i)
a[i] = (1+gasdev(&IDUM)*start_dev[i]*stepfac)*aprev[i];
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i])OMEGA_M = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])SPECTRAL_INDX = a[++j];
/* if(HOD.free[++i])SIGV = a[++j]; */
}
if(VBIAS_C<0)continue;
/* Hard-wire CVIR variation
*/
if(HOD.free[6])
CVIR_FAC = a[3];
/* Draw random value of cvir from prior.
*/
/* if(CVIR_FAC<0.3 || CVIR_FAC>1.2)continue; */
/* CVIR_FAC = 0.9*drand48()+0.3; */
/* GAMMA = gasdev(&IDUM)*0.02 + 0.15; */
chi2=chi2_wp_wrapper(a);
if(MCMC>2 && chi2<1.0E7)chi2+=chi2_zspace(a);
if(!ThisTask){
printf("TRY %d ",++count);
for(i=1;i<=n;++i)
printf("%.4e ",a[i]);
printf("%e\n",chi2);fflush(stdout);
}
pcheck[pcnt]=1;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
{
/* This for loop puts the prev element in the chain is
* the current trial point is rejected.
*/
/* For the initialization, don't use this: we need
* separate elements for estimating the covariance matrix.
*/
/*
for(i=1;i<=n;++i)
a[i] = aprev[i];
chi2 = chi2prev;
*/
if(USE_IWEIGHT)
iweight[niter+1]++;
pcheck[pcnt]=0;
continue;
}
niter++;
iweight[niter]++;
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask){
printf("ACCEPT %d %d ",niter,count);
for(i=1;i<=n;++i)
printf("%e ",a[i]);
printf("%e\n",chi2);fflush(stdout);
printf("HSTATS %d %e %e %e %e\n",niter,HOD.M_min,number_weighted_halo_mass(),
number_weighted_central_mass(),
qromo(func_satellite_density,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY);
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
printf("FSAT %d %e %e %e %e\n",niter,fsat,HOD.M_min,HOD.sigma_logM,CVIR_FAC);
}
}
RESTART_POINT:
stepfac=1.6/sqrt(n);
pcnt=-1;
t0 = second();
NSTEP = niter;
while(niter<imax_chain)
{
pcnt++;
if(pcnt==ptot)
{
for(j=i=0;i<ptot;++i)j+=pcheck[i];
//stepfac=1.6/sqrt(n);
//stepfac=0.6/sqrt(n);
// stepfac = stepfac*pow(0.9,6-j);
stepfac = 0.25;
stepfac = 0.5;
stepfac=1.6/sqrt(n);
if(!ThisTask)printf("STEPFAC %f %d %d\n",stepfac,j,count);
pcnt=0;
}
//stepfac = 0;
if(convergence)goto SKIP_MATRIX;
// if(niter>NSTEP_MAX && niter%NSTEP_MAX!=0)goto SKIP_MATRIX;
for(j=1;j<=n;++j)
{
avg1[j]=0;
for(k=1;k<=n;++k)
cov1[j][k]=0;
}
total_weight = 0;
for(i=1;i<=niter;++i)
{
for(j=1;j<=n;++j)
{
avg1[j]+=chain[i][j]*iweight[i];
for(k=1;k<=n;++k)
cov1[j][k]+=chain[i][j]*chain[i][k]*iweight[i];
}
total_weight+=iweight[i];
}
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
tmp[i][j] = cov1[i][j]/total_weight - avg1[i]*avg1[j]/(total_weight*total_weight);
jacobi(tmp,n,eval,evect,&nrot);
gaussj(evect,n,tmp1,1);
SKIP_MATRIX:
if(RESTART==4)convergence = 1;
for(i=1;i<=n;++i)
atemp[i] = gasdev(&IDUM)*sqrt(eval[i])*stepfac;
for(i=1;i<=n;++i)
for(a[i]=0,j=1;j<=n;++j)
a[i] += atemp[j]*evect[j][i];
for(i=1;i<=n;++i)
a[i] += aprev[i];
if(!ThisTask)
for(i=1;i<=n;++i)
{
printf("COV %d %d %e ",count,i,sqrt(eval[i]));
for(j=1;j<=n;++j)
printf("%e ",evect[j][i]);
printf("\n");
}
/* Using only variances
*/
//for(i=1;i<=n;++i)
// a[i] = aprev[i] + gasdev(&IDUM)*sqrt(tmp[i][i])*stepfac;
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i])OMEGA_M = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])SPECTRAL_INDX = a[++j];
/* if(HOD.free[++i])SIGV = a[++j]; */
}
if(VBIAS_C<0)continue;
/* Hard-wire CVIR variation
*/
if(HOD.free[6])
CVIR_FAC = a[3];
/* Draw random value of cvir from prior.
*/
/* CVIR_FAC = a[n]; */
/* if(CVIR_FAC<0.3 || CVIR_FAC>1.2)continue; */
/* CVIR_FAC = 0.7*drand48()+0.3; */
/* GAMMA = gasdev(&IDUM)*0.02 + 0.15; */
// printf("GAMMA %d %f %f\n",count+1,GAMMA,CVIR_FAC);
chi2=chi2_wp_wrapper(a);
if(MCMC>2 && chi2<1.0E7)chi2+=chi2_zspace(a);
tprev = t0;
t0 = second();
if(!ThisTask) {
printf("TRY %d ",++count);
for(i=1;i<=n;++i)
printf("%.4e ",a[i]);
printf("%e %.2f\n",chi2,timediff(tprev,t0));fflush(stdout);
}
pcheck[pcnt]=0;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
{
/*
for(i=1;i<=n;++i)
a[i] = aprev[i];
chi2 = chi2prev;
*/
if(USE_IWEIGHT)
iweight[niter+1]++;
continue;
}
pcheck[pcnt]=1;
// if(NSTEP<NSTEP_MAX)NSTEP++;
niter++;
if(!convergence)NSTEP = niter;
iweight[niter]++;
if(niter%NSTEP_MAX==0 && !convergence && niter>NSTEP_MAX)
{
convergence = 1;
for(i=1;i<=n;++i)
{
x1=fabs(eval[i]-eval_prev[i])/eval_prev[i];
if(x1>0.01)convergence = 0;
printf("CONVERGENCE CHECK %d %d %e %e %e\n",niter/NSTEP_MAX,i,x1,eval[i],eval_prev[i]);
}
for(i=1;i<=n;++i)
eval_prev[i] = eval[i];
convergence = 0;
if(convergence)
printf("CONVERGENCE ACCOMPLISHED %d %d \n",niter,count);
}
if(niter==NSTEP_MAX)
{
for(i=1;i<=n;++i)
eval_prev[i] = eval[i];
}
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask) {
printf("ACCEPT %d %d ",niter,count);
for(i=1;i<=n;++i)
printf("%e ",a[i]);
printf("%e\n",chi2);fflush(stdout);
printf("HSTATS %d %e %e %e %e\n",niter,HOD.M_min,number_weighted_halo_mass(),
number_weighted_central_mass(),
qromo(func_satellite_density,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY);
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
printf("FSAT %d %e %e %e %e\n",niter,fsat,HOD.M_min,HOD.sigma_logM,CVIR_FAC);
}
}
}
double chi2_wp_wrapper(double *a)
{
static int flag=1;
static double *b;
int i,j;
if(flag)
{
b=dvector(1,100);
flag=0;
}
for(j=0,i=1;i<=N_HOD_PARAMS;++i) {
if(HOD.free[i] && i!=5) {
if(a[++j]<=0) { printf("NEG %d %d %e\n",i,j,a[j]); return(1.0E7); } }
if(HOD.free[i] && i==5) {
++j; }
}
i=0;j=0;
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_min */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M1 */
if(HOD.free[++i]){j++;b[j]=a[j];} /* alpha */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_cut */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* sigma_logM */
if(HOD.free[++i]){j++;b[j]=a[j];} /* cvir_fac */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* MaxCen */
if(HOD.free[++i]){j++;b[j]=pow(10.0,a[j]);} /* M_sat_break */
if(HOD.free[++i]){j++;b[j]=a[j];} /* alpha1 */
return(chi2_wp(b));
}
double mcmc_initialize(double *a, double **cov1, double *avg1, double *start_dev)
{
int i,j=0;
double x1,x2;
long IDUM = -556;
i=0;j=0;
if(HOD.free[++i]){ a[++j]=log10(HOD.M_min);start_dev[j]=0.001; }
if(HOD.free[++i]){ a[++j]=log10(HOD.M1);start_dev[j]=0.001; } //.0005
if(HOD.free[++i]){ a[++j]=HOD.alpha;start_dev[j]=0.03; } //.005
if(HOD.free[++i]){ a[++j]=log10(HOD.M_cut);start_dev[j]=0.01; } //.001
if(HOD.free[++i]){ a[++j]=log10(HOD.sigma_logM);start_dev[j]=0.01; }
if(HOD.free[++i]){ a[++j]=CVIR_FAC;start_dev[j]=0.02; }
if(HOD.pdfc==7) {
if(HOD.free[++i])a[++j]=log10(HOD.M_cen_max); start_dev[j]=0.001; }
else {
if(HOD.free[++i])a[++j]=HOD.MaxCen; start_dev[j]=0.02; }
if(HOD.free[++i]){ a[++j]=log10(HOD.M_sat_break);start_dev[j]=0.001; }
if(HOD.free[++i]){ a[++j]=HOD.alpha1;start_dev[j]=0.02; }
if(MCMC>1)
{
if(HOD.free[++i])a[++j]=OMEGA_M;
if(HOD.free[++i])a[++j]=SIGMA_8;
if(HOD.free[++i])a[++j]=VBIAS;
if(HOD.free[++i])a[++j]=VBIAS_C;
if(HOD.free[++i])a[++j]=GAMMA;
if(HOD.free[++i])a[++j]=SPECTRAL_INDX;
}
printf("INITIAL VALUES: ");
for(i=1;i<=wp.ncf;++i)printf("%e ",a[i]);
printf("\n");
for(i=1;i<=wp.ncf;++i)
{
avg1[i]=a[i];
for(j=1;j<=wp.ncf;++j)
cov1[i][j]=a[i]*a[j];
}
if(MCMC>1)
{
RESET_COSMOLOGY++;
j=0;
for(i=1;i<=N_HOD_PARAMS;++i)if(HOD.free[i])j++;
i=N_HOD_PARAMS;
if(HOD.free[++i]){ OMEGA_M = a[++j]; start_dev[j] = 0.01; }
if(HOD.free[++i]){ SIGMA_8 = a[++j]; start_dev[j] = 0.01; }
if(HOD.free[++i]){ VBIAS = a[++j]; start_dev[j] = 0.01; }
if(HOD.free[++i]){ VBIAS_C = a[++j]; start_dev[j] = 0.02; }
if(HOD.free[++i]){ GAMMA = a[++j]; start_dev[j] = 0.015; }
if(HOD.free[++i]){ SPECTRAL_INDX = a[++j]; start_dev[j] = 0.02; }
}
x1=chi2_wp_wrapper(a);
if(MCMC>2)
x2=chi2_zspace(a);
else
x2=0;
if(!ThisTask) {
printf("TRY 0 ");
for(i=1;i<=wp.ncf;++i)
printf("%.4e ",a[i]);
printf("%e\n",x1+x2);fflush(stdout);
printf("INITIAL CHI2: %e %e\n",x1,x2);
fflush(stdout);
}
return(x1+x2);
}
/* This is to look at a chain and get the variances in each parameter.
*/
void mcmc_restart2(double *start_dev, int np)
{
int n,i,j,k,i1,i2;
FILE *fp;
char aa[100];
float xbar[10],xsqr[10],x;
fp = openfile(RESTART_FILE);
n = filesize(fp);
for(i=0;i<np;++i)
xbar[i] = xsqr[i] = 0;
for(i=1;i<=n;++i)
{
fscanf(fp,"%s %d %d",aa,&i1,&i2);
for(j=0;j<np;++j)
{
fscanf(fp,"%f",&x);
xbar[j] += x;
xsqr[j] += x*x;
}
fgets(aa,100,fp);
}
for(i=0;i<np;++i)
{
xbar[i]/=n;
xsqr[i]/=n;
xsqr[i] = sqrt(xsqr[i] - xbar[i]*xbar[i]);
start_dev[i+1] = xsqr[i];
if(!ThisTask)
fprintf(stdout,"RESTART_DEV %f %f\n",xbar[i],xsqr[i]);
}
}
int mcmc_restart3(double **chain, int n, double *chi2_prev, int *iweight)
{
FILE *fp;
char aa[100];
int niter,i,j,i1,i2,iprev;
double x,*a,chi2;
fp = openfile(RESTART_FILE);
niter = filesize(fp);
a = dvector(1,n);
fscanf(fp,"%s %d %d",aa,&i1,&i2);
rewind(fp);
iprev = i2 - 1;
for(i=1;i<=niter;++i)
{
fscanf(fp,"%s %d %d",aa,&i1,&i2);
if(USE_IWEIGHT)
iweight[i] = i2 - iprev;
else
iweight[i] = 1;
iprev = i2;
for(j=1;j<=n;++j)
fscanf(fp,"%lf",&chain[i][j]);
fscanf(fp,"%lf",&x);
}
*chi2_prev = 20000*x; // set it to automatically take the first element
//*chi2_prev = x;
/* Normalize all the masses by OMEGA_M
*/
for(i=1;i<=niter;++i)
{
chain[i][1] -= log10(chain[i][4]);
chain[i][3] -= log10(chain[i][4]);
}
return niter;
}

39
c_tools/hod/mstar.c
View file

@ -1,39 +0,0 @@
#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
View file

@ -1,204 +0,0 @@
#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;
}

231
c_tools/hod/one_halo_rspace.c
View file

@ -1,231 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "header.h"
/* These routines control the real-space one-halo term.
* For specifics, see:
*
* Berlind, A.\ A., \& Weinberg, D.\ H.\ 2002, \apj, 575, 587
* Zheng, Z. 2003, \apj, 610, 61
* Tinker, Weinberg, Zheng, Zehavi 2005 Apj 631 (App B)
*
*/
/* Local functions.
*/
void calc_real_space_one_halo(double *r, double *xi, int n);
double func1(double m);
double func1cs(double m);
double func1satsat(double m);
double func1_xcorr(double m);
double one_halo_ss(double r);
double one_halo_cs(double r);
double *xi_cs_g2,*xi_ss_g2,*xi_rad_g2;
/* These are the local globals to use during the qromo integration
*/
double r_g2;
/* This function tabulates the one-halo real-space term for spline interpolation.
* If the requested radius is out-of-bounds of the tabulated function, a value of
* zero is returned.
*/
double one_halo_real_space(double r)
{
static int flag=0;
static double *x,*y,*y2;
int i,n=100;
double a;
if(!HOD.pdfs)return(0);
if(!flag || RESET_FLAG_1H)
{
if(!flag)
{
x=dvector(1,n);
y=dvector(1,n);
y2=dvector(1,n);
}
flag=1;
RESET_FLAG_1H=0;
calc_real_space_one_halo(x,y,n);
spline(x,y,n,2.0E+30,2.0E+30,y2);
}
if(r>x[n])return(0);
if(r<x[1])return(0);
splint(x,y,y2,n,r,&a);
return(a);
}
/* Here we calculate the one-halo real space term
* logarithmically spaced in r. The minimum value of r = 0.01 Mpc/h. The maximum
* value of r is set to be approximately twice the virial radius of M_max.
*
* Only halos with virial radii greater than 1/2 the separation
* contribute to the 1-halo term.
* Terminate integrations when r>2*R_vir(M_max).
*/
void calc_real_space_one_halo(double *r, double *xi, int n)
{
static int ncnt=0;
double fac,s1,rhi=1,rlo=-2,dr,mlo,x1,x2;
int i,j;
FILE *fp;
char fname[100];
ncnt++;
rlo=log(0.01);
rhi=log(1.9*pow(3*HOD.M_max/(4*PI*DELTA_HALO*RHO_CRIT*OMEGA_M),1.0/3.0));
dr=(rhi-rlo)/(n-1);
if(OUTPUT>1)
printf("calc_one_halo> starting...\n");
if(!XCORR)
GALAXY_DENSITY2 = GALAXY_DENSITY;
for(i=1;i<=n;++i)
{
r_g2=r[i]=exp((i-1)*dr + rlo);
fac=1.0/(2*PI*r_g2*r_g2*GALAXY_DENSITY*GALAXY_DENSITY2);
mlo = 4./3.*PI*RHO_CRIT*DELTA_HALO*OMEGA_M*pow(r[i]*.5,3.0);
if(mlo<HOD.M_low)
mlo = HOD.M_low;
if(XCORR)
s1=fac*qromo(func1_xcorr,log(mlo),log(HOD.M_max),midpnt)*0.5;
else
s1=fac*qromo(func1,log(mlo),log(HOD.M_max),midpnt);
xi[i]=s1;
if(OUTPUT>1)
printf("calc_one_halo> %f %e %e\n",r[i],s1,fac);
continue;
}
}
/* This is the function passed to qromo in the above routine.
* It is the number density of
* galaxy pairs in halos of mass m at separation r_g2.
* See Equation (11) from Berlind & Weinberg.
*/
double func1(double m)
{
double N,n,fac2,rvir,f_ss,f_cs,cvir,x,rfof,ncen,nsat;
m=exp(m);
cvir=halo_concentration(m)*CVIR_FAC;
n=dndM_interp(m);
nsat=N_sat(m);
ncen=N_cen(m);
rvir=2*pow(3.0*m/(4*DELTA_HALO*PI*OMEGA_M*RHO_CRIT),1.0/3.0);
/* Break up the contribution of pairs into
* central-satellite (cs) and satellite-satellite (ss) pairs.
*/
f_ss=dFdx_ss(r_g2/rvir,cvir)*moment_ss(m)*0.5;
f_cs=dFdx_cs(r_g2/rvir,cvir)*nsat*ncen;
x=n*(f_ss+f_cs)/rvir*m;
return(x);
}
double func1satsat(double m)
{
double N,n,fac2,rvir,f_ss,f_cs,cvir,x,rfof,ncen,nsat;
m=exp(m);
cvir=halo_concentration(m)*CVIR_FAC;
n=dndM_interp(m);
nsat=N_sat(m);
ncen=N_cen(m);
rvir=2*pow(3.0*m/(4*DELTA_HALO*PI*OMEGA_M*RHO_CRIT),1.0/3.0);
/* Break up the contribution of pairs into
* central-satellite (cs) and satellite-satellite (ss) pairs.
*/
f_ss=dFdx_ss(r_g2/rvir,cvir)*moment_ss(m)*0.5;
x=n*(f_ss)/rvir*m;
return(x);
}
double func1cs(double m)
{
double N,n,fac2,rvir,f_ss,f_cs,cvir,x,rfof,ncen,nsat;
m=exp(m);
cvir=halo_concentration(m)*CVIR_FAC;
n=dndM_interp(m);
nsat=N_sat(m);
ncen=N_cen(m);
rvir=2*pow(3.0*m/(4*DELTA_HALO*PI*OMEGA_M*RHO_CRIT),1.0/3.0);
/* Break up the contribution of pairs into
* central-satellite (cs) and satellite-satellite (ss) pairs.
*/
f_cs=dFdx_cs(r_g2/rvir,cvir)*nsat*ncen;
x=n*(f_cs)/rvir*m;
return(x);
}
/* Same as above, only now we're calculating the number of pairs
* in the cross-correlation.
*
* NB! -- We're assuming that the satellite galaxy profiles
* of both HODs are the same.
*/
double func1_xcorr(double m)
{
double N,n,fac2,rvir,f_ss=0,f_cs=0,cvir,x,rfof,ncen1,nsat1,ncen2,nsat2;
m=exp(m);
cvir=halo_concentration(m)*CVIR_FAC;
n=dndM_interp(m);
nsat1=N_sat(m);
ncen1=N_cen(m);
nsat2=N_sat2(m);
ncen2=N_cen2(m);
rvir=2*pow(3.0*m/(4*DELTA_HALO*PI*OMEGA_M*RHO_CRIT),1.0/3.0);
/* Break up the contribution of pairs into
* central-satellite (cs) and satellite-satellite (ss) pairs.
* But for x-corr, we have c1-s2, c2-s1, s1-s2.
*/
f_ss=dFdx_ss(r_g2/rvir,cvir)*nsat1*nsat2;
f_cs=dFdx_cs(r_g2/rvir,cvir)*(nsat1*ncen2 + nsat2*ncen1);
x=n*(f_ss+f_cs)/rvir*m;
return(x);
}
/* The public version doesn't currently support cross-correlations.
*/
double N_sat2(double m)
{
return 0;
}
double N_cen2(double m)
{
return 0;
}

175
c_tools/hod/pair_density.c
View file

@ -1,175 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "header.h"
double m_g3,r_g3,*x_g3,*y_g3,*z_g3;
double n_g3;
int flag_g3=0;
double func_ng(double m);
double func_ng2(double m);
/* The restr
*/
double restricted_number_density(double r)
{
static int flag=1;
static double *x,*y,*y2;
int i,n=50,j;
double mlimit,dlogm,logm,mmin,sum=0,t0,t1,s1,s2,s3,m,r1,r2,ng2,rlim,rmin,rmax;
if(flag)
{
n_g3 = n;
x_g3=dvector(1,n);
y_g3=dvector(1,n);
z_g3=dvector(1,n);
flag=0;
}
/* Reset the static variables in this function.
*/
func_ng2(-1);
r_g3=r;
ng2=GALAXY_DENSITY*GALAXY_DENSITY;
/* Calculate the maximum allowable halo mass, which had
* rvir = r_g3 - rvir(M_low).
*/
r1=pow(3.*HOD.M_low/(4.*PI*DELTA_HALO*RHO_CRIT*OMEGA_M),1.0/3.0);
rlim = r_g3 - r1;
mlimit=log(4./3.*DELTA_HALO*RHO_CRIT*PI*rlim*rlim*rlim*OMEGA_M);
if(mlimit>log(HOD.M_max))mlimit=log(HOD.M_max);
mmin=log(HOD.M_low);
if(HOD.color==2)
{
dlogm=(mlimit-mmin)/(n-1);
m = mmin;
for(i=1;i<=n;++i)
{
if(N_avg(exp(m))>0)break;
m += dlogm;
}
mmin = m;
r1=pow(3.*exp(mmin)/(4.*PI*DELTA_HALO*RHO_CRIT*OMEGA_M),1.0/3.0);
rlim = r_g3 - r1;
mlimit=log(4./3.*DELTA_HALO*RHO_CRIT*PI*rlim*rlim*rlim*OMEGA_M);
}
if(EXCLUSION==2) {
dlogm=(mlimit-mmin)/(n-1);
x_g3[1] = mmin;
y_g3[1] = qromo(func_galaxy_density,mmin,mmin+dlogm,midpnt);
for(i=2;i<=n;++i)
{
x_g3[i] = i*dlogm+mmin;
y_g3[i] = y_g3[i-1] + qromo(func_galaxy_density,(i-1)*dlogm+mmin,mmin+i*dlogm,midpnt);
}
spline(x_g3,y_g3,n,1.0E+30,1.0E+30,z_g3);
s1 = qromo(func_ng2,mmin,mlimit,midpnt);
return(sqrt(s1));
}
/* Calculate the double integral at specified masses.
*/
dlogm=(mlimit-mmin)/(n-1);
for(i=1;i<=n;++i)
{
logm=(i-0.5)*dlogm+mmin;
m_g3=exp(logm);
r2 = pow(3*m_g3/(4*PI*DELTA_HALO*RHO_CRIT*OMEGA_M),1.0/3.0);
if(EXCLUSION==3) {
if(ellipsoidal_exclusion_probability(r1/r2,r_g3/(r1+r2))==0)break; }
else {
if(r1+r2>r_g3)break; }
s1=qtrap(func_ng,mmin,mlimit,1.0E-4);
sum+=s1*m_g3*dlogm;
if(s1==0)break;
if(sum>=ng2)break;
}
return sqrt(sum);
}
double func_ng2(double m)
{
static double fac2=-1,fac1=-1;
double s1,rv1,n,N,m1,mx;
if(m<0)
{
fac1=fac2=-1;
return(0);
}
m1=exp(m);
if(fac2<0)
fac2=pow(3.0/(4.*PI*DELTA_HALO*RHO_CRIT*OMEGA_M),1.0/3.0);
if(fac1<0)
fac1=4./3.*PI*RHO_CRIT*DELTA_HALO*OMEGA_M;
rv1 = r_g3 - pow(m1,1.0/3.0)*fac2;
rv1 = rv1;
mx = fac1*rv1*rv1*rv1;
n=dndM_interp(m1);
N=N_avg(m1);
splint(x_g3,y_g3,z_g3,n_g3,log(mx),&s1);
return(n*N*s1*m1);
}
double func_ng(double m)
{
static double fac2=-1;
double s1,rv1,rv2,exfac=1,n,N;
m=exp(m);
if(fac2<0)
fac2=pow(3.0/(4*DELTA_HALO*PI*RHO_CRIT*OMEGA_M),1.0/3.0);
rv1=pow(m_g3,1.0/3.0)*fac2;
rv2=pow(m,1.0/3.0)*fac2;
if(EXCLUSION==3)
{
if(0.5*(rv1+rv2)>r_g3)return(0);
if(1.5*(rv1+rv2)>r_g3)exfac=ellipsoidal_exclusion_probability(rv2/rv1,r_g3/(rv2+rv1));
}
else
{
if(rv1+rv2>r_g3)return(0);
}
n=dndM_interp(m)*dndM_interp(m_g3);
N=N_avg(m)*N_avg(m_g3);
return(exfac*n*N*m);
}
/* This is the probability that two halos do not overlap, given their
* radii and separation. Of course, for spherical halos P(x) is a step function
* at x = (r1+r2)/r_sep = 1, but for ellipsoidal halos there is a chance
* that they could be closer. In detail, P(x) changes depending on the mass
* ratio of the halos, but using tabulated values does not appear to make
* significant difference in the results for xi_2h(r). The function below is
* a fit to Monte Carlo results for a halos with a distribution of axis ratios
* which is lognormal in e_b = (1-b/a) and e_c = (1-c/a) with dispersions of 0.2
* mean <b/a>=0.9 and <c/a>=0.8 (pretty reasonable values).
*/
double ellipsoidal_exclusion_probability(double rv, double r)
{
static int flag=0,nr=101,nratio=31;
static double **xprob,*rad,*ratio,rhi,rlo,mhi,mlo,dr,dm;
float x1,x2,x3;
int i,j,im,ir;
FILE *fp;
if(rv<1)rv=1.0/rv;
r=(r-0.8)/0.29;
if(r>1)return(1.0);
if(r<0)return(0.0);
return(3*r*r-2*r*r*r);
}

5
c_tools/hod/test.c
View file

@ -1,5 +0,0 @@
#include "header.h"
void test(int argc, char **argv)
{
}

140
c_tools/hod/tf_eisenstein_hu.c
View file

@ -1,140 +0,0 @@
#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
View file

@ -1,70 +0,0 @@
/* 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
View file

@ -1,59 +0,0 @@
#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));
}

511
c_tools/hod/two_halo_rspace.c
View file

@ -1,511 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "header.h"
/* This is the two-halo term of the correlation function in both real space.
* This is solved in Fourier space and requires the
* non-linear power spectrum P_nl(k), which is from the model of Smith et al.
* (see Smith et al, astro-ph/0207664).
*
* For specifics see Berlind & Weinberg (2002), Zheng (2003), Tinker et al (2005)
*
*/
/* Internal functions.
*/
double func2(double m);
double func2_cen(double m);
double func2_sat(double m);
double func5(double xk);
double func_mlimit(double m);
void calc_real_space_two_halo(double *r, double *xi, int *n);
double HOD2_two_halo_real_space(double r);
/* the restricted number density.
*/
double NG_MINUS2;
/* Globals needed for the qromo functions.
*/
double r_g1,
k1;
/* Global checkflag to stop doing the restricted
* number density.
*/
int checkflag;
/* This tabulates the two-halo real-space term for future
* spline interpolation.
*/
double two_halo_real_space(double r)
{
static int flag=0,n=45;
static double *x,*y,*y2;
int i;
double max=16,min=9,a,rvir_min;
float x1,x2;
FILE *fp;
// if(r<25)
// return(nbody_two_halo(r));
if(!flag || RESET_FLAG_2H)
{
if(!LINEAR_PSP)
nonlinear_sigmac(8.0);
n=30;
if(!flag)
{
x=dvector(1,n);
y=dvector(1,n);
y2=dvector(1,n);
}
RESET_FLAG_2H=0;
flag=1;
if(OUTPUT>1)
fprintf(stdout,"Calculating real-space two-halo term...\n");
calc_real_space_two_halo(x,y,&n);
if(!HOD.color==2)
check_for_smoothness(x,y,n,1.5);
XI_MAX_RADIUS=x[n];
if(XCORR)
for(i=1;i<=n;++i)
{
a = HOD2_two_halo_real_space(x[i]);
y[i] = a*y[i];
if(y[i]<0)y[i]=-1;
else y[i] = sqrt(y[i]);
}
spline(x,y,n,2.0E+30,2.0E+30,y2);
}
if(r<R_MIN_2HALO)return(-1);
if(r>XI_MAX_RADIUS)return(-1);
splint(x,y,y2,n,r,&a);
/* This check is against the spline interpolation, which
* sometimes has (smoothness) problems with the sharp truncation of xi_2h
* at small separations.
*/
if(a<-1)return(-1);
return(a);
}
void calc_real_space_two_halo(double *r, double *xi, int *nn)
{
double xtemp[200],rtemp[200];
double mlimit,s2,rlo,rhi=90,dr,klo,tolerance=1.0e-7,s1;
int i,j,imin=0,n;
double t0,t1,t2,t1s=0,t2s=0;
n=*nn;
/* Set the minimum separation of two-halo pairs.
*/
if(HOD.color)
{
while(N_avg(HOD.M_low)<0.001)
HOD.M_low*=1.01;
}
rlo = 2.2*pow(3*HOD.M_low/(4*DELTA_HALO*PI*OMEGA_M*RHO_CRIT),1.0/3.0);
if(EXCLUSION==3)
rlo = rlo/1.3;
if(EXCLUSION==4)
rlo = rlo/2.1;
R_MIN_2HALO = rlo;
rlo = log(rlo);
dr=(log(rhi)-rlo)/(n-1);
checkflag=0;
for(i=1;i<=n;++i)
{
r_g1=r[i]=exp(rlo+dr*(i-1));
if(r_g1>30)
{
GALAXY_BIAS = qromo(func_galaxy_bias,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
xi[i]=xi_interp(r_g1)*GALAXY_BIAS*GALAXY_BIAS;
continue;
}
/* Fourier transform. If BOX_SIZE specified, then truncate
* transform at that k. If no box, take integral to k=0.
*/
klo = 0;
if(BOX_SIZE)klo=1/BOX_SIZE;
j=16;
s1 = qromo(func5,klo,j*TWOPI/r_g1,midpnt);
s2 = s1;
klo = j*TWOPI/r_g1;
while(fabs(s1)>tolerance*s2) {
j+=16;
s1 = qromo(func5,klo,j*TWOPI/r_g1,midpnt);
s2 += s1;
klo = j*TWOPI/r_g1;
}
/* Divide the integral by the restricted number density.
* (This is calculated in the pk_gg_2h function.)
*/
s2=s2/NG_MINUS2;
/* Correction factor b/c we haven't been
* including all halos.
*/
xi[i]=(NG_MINUS2/GALAXY_DENSITY/GALAXY_DENSITY)*(1+s2)-1;
if(isnan(xi[i]))xi[i]=-1;
if(xi[i]==-1)imin=i;
if(OUTPUT>1)
printf("calc_2halo> %f %f %d\n",r[i],xi[i],checkflag);
//printf("calc_2halo> %f %f %d\n",r[i],xi[i],checkflag);
}
/* Eliminate all entries which have xi=-1
* (to eliminate a "wavy" spline fit at small r
*/
for(j=0,i=imin+1;i<=n;++i)
{
j++;
xtemp[j]=xi[i];
rtemp[j]=r[i];
}
n=j;
for(i=1;i<=n;++i)
{
r[i]=rtemp[i];
xi[i]=xtemp[i];
}
*nn=n;
R_MIN_2HALO=r[1];
/*printf("calc_2halo> %d %d %f\n",imin,n,r[1]);*/
}
/* This calculates and tabulates the galaxy power spectrum. This is done
* by taking the galaxy-number-weighted halo bias (with scale-dependent
* halo bias and extended structure of halos taken into account) and multiplying
* it by the non-linear matter power spectrum.
*
* The galaxy-number-weighted average is an integral over the halo mass function
* (see Equations B10 & B12 in Appendix B of Tinker et al.) over all allowed halo
* pairs.
*
* The allowed number of halo pairs is controlled by the type of halo exclusion used:
* EXCLUSION = 1 --> only halos with Rvir < r/2
* EXCLUSION = 2 --> only halo PAIRS with R1+R2 < r ("spherical eclusion")
* EXCLUSION = 3 --> same as 3
* EXCLUSION = 4 --> only halos with Rvir < r
*
* 1 and 4 are fast b/c the integral is separable, calculated once and squared. Can drastically
* underestimate the number of small-sep pairs.
*
* for 2/3, I've incorporated the approximation in Tinker etal 2005 called the n_g-matched
* method, where instead of actually doing the double integral over halo_mass_1 and halo_mass_2,
* i calculate the number density of that type of exclusion, and calculate the effective
* radius for the 1/4-type exlcusion-- ie, halos with Rvir < [x]*r where 0.5<[x]<1.0.
*/
double psp_gg_2h(double k, double r)
{
static double rp=-1,*x,*y,*y2;
static int flag=0;
int n=60,i;
double a,dk,klo=-3,khi=3,xk,s1,s2,mhi,mlo,mhi1,mlo1;
double t1,t0,ttot1,ttot2;
/* The fourier transform is done at each r, so tabulate the power spectrum
* at all k for a given r. If r has changed, re-tabulate.
*/
if(rp!=r)
{
mlo=HOD.M_low;
if(!flag)
{
flag=1;
x=dvector(1,n);
y=dvector(1,n);
y2=dvector(1,n);
}
dk=(khi-klo)/n;
switch(EXCLUSION) {
case 1:
mhi=4./3.*DELTA_HALO*RHO_CRIT*PI*r*r*r*OMEGA_M*0.125;
if(mhi>HOD.M_max)mhi=HOD.M_max;
NG_MINUS2 = qromo(func_galaxy_density,log(mlo),log(mhi),midpnt);
NG_MINUS2*=NG_MINUS2;
for(i=n;i>=1;--i)
{
xk=pow(10.0,klo+dk*i);
k1=xk;
if(nfw_transform(xk,mhi)<0.999)
a=qromo(func2,log(mlo),log(mhi),midpnt);
x[i]=log(xk);
if(LINEAR_PSP)
y[i]=log(linear_power_spectrum(xk)*a*a);
else
y[i]=log(nonlinear_power_spectrum(xk)*a*a);
}
break;
case 2:
case 3:
if(!checkflag)
s2=restricted_number_density(r);
else
s2=GALAXY_DENSITY;
mhi=4./3.*DELTA_HALO*RHO_CRIT*PI*r*r*r*OMEGA_M*0.125;
if(s2>=GALAXY_DENSITY || mhi>HOD.M_max)
{
s2=GALAXY_DENSITY;
mhi=log(HOD.M_max);
checkflag=1;
}
else
{
mhi=0;
NG_MINUS2=s2;
if(N_avg(HOD.M_low*1.001)>0)
{
if(qromo(func_galaxy_density,log(HOD.M_low),log(HOD.M_low*1.0001),midpnt)
>NG_MINUS2)mhi=log(HOD.M_low*1.0001);
}
if(qromo(func_galaxy_density,log(HOD.M_low),log(HOD.M_max*4.0),midpnt)
<NG_MINUS2)mhi=log(HOD.M_max);
if(!mhi)
mhi=zbrent(func_mlimit,log(HOD.M_low*1.0001),log(HOD.M_max*4.0),1.0E-4);
if(mhi>log(HOD.M_max))mhi=log(HOD.M_max);
}
NG_MINUS2=s2*s2;
for(i=n;i>=1;--i)
{
xk=pow(10.0,klo+dk*i);
k1=xk;
if(nfw_transform(xk,exp(mhi))<0.999) {
if(HOD.pdfc==7)/* || HOD.pdfc==6 || HOD.pdfc==9 || HOD.pdfc==8)*/
{
a = qromo(func2_cen,log(HOD.M_low),log(HOD.M_cen_max),midpnt);
a += qromo(func2_sat,log(HOD.M_low),log(HOD.M_max),midpnt);
}
else {
MAG_21_CHECK:
a=qromo(func2,log(mlo),mhi,midpnt);
}
}
x[i]=log(xk);
if(LINEAR_PSP)
y[i]=log(linear_power_spectrum(xk)*a*a);
else
y[i]=log(nonlinear_power_spectrum(xk)*a*a);
}
break;
case 4: // where max halo mass is M(R=r) rather than M(R=r/2)
mhi=4./3.*DELTA_HALO*RHO_CRIT*PI*r*r*r*OMEGA_M;
if(mhi>HOD.M_max)mhi=HOD.M_max;
NG_MINUS2 = qromo(func_galaxy_density,log(mlo),log(mhi),midpnt);
NG_MINUS2*=NG_MINUS2;
for(i=n;i>=1;--i)
{
xk=pow(10.0,klo+dk*i);
k1=xk;
if(nfw_transform(xk,mhi)<0.999)
a=qromo(func2,log(mlo),log(mhi),midpnt);
x[i]=log(xk);
if(LINEAR_PSP)
y[i]=log(linear_power_spectrum(xk)*a*a);
else
y[i]=log(nonlinear_power_spectrum(xk)*a*a);
}
break;
default:
endrun("Error: invalid choice of EXCLUSION");
}
spline(x,y,n,2.0E+30,2.0E+30,y2);
rp=r;
}
k=log(k);
splint(x,y,y2,n,k,&a);
return(exp(a));
}
/* The routine called by qromo for the integral in Zheng's Eq [7]
*/
double func2(double m)
{
double n,N,b,yg,x,Ncen,Nsat;
m=exp(m);
n=dndM_interp(m);
b=bias_interp(m,r_g1);
Ncen=N_cen(m);
Nsat=N_sat(m);
yg=nfw_transform(k1,m);
x = n*(Ncen + Nsat*yg)*b*m;
return(x);
}
/* The routine called by qromo for the integral in Zheng's Eq [7]
* FOR CENTRAL GALAXIES ONLY: This means no Fourier transform of NFW profile.
*/
double func2_cen(double m)
{
double n,N,b,yg,x;
m=exp(m);
n=dndM_interp(m);
b=bias_interp(m,r_g1);
N=N_cen(m);
x=n*N*b*m;
return(x);
}
/* The routine called by qromo for the integral in Zheng's Eq [7]
* FOR SATELLITE GALAXIES ONLY.
*/
double func2_sat(double m)
{
double n,N,b,yg,x;
m=exp(m);
n=dndM_interp(m);
b=bias_interp(m,r_g1);
N=N_sat(m);
yg=nfw_transform(k1,m);
x=n*N*b*yg*m;
return(x);
}
/* This is the integrand of the Fourier transform
* of the two-halo power spectrum to correlation function
*/
double func5(double xk)
{
double psp1;
double xk1,xk2,x3;
if(xk==0)return(0);
psp1=psp_gg_2h(xk,r_g1);
xk1=r_g1*xk;
psp1*=sin(xk1)/xk1/xk;
return(psp1);
}
/* This is the function sent to zbrent to calculate the halo mass
* which gives the matched restricted number density.
*/
double func_mlimit(double m)
{
double s1;
if(N_avg(exp(m))>0)
s1=qromo(func_galaxy_density,log(HOD.M_low),m,midpnt);
else
s1=0;
return(s1-NG_MINUS2);
}
/* This tabulates the two-halo real-space term
* for the second HOD function for future
* spline interpolation.
*/
double HOD2_two_halo_real_space(double r)
{
static int flag=0,n=45;
static double *x,*y,*y2;
int i;
double max=16,min=9,a,rvir_min,galtemp;
float x1,x2;
FILE *fp;
if(!flag || RESET_FLAG_2H)
{
/* Switch HODs temporarily
*/
HODt.M_min = HOD.M_min;
HODt.M_low = HOD.M_low;
HODt.M1 = HOD.M1;
HODt.alpha = HOD.alpha;
HODt.M_cen_max = HOD.M_cen_max;
HODt.sigma_logM = HOD.sigma_logM;
HODt.M_max = HOD.M_max;
HODt.pdfc = HOD.pdfc;
HODt.pdfs = HOD.pdfs;
galtemp = GALAXY_DENSITY;
HOD.M_min = HOD2.M_min;
HOD.M_low = HOD2.M_low;
HOD.M1 = HOD2.M1;
HOD.alpha = HOD2.alpha;
HOD.M_cen_max = HOD2.M_cen_max;
HOD.sigma_logM = HOD2.sigma_logM;
HOD.M_max = HOD2.M_max;
HOD.pdfc = HOD2.pdfc;
HOD.pdfs = HOD2.pdfs;
GALAXY_DENSITY = GALAXY_DENSITY2;
set_HOD_params();
if(!LINEAR_PSP)
nonlinear_sigmac(8.0);
n=30;
if(!flag)
{
x=dvector(1,n);
y=dvector(1,n);
y2=dvector(1,n);
}
RESET_FLAG_2H=0;
flag=1;
if(OUTPUT)
fprintf(stdout,"Calculating HOD2 real-space two-halo term...\n");
calc_real_space_two_halo(x,y,&n);
check_for_smoothness(x,y,n,1.5);
XI_MAX_RADIUS=x[n];
spline(x,y,n,2.0E+30,2.0E+30,y2);
/* Switch HODs back
*/
HOD.M_min = HODt.M_min;
HOD.M_low = HODt.M_low;
HOD.M1 = HODt.M1;
HOD.alpha = HODt.alpha;
HOD.M_cen_max = HODt.M_cen_max;
HOD.sigma_logM = HODt.sigma_logM;
HOD.M_max = HODt.M_max;
HOD.pdfc = HODt.pdfc;
HOD.pdfs = HODt.pdfs;
GALAXY_DENSITY = galtemp;
set_HOD_params();
}
if(r>XI_MAX_RADIUS)return(-1);
if(r<R_MIN_2HALO)return(-1);
splint(x,y,y2,n,r,&a);
/* This check is against the spline interpolation, which
* sometimes has (smoothness) problems with the sharp truncation of xi_2h
* at small separations.
*/
if(a<-1)return(-1);
return(a);
}

745
c_tools/hod/wp_minimization.c
View file

@ -1,745 +0,0 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#ifdef PARALLEL
#include <mpi.h>
#endif
#include "header.h"
/* This is a series of routines to calculate the projected correlation
* function from the real-space one (HOD calculation) and then chi2
* minimize the parameters based on the SDSS data.
*/
int FIT_WITH_COSMO = 0,
FIRST_CALL = 0;
double rp_g1;
double func_wp(double z);
double func_wp_rspace(double z);
double func_wp_matter(double z);
double chi2_wp(double *a);
void wp_input(void);
void initial_wp_values(double *a, double **pp, double *yy);
/* This integrates the calculated real-space correlation function
* along the line-of-sight to get the projected correlation function (wp(rp).
* The value passed if the projected separation.
*
* In most SDSS work, the maximum line of sight separation considered when
* claculating wp is pi=40 Mpc/h, which is the default, but can be set
* in the parameter file if desired.
*
* NB! - note that this routine uses a correction for redshift space distortions.
*/
double projected_xi(double r)
{
double x,zmax;
rp_g1=r*r;
two_halo_real_space(1.0);
zmax=sqrt(XI_MAX_RADIUS*XI_MAX_RADIUS-rp_g1);
if(zmax>wp.pi_max)zmax=wp.pi_max;
zmax = wp.pi_max;
x=2*qromo(func_wp,log(0.001),log(zmax),midpnt);
return(x);
}
/* Same as above, but it projects the real-space correlation function
* directly, with no correction for redshift-space effects.
*/
double projected_xi_rspace(double r)
{
double x,zmax;
rp_g1=r*r;
two_halo_real_space(1.0);
zmax=sqrt(XI_MAX_RADIUS*XI_MAX_RADIUS-rp_g1);
if(zmax>wp.pi_max)zmax=wp.pi_max;
zmax = wp.pi_max;
x=2*qromo(func_wp_rspace,log(0.001),log(zmax),midpnt);
return(x);
}
/* Same as above but for matter.
* I've set the maximum line-of-sight separation to be 50 Mpc/h,
* which should be good for visual comparison purposes.
*/
double projected_xi_matter(double r)
{
double x,zmax;
rp_g1=r*r;
zmax=50.0;
x=2*qtrap(func_wp_matter,log(0.001),log(zmax),1.0E-3);
return(x);
}
/* Function called from qromo/qtrap to get xi->wp
* Note that the two-halo term has the linear Kaiser distortion correction.
*/
double func_wp(double z)
{
double r;
z=exp(z);
r=sqrt(rp_g1 + z*z);
return(z*(one_halo_real_space(r)+linear_kaiser_distortion(r,z)));
}
/* Function called from qromo/qtrap to get xi->wpm but without
* correction for redshift-space distortions in the two-halo term
*/
double func_wp_rspace(double z)
{
double r;
z=exp(z);
r=sqrt(rp_g1 + z*z);
return(z*(one_halo_real_space(r)+two_halo_real_space(r)));
}
/* Function called from qromo/qtrap to get xi->wp (dark matter)
*/
double func_wp_matter(double z)
{
double r;
z=exp(z);
r=sqrt(rp_g1 + z*z);
return(z*(xi_interp(r)));
}
/********************************************************************
* Below is the actual minimization of the wp data.
* Relevant variables: [default]
*
* COVAR -> [1]=use covariance matrixl; 0=diagonal error bars
* DEPROJECTED -> [1]=we're fitting a real-space xi(r) [0]= fitting wp(rp)
*
* HOD.free[] is a vector which holds 1/0 as to whether or not a parameter is going
* to be help constant during the chi^2 minimization. 1==vary 0==constant. The default
* on any will be [0].
*
* i variable
* --- --------
* [1] -> M_min
* [2] -> M1
* [3] -> alpha
* [4] -> M_cut
* [5] -> sigmaM
* [6] -> CVIR_FAC
* [7] -> MaxCen (M_cen_max)
*
* Currently I have no checks on whether the values of this vector line
* up correctly with the specific HOD pdfs, so double-check hod.bat files.
*
* Once the code is finished, it will output the values of the HOD parameters
* to a file called [filename].fit (+ the bias, satellite fraction, and chi^2).
* Then it outputs the mean <N>_M to a file called [filename].HOD.
* Then it will go through all the TASKS asked for in the hod.bat file.
*/
void wp_minimization(char *fname)
{
int n,niter,i,j;
double *a,**pp,*yy,FTOL=1.0E-3,chi2min,s1,dlogm,m;
FILE *fp;
char aa[1000];
fprintf(stderr,"\n\nCHI2 MINIMIZATION OF W_P(R_P) DATA..........\n");
fprintf(stderr, "--------------------------------------------\n\n");
OUTPUT = 0;
FIRST_CALL = 1;
if(POWELL)
FTOL=1.0E-3;
else
FTOL=1.0E-4;
for(n=0,i=1;i<=N_HOD_PARAMS;++i)
{
n+=HOD.free[i];
if(!OUTPUT)continue;
printf("wp_min> free[%i] = %d\n",i,HOD.free[i]);
}
if(XCORR)n*=2;
if(OUTPUT)printf("wp_min> Number of free parameters: %d\n",n);
printf("FNAME %s\n",Task.root_filename);
wp_input();
wp.ncf=n;
a=dvector(1,n);
if(POWELL)
pp=dmatrix(1,n,1,n);
else
pp=dmatrix(1,n+1,1,n);
yy=dvector(1,n+1);
initial_wp_values(a,pp,yy);
printf("IVALS %e %e %e %e\n",a[1],a[2],a[3],a[4]);
if(POWELL)
{
if(OUTPUT)printf("wp_min> starting powell.\n");
powell(a,pp,n,FTOL,&niter,&chi2min,chi2_wp);
chi2min = chi2_wp(a);
}
else
{
if(OUTPUT)printf("wp_min> starting amoeba.\n");
amoeba(pp,yy,n,FTOL,chi2_wp,&niter);
for(i=1;i<=n;++i)a[i]=pp[1][i];
chi2min = chi2_wp(a);
}
s1=qromo(func_galaxy_bias,log(HOD.M_low),log(HOD.M_max),midpnt);
GALAXY_BIAS=s1/GALAXY_DENSITY;
printf("POWELL %e %e ",chi2min,HOD.M_min);
for(i=1;i<=n;++i)printf("%e ",a[i]);
printf(" %f\n",GALAXY_BIAS);
/* These outputs are for easy cut & paste into
* another batch file.
*/
//output_parameter_file(fname);
/* Output the fit and the HOD curve.
*/
printf("FNAME2 %s\n",Task.root_filename);
sprintf(aa,"%s.fit",Task.root_filename);
fp=fopen(aa,"w");
fprintf(fp,"%e %e ",chi2min,HOD.M_min);
for(i=1;i<=n;++i)fprintf(fp,"%e ",a[i]);
fprintf(fp," %f\n",GALAXY_BIAS);
fclose(fp);
sprintf(aa,"%s.HOD",Task.root_filename);
fp=fopen(aa,"w");
dlogm=(log(HOD.M_max)-log(HOD.M_low))/99;
for(i=1;i<=100;++i)
{
m=exp((i-1)*dlogm)*HOD.M_low;
fprintf(fp,"%e %e %e %e\n",m,N_cen(m),N_sat(m),N_avg(m));
}
fclose(fp);
fprintf(stderr,"here\n");
free_dvector(a,1,n);
if(POWELL)
free_dmatrix(pp,1,n,1,n);
else
free_dmatrix(pp,1,n+1,1,n);
free_dvector(yy,1,n+1);
fprintf(stderr,"here\n");
free_dvector(wp.r,1,wp.np);
free_dvector(wp.x,1,wp.np);
free_dvector(wp.e,1,wp.np);
if(COVAR)
free_dmatrix(wp.covar,1,wp.np,1,wp.np);
fprintf(stderr,"done in wp_min\n");
}
double integrated_wp_bin(double r)
{
return(projected_xi(r)*r);
}
double integrated_xi_bin(double r)
{
//return((1+one_halo_real_space(r))*r*r);
return((one_halo_real_space(r)+two_halo_real_space(r))*r*r);
}
double chi2_wp(double *a)
{
static int flag=1,niter=0,ichi=-1;
static double *x,mmin_prev=0,t0=-1,t1,sig_prev=0,chi2_prev,chi2_array[10];
double **tmp,**tmp2,chi2,x1,ta1,ta2,dt1h,dt2h,par_chi,chi2ngal;
int i,j,k,ncf_hod;
double rlo,rhi,rmin,rmax,dlogr,integrated_wp_bin();
if(FIRST_CALL)
{
flag = 1;
FIRST_CALL = 0;
}
t0 = clock();
if(HOD.free[1])FIX_PARAM = 0;
wp.iter=niter;
for(j=0,i=1;i<=N_HOD_PARAMS;++i)
if(HOD.free[i])
if(a[++j]<=0) { printf("%d %e\n",j,a[j]); return(1.0E7); }
ncf_hod = j;
RESET_FLAG_1H=1;
RESET_FLAG_2H=1;
RESET_KAISER++;
i=0;j=0;
if(HOD.free[++i])HOD.M_min=a[++j];
if(HOD.free[++i])HOD.M1=a[++j];
if(HOD.free[++i])HOD.alpha=a[++j];
if(HOD.free[++i])HOD.M_cut=a[++j];
if(HOD.free[++i])HOD.sigma_logM=a[++j];
if(HOD.pdfc!=9) {
if(HOD.free[++i])CVIR_FAC=a[++j];
if(HOD.pdfc>=7) {
if(HOD.free[++i])HOD.M_cen_max=a[++j]; }
else {
if(HOD.free[++i])HOD.MaxCen=a[++j]; }
}
if(HOD.free[++i])HOD.M_sat_break=a[++j];
if(HOD.free[++i])HOD.alpha1=a[++j];
if(XCORR) {
i=0;
if(HOD.free[++i])HOD2.M_min=a[++j];
if(HOD.free[++i])HOD2.M1=a[++j];
if(HOD.free[++i])HOD2.alpha=a[++j];
if(HOD.free[++i])HOD2.M_cut=a[++j];
if(HOD.free[++i])HOD2.sigma_logM=a[++j];
if(HOD2.pdfc!=9) {
if(HOD.free[++i])CVIR_FAC=a[++j];
if(HOD2.pdfc>=7) {
if(HOD2.free[++i])HOD2.M_cen_max=a[++j]; }
else {
if(HOD2.free[++i])HOD2.MaxCen=a[++j]; }
}
}
if(!ThisTask) {
printf("START %d ",niter);
for(i=1;i<=ncf_hod;++i)printf("%e ",a[i]);
printf("\n");
}
if(HOD.pdfs==2 && HOD.M_cut<1.0e7)return(1.0e7);
if(HOD.pdfs==2 && HOD.M_cut>1.0e15)return(1.0e7);
/* if(HOD.M_min>HOD.M_max)return(1.0e7); */
/* I've noticed some problems when sigma_logM gets to be
* unresonably high or low, so I've put some limits on the
* values they can have when doing mag-bin fitting.
*/
if(HOD.pdfc==6) {
if(HOD.sigma_logM<0.07)return(1.0e7);
if(HOD.sigma_logM>1.2)return(1.0e7);
}
if(HOD.pdfc==2 || HOD.pdfc==9) {
if(HOD.sigma_logM>1.8)return(1.0e7);
if(HOD.sigma_logM<0.05)return(1.0e7);
}
if(HOD.M1>1.0e17)return(1.0e7);
if(FIX_PARAM==2)
{
HOD.M1=HOD.M_low;
x1=qromo(func_galaxy_density,log(HOD.M_low),log(HOD.M_max),midpnt);
if(x1<GALAXY_DENSITY)return(1.0e7);
HOD.M1=pow(10.0,14.8);
x1=qromo(func_galaxy_density,log(HOD.M_low),log(HOD.M_max),midpnt);
if(x1>GALAXY_DENSITY)return(1.0e7);
HOD.M1=0;
}
/* Check the make sure these are reasonable parameters
* (Assuming that M_min is NOT a FREE parameter but is
* calculated from the GALAXY_DENSITY.)
*/
if(FIX_PARAM==1 && !HOD.color)
{
HOD.M_min=pow(10.0,8.0);
HOD.M_low=set_low_mass();
if(HOD.M_low<1.0e8)HOD.M_low=1.0e8;
x1=qromo(func_galaxy_density,log(HOD.M_low),log(HOD.M_max),midpnt);
// fprintf(stderr,"PC1 %e %e\n",x1,GALAXY_DENSITY);
if(x1<GALAXY_DENSITY) {
fprintf(stdout,"PCHECK %e %e %e\n",x1,GALAXY_DENSITY,HOD.M_low);
return(1.0e7); }
HOD.M_min=pow(10.0,14.8);
if(HOD.pdfc==7 && HOD.pdfc==8)
HOD.M_min=HOD.M_cen_max*0.99;
HOD.M_low=set_low_mass();
x1=qromo(func_galaxy_density,log(HOD.M_low),log(HOD.M_max),midpnt);
/* fprintf(stderr,"PC2 %e %e %e %e\n",HOD.M_min,HOD.M_low,x1,GALAXY_DENSITY); */
if(x1>GALAXY_DENSITY) {
fprintf(stdout,"PCHECK %e %e\n",x1,GALAXY_DENSITY);
return(1.0e7); }
HOD.M_min=0;
}
if(ERROR_FLAG)
{
ERROR_FLAG=0;
return(1e7);
}
if(HOD.free[0] || HOD.free[1])
GALAXY_DENSITY=0;
if(!HOD.color)
set_HOD_params();
/*
if(XCORR)
set_HOD2_params();
*/
if(HOD.free[0])
{
chi2ngal = (GALAXY_DENSITY-wp.ngal)*(GALAXY_DENSITY-wp.ngal)/wp.ngal_err/wp.ngal_err;
if(chi2ngal>1.0E3)return(chi2ngal);
}
if(ERROR_FLAG)
{
ERROR_FLAG=0;
return(1e7);
}
/* if(HOD.pdfs==3 && HOD.M_cut<HOD.M_low)return(1.1e7); */
if(HOD.M_min>HOD.M1)return(1.0e7);
/* if(HOD.pdfs==3)if(HOD.M_cut<HOD.M_min*0.9)return(1.2e7); */
mmin_prev=HOD.M_min;
sig_prev=HOD.sigma_logM;
if(HOD.color>=niter)
flag = 1;
if(flag && COVAR)
{
// printf("INVERTING COVARIANCE MATRIX\n");
flag=0;
tmp=dmatrix(1,wp.np,1,1);
tmp2=dmatrix(1,wp.np,1,wp.np);
for(i=1;i<=wp.np;++i)
for(j=1;j<=wp.np;++j)
tmp2[i][j]=wp.covar[i][j];
gaussj(tmp2,wp.np,tmp,1);
for(i=1;i<=wp.np;++i)
for(j=1;j<=wp.np;++j)
wp.covar[i][j]=tmp2[i][j];
free_dmatrix(tmp,1,wp.np,1,1);
free_dmatrix(tmp2,1,wp.np,1,wp.np);
x=dvector(1,wp.np);
}
if(!COVAR)
x=dvector(1,wp.np);
rmax = wp.r[wp.np];
rmin = wp.r[1];
dlogr = (log(rmax)-log(rmin))/(wp.np-1);
BETA = pow(OMEGA_M,0.6)/qromo(func_galaxy_bias,log(HOD.M_low),log(HOD.M_max),midpnt)*
GALAXY_DENSITY;
if(OUTPUT)
printf("BETA = %f\n",BETA);
rlo = exp(log(rmin) - 0.5*dlogr);
for(i=1;i<=wp.np;++i)
{
rhi = exp(dlogr)*rlo;
if(DEPROJECTED)
x[i]=one_halo_real_space(wp.r[i])+two_halo_real_space(wp.r[i]);
else
{
x[i]=projected_xi(wp.r[i]);
if(wp.format==3)
x[i]/=wp.r[i];
}
if(OUTPUT && !ThisTask)
printf("WP%d %f %e %e %e %e\n",niter+1,wp.r[i],wp.x[i],x[i],rlo,rhi);
rlo=rhi;
}
if(ERROR_FLAG)
{
ERROR_FLAG=0;
return(1e7);
}
chi2=0;
if(COVAR)
{
for(i=1;i<=wp.np;++i)
for(j=1;j<=wp.np;++j)
{
x1=(x[i]-wp.x[i])*(x[j]-wp.x[j])*wp.covar[i][j];
chi2+=x1;
}
}
else
{
if(!PCA)
{
for(i=1;i<=wp.np;++i)
{
x1=(x[i]-wp.x[i])*(x[i]-wp.x[i])/
(wp.e[i]*wp.e[i] + wp.esys*wp.esys*x[i]*x[i]);
chi2+=x1;
}
}
else
{
chi2=0;
for(i=1;i<=wp.npca;++i)
{
par_chi=0;
for(j=1;j<=wp.np;++j)
par_chi+=wp.covar[j][i]*(x[j]-wp.x[j])/wp.e[j];
chi2+=par_chi*par_chi/wp.eigen[i];
}
}
}
/*
* From Peder Norberg's instructions for use of PCA:
do i=1,npca
par_chi=0.d0
do j=1,npoints
par_chi=par_chi+pca(j,i)*(xi_teo(j)-xi_data(j))/err_data(j)
enddo
chi2=chi2+(par_chi**2)/ev(i) ! (****)
enddo
*/
/* Add in the error on the galaxy density
*/
if(HOD.free[0])
chi2+=chi2ngal;
t1 = clock();
t0 = difftime(t1,t0)/CLOCKS_PER_SEC;
niter++;
if(!ThisTask){
printf("ITER %7d %e ",niter,chi2);
for(i=1;i<=ncf_hod;++i)printf("%e ",a[i]);
printf(" %.2f\n",t0);
fflush(stdout);
if(HOD.free[0])
printf("NGAL %e %e %e\n",chi2ngal,
GALAXY_DENSITY,wp.ngal);
}
chi2_prev=chi2;
chi2_array[ichi]=chi2;
ichi++;
if(ichi==10)ichi=0;
fflush(stdout);
return(chi2);
}
void initial_wp_values(double *a, double **pp, double *yy)
{
static int flag=1;
int i,j;
double d[100];
if(flag) {
i=0;j=0;
if(HOD.free[++i])a[++j]=HOD.M_min;
if(HOD.free[++i])a[++j]=HOD.M1;
if(HOD.free[++i])a[++j]=HOD.alpha;
if(HOD.free[++i])a[++j]=HOD.M_cut;
if(HOD.free[++i])a[++j]=HOD.sigma_logM;
if(HOD.free[++i])a[++j]=CVIR_FAC;
if(HOD.pdfc>=7){
if(HOD.free[++i])a[++j]=HOD.M_cen_max; }
else {
if(HOD.free[++i])a[++j]=HOD.MaxCen; }
if(HOD.free[++i])a[++j]=HOD.M_sat_break;
if(HOD.free[++i])a[++j]=HOD.alpha1;
if(XCORR){
i=0;
if(HOD.free[++i])a[++j]=HOD2.M_min;
if(HOD.free[++i])a[++j]=HOD2.M1;
if(HOD.free[++i])a[++j]=HOD2.alpha;
if(HOD.free[++i])a[++j]=HOD2.M_cut;
if(HOD.free[++i])a[++j]=HOD2.sigma_logM;
if(HOD.free[++i])a[++j]=CVIR_FAC;
if(HOD2.pdfc>=7){
if(HOD.free[++i])a[++j]=HOD2.M_cen_max; }
else {
if(HOD.free[++i])a[++j]=HOD2.MaxCen; }
}
printf("INITIAL VALUES: ");
for(i=1;i<=wp.ncf;++i)printf("%e ",a[i]);
printf("\n");
}
//flag++;
/* Make the starting stepsize 10% of the initial values.
*/
for(i=1;i<=wp.ncf;++i)
d[i]=a[i]*0.25/flag;
if(POWELL)
{
for(i=1;i<=wp.ncf;++i)
{
for(j=1;j<=wp.ncf;++j)
{
pp[i][j]=0;
if(i==j)pp[i][j]+=d[j];
}
}
}
else
{
for(j=1;j<=wp.ncf;++j)
pp[1][j]=a[j];
yy[1]=chi2_wp(a);
for(i=1;i<=wp.ncf;++i)
{
a[i]+=d[i];
if(i>1)a[i-1]-=d[i-1];
yy[i+1]=chi2_wp(a);
for(j=1;j<=wp.ncf;++j)
pp[i+1][j]=a[j];
}
a[wp.ncf]-=d[wp.ncf];
}
}
/* This routine reads in the wp data and covariance matrix from
* the filenames specified.
*
* FILE FORMATS:
*
* - fname_wp -> r xi e_xi
* - fname_covar -> (i=1,np)(j=1,np) read(covar[i][j])
*
*/
void wp_input()
{
float x1,x2,x3;
FILE *fp;
int i,j,n;
char a[1000];
if(!(fp=fopen(wp.fname_wp,"r")))
{
fprintf(stdout,"ERROR opening [%s]\n",wp.fname_wp);
endrun("error in wp_input");
}
wp.np=filesize(fp);
/* [wp.format==2] means that there are two header lines at
* the top of the file.
*/
if(wp.format==2)
{
wp.np-=2;
fgets(a,1000,fp);
fgets(a,1000,fp);
}
/* [wp.format==3] means that there is one header lines at
* the top of the file.
*/
if(wp.format==3)
{
wp.np-=1;
fgets(a,1000,fp);
}
wp.r=dvector(1,wp.np);
wp.x=dvector(1,wp.np);
wp.e=dvector(1,wp.np);
if(PCA)
{
wp.eigen=dvector(1,wp.np);
wp.covar=dmatrix(1,wp.np,1,wp.np);
}
/* Read in the projected correlation function data.
* Standard format [wp.format==1] is linear r, linear wp, linear err.
* [wp.format==2] is log10 r, log10 wp, linear err.
* NB! Peder's format is to list the inner edge of the bin, so we're adding 0.1 to each number.
*/
for(i=1;i<=wp.np;++i)
{
fscanf(fp,"%f %f %f",&x1,&x2,&x3);
wp.r[i]=x1;
wp.x[i]=x2;
wp.e[i]=x3;
if(wp.format==2){
wp.r[i] = pow(10.0,wp.r[i]+0.1);
wp.x[i] = pow(10.0,wp.x[i])*wp.r[i];
wp.e[i] = wp.e[i]*wp.r[i];
}
if(wp.format==3){
fscanf(fp,"%f",&x3);
wp.e[i]=x3;
wp.r[i] = pow(10.0,wp.r[i]+0.1);
// wp.x[i] = wp.x[i]*wp.r[i];
// wp.e[i] = wp.e[i]*wp.r[i];
}
if(wp.format==3 && PCA) {
fscanf(fp,"%lf",&wp.eigen[i]);
for(j=1;j<=wp.np;++j)
fscanf(fp,"%lf",&wp.covar[j][i]);
if(wp.npca==0)
wp.npca = wp.np;
}
// if(wp.format==3 && !PCA)
fgets(a,1000,fp);
}
fclose(fp);
fprintf(stderr,"Done reading %d lines from [%s]\n",wp.np,wp.fname_wp);
if(!COVAR || PCA)
return;
/*
if(wp.format==1)
{
Work.SDSS_bins=1;
for(i=1;i<=40;++i)
Work.rad[i] = i-0.5;
}
*/
if(!(fp=fopen(wp.fname_covar,"r")))
{
fprintf(stdout,"ERROR opening [%s]\n",wp.fname_covar);
endrun("error in wp_input");
}
wp.covar=dmatrix(1,wp.np,1,wp.np);
for(i=1;i<=wp.np;++i)
for(j=1;j<=wp.np;++j)
{
fscanf(fp,"%lf",&(wp.covar[i][j]));
/* printf("COVAR %d %d %e\n",i,j,wp.covar[i][j]); */
}
fclose(fp);
if(!ThisTask)
fprintf(stdout,"Done reading %d lines from [%s]\n",wp.np,wp.fname_covar);
}

180
c_tools/hod/xi_matter.c
View file

@ -1,180 +0,0 @@
#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);
}
}