Skip to content

Commit

Permalink
Merge branch 'jhod0-upstream-version'
Browse files Browse the repository at this point in the history
  • Loading branch information
@tmcclintock
tmcclintock committed Jun 18, 2019
2 parents 9d4a45a + 752addd commit 9dae82a
Show file tree
Hide file tree
Showing 7 changed files with 141 additions and 207 deletions.
18 changes: 8 additions & 10 deletions src/C_averaging.c
View file
Expand Up @@ -37,23 +37,22 @@ double ave_integrand(double lR, void*params){

int average_profile_in_bins(double*Redges, int Nedges, double*R, int NR,
double*profile, double*ave_profile){
int i;
gsl_spline*spline = gsl_spline_alloc(gsl_interp_cspline, NR);
gsl_interp_accel*acc= gsl_interp_accel_alloc();
gsl_integration_workspace * ws = gsl_integration_workspace_alloc(workspace_size);
integrand_params*params=malloc(sizeof(integrand_params));
gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, NR);
gsl_interp_accel *acc = gsl_interp_accel_alloc();
gsl_integration_workspace *ws = gsl_integration_workspace_alloc(workspace_size);
gsl_function F;
double result, err;

gsl_spline_init(spline, R, profile, NR);

params->acc = acc;
params->spline = spline;
F.params = params;
integrand_params params;
params.acc = acc;
params.spline = spline;
F.params = &params;
F.function = &ave_integrand;

//Loop over bins and compute the average
for(i = 0; i < Nedges-1; i++){
for(int i = 0; i < Nedges-1; i++){
gsl_integration_qag(&F, log(Redges[i]), log(Redges[i+1]), ABSERR, RELERR,
workspace_size, 6, ws, &result, &err);
ave_profile[i] = 2*result/(Redges[i+1]*Redges[i+1]-Redges[i]*Redges[i]);
Expand All @@ -63,6 +62,5 @@ int average_profile_in_bins(double*Redges, int Nedges, double*R, int NR,
gsl_spline_free(spline);
gsl_interp_accel_free(acc);
gsl_integration_workspace_free(ws);
free(params);
return 0;
}
29 changes: 14 additions & 15 deletions src/C_concentration.c
View file
Expand Up @@ -78,7 +78,7 @@ double DK15_concentration_at_Mmean(double Mass, double*k, double*Plin, int Nk, i
double n_s,double Omega_b, double Omega_m,
double h, double T_CMB){
int status;
mc_params*pars = (mc_params*)malloc(sizeof(mc_params));
mc_params pars;
double R = pow(Mass/(1.33333333333*M_PI*rhocrit*Omega_m*delta), 0.33333333); //R200m
double M_lo = Mass/10;
double M_hi = Mass*10;
Expand All @@ -88,20 +88,20 @@ double DK15_concentration_at_Mmean(double Mass, double*k, double*Plin, int Nk, i
gsl_root_fsolver*s = gsl_root_fsolver_alloc(T);
gsl_function F;
int iter = 0, max_iter = 100;
pars->Mm = Mass;
pars->Rm = R;
pars->k = k;
pars->Plin = Plin;
pars->Nk = Nk;
pars->delta = delta;
pars->h = h;
pars->n_s = n_s;
pars->Omega_b = Omega_b;
pars->Omega_m = Omega_m;
pars->T_CMB = T_CMB;
pars.Mm = Mass;
pars.Rm = R;
pars.k = k;
pars.Plin = Plin;
pars.Nk = Nk;
pars.delta = delta;
pars.h = h;
pars.n_s = n_s;
pars.Omega_b = Omega_b;
pars.Omega_m = Omega_m;
pars.T_CMB = T_CMB;

F.function = &Mm_from_Mc;
F.params = pars;
F.params = &pars;

status = gsl_root_fsolver_set(s, &F, M_lo, M_hi);

Expand All @@ -115,9 +115,8 @@ double DK15_concentration_at_Mmean(double Mass, double*k, double*Plin, int Nk, i
status = gsl_root_test_interval(M_lo, M_hi, 0, 0.001);
}while(status == GSL_CONTINUE && iter < max_iter);

cm = pars->cm;
cm = pars.cm;

free(pars);
gsl_root_fsolver_free(s);
return cm;
}
Expand Down
57 changes: 24 additions & 33 deletions src/C_deltasigma.c
View file
Expand Up @@ -36,14 +36,8 @@
* Note: all distances are comoving.
*/
double Sigma_nfw_at_R(double R, double M, double c, int delta, double om){
double*Rs = (double*)malloc(sizeof(double));
double*Sigma = (double*)malloc(sizeof(double));
double result;
Rs[0] = R;
Sigma_nfw_at_R_arr(Rs, 1, M, c, delta, om, Sigma);
result = Sigma[0];
free(Rs);
free(Sigma);
double result = 0.0;
Sigma_nfw_at_R_arr(&R, 1, M, c, delta, om, &result);
return result;
}

Expand Down Expand Up @@ -152,7 +146,7 @@ int Sigma_at_R_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double M, d

gsl_interp_accel*acc= gsl_interp_accel_alloc();
gsl_integration_workspace*workspace = gsl_integration_workspace_alloc(workspace_size);
integrand_params*params = malloc(sizeof(integrand_params));
integrand_params params;
gsl_function F;
double result1, err1, result2, err2;
int i;
Expand All @@ -161,17 +155,17 @@ int Sigma_at_R_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double M, d
lnRxi[i] = log(Rxi[i]);
}
gsl_spline_init(spline, lnRxi, xi, Nxi);
params->acc = acc;
params->spline = spline;
params->M = M;
params->conc= conc;
params->delta = delta;
params->om = om;
F.params = params;
params.acc = acc;
params.spline = spline;
params.M = M;
params.conc= conc;
params.delta = delta;
params.om = om;
F.params = &params;
int status;
for(i = 0; i < NR; i++){
ln_z_max = log(sqrt(Rxi_max*Rxi_max - R[i]*R[i])); //Max distance to integrate to
params->Rp = R[i];
params.Rp = R[i];
if(R[i] < Rxi0){
F.function = &integrand_small_scales;
status = gsl_integration_qag(&F, log(Rxi0)-10, log(sqrt(Rxi0*Rxi0-R[i]*R[i])), ABSERR, RELERR, workspace_size, KEY, workspace, &result1, &err1);
Expand All @@ -193,7 +187,6 @@ int Sigma_at_R_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double M, d
gsl_interp_accel_free(acc);
gsl_integration_workspace_free(workspace);
free(lnRxi);
free(params);
return 0;
}

Expand All @@ -213,22 +206,21 @@ int Sigma_at_R_full_arr(double*R, int NR, double*Rxi, double*xi, int Nxi, double
double ln_z_max;
double result3, err3;
gsl_integration_workspace*workspace = gsl_integration_workspace_alloc(workspace_size);
integrand_params*params=malloc(sizeof(integrand_params));
integrand_params params;
gsl_function F;
int i;
Sigma_at_R_arr(R, NR, Rxi, xi, Nxi, M, conc, delta, om, Sigma);
params->slope = log(xi[Nxi-1]/xi[Nxi-2])/log(Rxi[Nxi-1]/Rxi[Nxi-2]);
params->intercept = xi[Nxi-1]/pow(Rxi[Nxi-1], params->slope);
F.params=params;
params.slope = log(xi[Nxi-1]/xi[Nxi-2])/log(Rxi[Nxi-1]/Rxi[Nxi-2]);
params.intercept = xi[Nxi-1]/pow(Rxi[Nxi-1], params.slope);
F.params = &params;
F.function = &integrand_large_scales;
for(i = 0; i < NR; i++){
ln_z_max = log(sqrt(Rxi_max*Rxi_max - R[i]*R[i]));
params->Rp = R[i];
params.Rp = R[i];
gsl_integration_qag(&F, ln_z_max, ln_z_max+ulim, ABSERR, RELERR, workspace_size, KEY, workspace, &result3, &err3);
Sigma[i] += (result3*rhom*2);
}
gsl_integration_workspace_free(workspace);
free(params);
return 0;
}

Expand All @@ -255,7 +247,7 @@ int DeltaSigma_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, doubl
gsl_spline*spline = gsl_spline_alloc(gsl_interp_cspline, Ns);
gsl_interp_accel*acc = gsl_interp_accel_alloc();
gsl_integration_workspace* workspace = gsl_integration_workspace_alloc(workspace_size);
integrand_params*params=malloc(sizeof(integrand_params));
integrand_params params;
double result1, result2, err1, err2;
gsl_function F;
int i;
Expand All @@ -264,13 +256,13 @@ int DeltaSigma_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, doubl
lnRs[i] = log(Rs[i]);
}
gsl_spline_init(spline, lnRs, Sigma, Ns);
params->spline = spline;
params->acc = acc;
params->M = M;
params->conc = conc;
params->delta = delta;
params->om = om;
F.params = params;
params.spline = spline;
params.acc = acc;
params.M = M;
params.conc = conc;
params.delta = delta;
params.om = om;
F.params = &params;
F.function = &DS_integrand_small_scales;
gsl_integration_qag(&F, lrmin-10, lrmin, ABSERR, RELERR, workspace_size, KEY, workspace, &result1, &err1);
F.function = &DS_integrand_medium_scales;
Expand All @@ -282,6 +274,5 @@ int DeltaSigma_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns, doubl
gsl_interp_accel_free(acc);
gsl_integration_workspace_free(workspace);
free(lnRs);
free(params);
return 0;
}
87 changes: 42 additions & 45 deletions src/C_miscentering.c
View file
Expand Up @@ -111,38 +111,37 @@ int Sigma_mis_single_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns,
static gsl_spline*spline = NULL;
static gsl_interp_accel*acc = NULL;
static gsl_integration_workspace*workspace = NULL;
static integrand_params*params = NULL;
static integrand_params params;
if (init_flag == 0){
init_flag = 1;
spline = gsl_spline_alloc(gsl_interp_cspline, Ns);
acc = gsl_interp_accel_alloc();
workspace = gsl_integration_workspace_alloc(workspace_size);
params = malloc(sizeof(integrand_params));
}

gsl_spline_init(spline, lnRs, Sigma, Ns);

params->acc = acc;
params->spline = spline;
params->workspace = workspace;
params->M = M;
params->conc = conc;
params->delta = delta;
params->Omega_m = Omega_m;
params->Rmis = Rmis;
params->Rmis2 = Rmis*Rmis;
params->rmin = Rs[0];
params->rmax = Rs[Ns-1];
params->lrmin = log(Rs[0]);
params->lrmax = log(Rs[Ns-1]);
params.acc = acc;
params.spline = spline;
params.workspace = workspace;
params.M = M;
params.conc = conc;
params.delta = delta;
params.Omega_m = Omega_m;
params.Rmis = Rmis;
params.Rmis2 = Rmis*Rmis;
params.rmin = Rs[0];
params.rmax = Rs[Ns-1];
params.lrmin = log(Rs[0]);
params.lrmax = log(Rs[Ns-1]);

//Angular integral
F.function=&single_angular_integrand;
F.params=params;
F.function = &single_angular_integrand;
F.params = &params;

for(i = 0; i < NR; i++){
params->Rp = R[i];
params->Rp2 = R[i] * R[i];
params.Rp = R[i];
params.Rp2 = R[i] * R[i];
gsl_integration_qag(&F, 0, M_PI, ABSERR, RELERR, workspace_size,
KEY, workspace, &result, &err);
Sigma_mis[i] = result/M_PI;
Expand Down Expand Up @@ -272,32 +271,31 @@ int Sigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns,
static gsl_interp_accel*acc = NULL;
static gsl_integration_workspace*workspace = NULL;
static gsl_integration_workspace*workspace2 = NULL;
static integrand_params*params = NULL;
static integrand_params params;
if (init_flag == 0){
init_flag = 1;
spline = gsl_spline_alloc(gsl_interp_cspline, Ns);
acc = gsl_interp_accel_alloc();
workspace = gsl_integration_workspace_alloc(workspace_size);
workspace2 = gsl_integration_workspace_alloc(workspace_size);
params = malloc(sizeof(integrand_params));
}

gsl_spline_init(spline, lnRs, Sigma, Ns);

params->spline = spline;
params->acc = acc;
params->workspace = workspace;
params->workspace2 = workspace2;
params->M = M;
params->conc = conc;
params->delta = delta;
params->Omega_m = Omega_m;
params->Rmis = Rmis;
params->Rmis2 = Rmis*Rmis;
params->rmin = Rs[0];
params->rmax = Rs[Ns-1];
params->lrmin = log(Rs[0]);
params->lrmax = log(Rs[Ns-1]);
params.spline = spline;
params.acc = acc;
params.workspace = workspace;
params.workspace2 = workspace2;
params.M = M;
params.conc = conc;
params.delta = delta;
params.Omega_m = Omega_m;
params.Rmis = Rmis;
params.Rmis2 = Rmis*Rmis;
params.rmin = Rs[0];
params.rmax = Rs[Ns-1];
params.lrmin = log(Rs[0]);
params.lrmax = log(Rs[Ns-1]);

//Angular integral
F.function = &angular_integrand;
Expand All @@ -313,14 +311,14 @@ int Sigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma, int Ns,
}

//Assign the params struct to the GSL functions.
F_radial.params = params;
params->F_radial = F_radial;
F.params = params;
F_radial.params = &params;
params.F_radial = F_radial;
F.params = &params;

//Angular integral first
for(i = 0; i < NR; i++){
params->Rp = R[i];
params->Rp2 = R[i] * R[i]; //Optimization
params.Rp = R[i];
params.Rp2 = R[i] * R[i]; //Optimization

gsl_integration_qag(&F, 0, M_PI, ABSERR, RELERR, workspace_size,
KEY, workspace, &result, &err);
Expand Down Expand Up @@ -384,20 +382,19 @@ int DeltaSigma_mis_at_R_arr(double*R, int NR, double*Rs, double*Sigma_mis, int N
static gsl_spline*spline = NULL;
static gsl_interp_accel*acc = NULL;
static gsl_integration_workspace*workspace = NULL;
static integrand_params*params = NULL;
static integrand_params params;
if (init_flag == 0){
init_flag = 1;
spline = gsl_spline_alloc(gsl_interp_cspline, Ns);
acc = gsl_interp_accel_alloc();
workspace = gsl_integration_workspace_alloc(workspace_size);
params = malloc(sizeof(integrand_params));
}

gsl_spline_init(spline, Rs, Sigma_mis, Ns);

params->spline = spline;
params->acc = acc;
F.params = params;
params.spline = spline;
params.acc = acc;
F.params = &params;
F.function = &DS_mis_integrand;

for(i = 0; i < NR; i++){
Expand Down

0 comments on commit 9dae82a

Please sign in to comment.