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
vide_public/c_tools/hod/aspen_breakout.c

301 lines
6.6 KiB
C
Raw Blame History

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#ifdef PARALLEL
#include <mpi.h>
#endif
#include "header.h"
void output_hod(char fname[]);
void output_wp(char fname[]);
double func_find_mone(double m);
double func_find_mone2(double m);
double func_find_alpha(double a);
double fsat_g1;
void aspen_breakout()
{
int i,j,k,i1;
double x1,x2,x3,x4,x5,x6,r,dr,fsat,m1,mcut,mmin,dlogm;
char fname[100];
FILE *fp;
//goto loop_m1;
// goto loop_cvir;
//goto loop_alpha;
//goto loop_mcut2;
/* Vary M_min, find M1 such that fsat is fixed.
*/
OUTPUT=3;
HOD.alpha = 1.0;
for(i=1;i<=5;++i)
{
muh(0);
RESET_FLAG_1H++;
RESET_FLAG_2H++;
RESET_KAISER++;
sprintf(Task.root_filename,"fsat.%d",i);
HOD.M_min = pow(10.0,i+9.0);
HOD.M_cut = 4*HOD.M_min;
if(i==1) {
HOD.M1 = 30*HOD.M_min;
set_HOD_params();
fsat = fsat_g1 = qromo(func_satfrac,log(HOD.M_min),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
fprintf(stderr,"%e %e %e %f %f\n",HOD.M_min,HOD.M1,GALAXY_DENSITY,fsat,HOD.M1/HOD.M_min);
muh(2);
}
set_HOD_params();
HOD.M1 = exp(zbrent(func_find_mone2,log(HOD.M_low/10),log(HOD.M_max),1.0E-4));
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
fprintf(stdout,"MUH %e %e %e %f %f\n",HOD.M_min,HOD.M1,GALAXY_DENSITY,fsat,HOD.M1/HOD.M_min);
sprintf(fname,"fsat_1h.%d",i);
output_wp(fname);
sprintf(fname,"fsat_wp.%d",i);
output_wp(fname);
sprintf(fname,"fsat_hod.%d",i);
output_hod(fname);
}
exit(0);
/* Vary sigma8 at fixed HOD/ngal (so change mmin)
*/
for(i=0;i<=3;++i)
{
SIGMA_8 = 0.9*growthfactor((double)i);
fprintf(stderr,"z=%d, sigma8= %f\n",i,SIGMA_8);
RESET_COSMOLOGY++;
RESET_FLAG_1H++;
RESET_FLAG_2H++;
RESET_KAISER++;
//HOD.M_min = 0;
set_HOD_params();
sprintf(fname,"SIGMA8a.%d",i);
output_wp(fname);
}
exit(0);
SIGMA_8=0.9;
RESET_COSMOLOGY++;
loop_m1:
/* Vary fsat by varying M1
*/
m1 = HOD.M1;
i1 = 0;
for(i=-2;i<=2;++i)
{
HOD.M1 = m1*pow(10.0,i/10.0);
RESET_FLAG_1H++;
RESET_FLAG_2H++;
RESET_KAISER++;
HOD.M_min = HOD.M_low = 0;
set_HOD_params();
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
fprintf(stderr,"M1= %e fsat=%f\n",HOD.M1,fsat);
sprintf(fname,"M1_wp.%d",++i1);
output_wp(fname);
sprintf(fname,"M1_hod.%d",i1);
output_hod(fname);
}
HOD.M1 = m1;
exit(0);
loop_mcut1:
/* Vary M_cut, fix M1 and fsat
*/
mcut = HOD.M_cut;
HOD.M_min = 0;
set_HOD_params();
dlogm = (log(HOD.M1) - log(HOD.M_min))/4;
mmin = log(HOD.M_min);
fsat_g1 = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
i1 = 0;
for(i=0;i<5;++i)
{
HOD.M_cut = exp(dlogm*i+mmin);
HOD.M1 = exp(zbrent(func_find_mone,log(HOD.M_min),log(HOD.M_max),1.0E-4));
RESET_FLAG_1H++;
RESET_FLAG_2H++;
RESET_KAISER++;
HOD.M_min = HOD.M_low = 0;
set_HOD_params();
fprintf(stderr,"M_cut= %e M1= %e %f\n",HOD.M_cut,HOD.M1,HOD.M1/HOD.M_cut);
sprintf(fname,"Mcut1_wp.%d",++i1);
output_wp(fname);
sprintf(fname,"Mcut1_hod.%d",i1);
output_hod(fname);
}
loop_mcut2:
/* Vary M_cut/fsat, keep M_cut/M1 = 1
*/
mcut = 3.0e13;
m1 = 3.0e13;
i1 = 0;
dlogm = (log(3.0e13) - log(10.0e12))/4;
for(i=0;i<5;++i)
{
HOD.M_cut = HOD.M1 = exp(dlogm*i)*10.0e12;
RESET_FLAG_1H++;
RESET_FLAG_2H++;
RESET_KAISER++;
HOD.M_min = HOD.M_low = 0;
set_HOD_params();
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
fprintf(stderr,"M_min= %e M1= %e fsat=%f\n",HOD.M_min,HOD.M1,fsat);
sprintf(fname,"Mcut2_wp.%d",++i1);
output_wp(fname);
sprintf(fname,"Mcut2_hod.%d",i1);
output_hod(fname);
}
loop_alpha:
/* Vary Mcut as above, but fix f_sat by varying alpha
*/
mcut = 3.0e13;
m1 = 3.0e13;
i1 = 0;
mmin = log(6.0e12);
dlogm = (log(3.0e13) - mmin)/4;
HOD.M_cut = HOD.M1 = exp(mmin);
HOD.alpha = 0.1;
HOD.M_min = HOD.M_low = 0;
set_HOD_params();
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
fprintf(stderr,"fsat = %f %e\n",fsat,HOD.M_min);
for(i=0;i<5;++i)
{
HOD.M_cut = HOD.M1 = exp(dlogm*i + mmin);
RESET_FLAG_1H++;
RESET_FLAG_2H++;
RESET_KAISER++;
HOD.alpha = zbrent(func_find_alpha,0.05,2.0,1.0E-4);
fprintf(stderr,"M_cut= %e alpha = %f\n",HOD.M_cut,HOD.alpha);
sprintf(fname,"alpha_wp.%d",++i1);
output_wp(fname);
sprintf(fname,"alpha_hod.%d",i1);
output_hod(fname);
}
loop_cvir:
/* Vary cvir with central one above
*/
HOD.M_cut = HOD.M1 = exp(dlogm*2 + mmin);
HOD.alpha = zbrent(func_find_alpha,0.05,2.0,1.0E-4);
fprintf(stderr,"M_cut= %e alpha = %f\n",HOD.M_cut,HOD.alpha);
i1 = 0;
for(i=0;i<5;++i)
{
CVIR_FAC = pow(3.0,i-2);
RESET_FLAG_1H++;
RESET_FLAG_2H++;
RESET_KAISER++;
sprintf(fname,"cvir_wp.%d",++i1);
output_wp(fname);
sprintf(fname,"cvir_hod.%d",i1);
output_hod(fname);
}
exit(0);
}
double func_find_alpha(double a)
{
HOD.alpha = a;
return qromo(func_galaxy_density,log(HOD.M_low),log(HOD.M_max),midpnt) - GALAXY_DENSITY;
}
double func_find_mone(double m)
{
double fsat;
HOD.M1=exp(m);
HOD.M_min = 0;
set_HOD_params();
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
return fsat-fsat_g1;
}
double func_find_mone2(double m)
{
double fsat;
HOD.M1=exp(m);
set_HOD_params();
fsat = qromo(func_satfrac,log(HOD.M_low),log(HOD.M_max),midpnt)/GALAXY_DENSITY;
return fsat-fsat_g1;
}
void output_hod(char fname[])
{
FILE *fp;
double dlogm,sig,m;
int i;
fp=fopen(fname,"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;
sig = N_sat(m)*N_sat(m) + N_cen(m)*(1-N_cen(m));
fprintf(fp,"%e %e %e %e %e %e\n",
m,N_cen(m),N_sat(m),N_avg(m),sig,sig/(N_avg(m)*N_avg(m)));
}
fclose(fp);
}
void output_wp(char fname[])
{
FILE *fp;
double r,dr,x1,x2,x3,x4,x5;
int j;
fp=fopen(fname,"w");
dr=(log(70.0)-log(0.05))/49.0;
for(j=0;j<50;++j)
{
r=exp(j*dr+log(0.05));
x1=one_halo_real_space(r);
x2=two_halo_real_space(r);
x3=projected_xi(r);
//x4 = projected_xi1h(r);
//x5 = projected_xi2h(r);
fprintf(fp,"%f %e %e %e %e\n",r,x1,x2,x1+x2,x3);
fflush(fp);
}
fclose(fp);
}
Reference in a new issue View git blame Copy permalink