Merge branch 'jhod0-upstream-version'

Author: tmcclintock

src/C_averaging.c

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

src/C_concentration.c

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

src/C_deltasigma.c

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

src/C_miscentering.c

@@ -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;
@@ -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;
@@ -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);
@@ -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++){