Skip to content

Commit

Permalink
added selection and evolution function read from file
Browse files Browse the repository at this point in the history
  • Loading branch information
Julien Lesgourgues committed Mar 5, 2014
1 parent 90e86ec commit 300a50c
Show file tree
Hide file tree
Showing 4 changed files with 265 additions and 6 deletions.
9 changes: 9 additions & 0 deletions include/transfer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,10 +59,19 @@ struct transfers {
short has_nz_file; /**< Has dN/dz (selection function) input file? */
short has_nz_analytic; /**< Use analytic form for dN/dz (selection function) distribution? */
FileName nz_file_name; /**< dN/dz (selection function) input file name */
int nz_size;
double * nz_z;
double * nz_nz;
double * nz_ddnz;

short has_nz_evo_file; /**< Has dN/dz (evolution function) input file? */
short has_nz_evo_analytic; /**< Use analytic form for dN/dz (evolution function) distribution? */
FileName nz_evo_file_name; /**< dN/dz (evolution function) input file name */
int nz_evo_size;
double * nz_evo_z;
double * nz_evo_nz;
double * nz_evo_dlog_nz;
double * nz_evo_dd_dlog_nz;

//@}

Expand Down
46 changes: 46 additions & 0 deletions myevolution.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
0.1 0.0010000000000000002
0.2 0.008
0.3 0.02699999999999939
0.4 0.0639999999999378
0.5 0.1249999999977896
0.6 0.21599999995913313
0.7 0.34299999951858506
0.8 0.5119999959225299
0.9 0.728999973156994
1. 0.999999855139168
1.1 1.330999334368509
1.2 1.7279973217531261
1.3 2.196990360740088
1.4 2.743968450165843
1.5 3.3749048508244637
1.6 4.095732787627785
1.7 4.912295140442309
1.8 5.830241082017937
1.9 6.854822885825836
2. 7.990512030628799
2.1 9.240299830454715
2.2 10.604466399010644
2.3 12.078486025087846
2.4 13.649587865208744
2.5 15.291333774900915
2.6 16.955497974925894
2.7 18.56075979355865
2.8 19.97872820320106
2.9 21.02045222024772
3. 21.431930578229572
3.1 20.915524906920037
3.2 19.201919576687313
3.3 16.190881291907615
3.4 12.133254618683612
3.5 7.729006674620789
3.6 3.942205156623325
3.7 1.4868118911900279
3.8 0.3732807590089654
3.9 0.05436027966326381
4. 0.003838922724865479
4.1 0.00010430820255980766
4.2 8.102007168749191e-7
4.3 1.2314496393552691e-9
4.4 2.2644902011221728e-13
4.5 2.7467569667648903e-18
4.6 1.0265616267989809e-24
11 changes: 11 additions & 0 deletions myselection.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
0.0 0.0
0.1 0.030591786603680816
0.2 0.10619451692325403
0.30000000000000004 0.19886750367043496
0.4 0.28446992920441827
0.5 0.34735619337280493
0.6 0.38083288182668446
0.7000000000000001 0.3853875240364537
0.8 0.3660955463257228
0.9 0.330122848124979
1. 0.2848017764321144
205 changes: 199 additions & 6 deletions source/transfer.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,11 @@ int transfer_init(
HyperInterpStruct BIS;
double xmax;

/* for reading selection function */
FILE * input_file;
int row,status;
double tmp1,tmp2;

/* This code can be optionally compiled with the openmp option for parallel computation.
Inside parallel regions, the use of the command "return" is forbidden.
For error management, instead of "return _FAILURE_", we will set the variable below
Expand Down Expand Up @@ -267,6 +272,102 @@ int transfer_init(
);
*/

/** - eventually read the selection and evolution functions */

ptr->nz_size = 0;

if (ptr->has_nz_file == _TRUE_) {

input_file = fopen(ptr->nz_file_name,"r");
class_test(input_file == NULL,
ptr->error_message,
"Could not open file %s!",ptr->nz_file_name);

/* Find size of table */
for (row=0,status=2; status==2; row++){
status = fscanf(input_file,"%lf %lf",&tmp1,&tmp2);
}
rewind(input_file);
ptr->nz_size = row-1;

/* Allocate room for interpolation table */
class_alloc(ptr->nz_z,sizeof(double)*ptr->nz_size,ptr->error_message);
class_alloc(ptr->nz_nz,sizeof(double)*ptr->nz_size,ptr->error_message);
class_alloc(ptr->nz_ddnz,sizeof(double)*ptr->nz_size,ptr->error_message);

for (row=0; row<ptr->nz_size; row++){
status = fscanf(input_file,"%lf %lf",
&ptr->nz_z[row],&ptr->nz_nz[row]);
//printf("%d: (z,dNdz) = (%g,%g)\n",row,ptr->nz_z[row],ptr->nz_nz[row]);
}
fclose(input_file);

/* Call spline interpolation: */
class_call(array_spline_table_lines(ptr->nz_z,
ptr->nz_size,
ptr->nz_nz,
1,
ptr->nz_ddnz,
_SPLINE_EST_DERIV_,
ptr->error_message),
ptr->error_message,
ptr->error_message);
}

ptr->nz_evo_size = 0;

if (ptr->has_nz_evo_file == _TRUE_) {

input_file = fopen(ptr->nz_evo_file_name,"r");
class_test(input_file == NULL,
ptr->error_message,
"Could not open file %s!",ptr->nz_evo_file_name);

/* Find size of table */
for (row=0,status=2; status==2; row++){
status = fscanf(input_file,"%lf %lf",&tmp1,&tmp2);
}
rewind(input_file);
ptr->nz_evo_size = row-1;

/* Allocate room for interpolation table */
class_alloc(ptr->nz_evo_z,sizeof(double)*ptr->nz_evo_size,ptr->error_message);
class_alloc(ptr->nz_evo_nz,sizeof(double)*ptr->nz_evo_size,ptr->error_message);
class_alloc(ptr->nz_evo_dlog_nz,sizeof(double)*ptr->nz_evo_size,ptr->error_message);
class_alloc(ptr->nz_evo_dd_dlog_nz,sizeof(double)*ptr->nz_evo_size,ptr->error_message);

for (row=0; row<ptr->nz_evo_size; row++){
status = fscanf(input_file,"%lf %lf",
&ptr->nz_evo_z[row],&ptr->nz_evo_nz[row]);
//printf("%d: (z,dNdz) = (%g,%g)\n",row,ptr->nz_evo_z[row],ptr->nz_evo_nz[row]);
}
fclose(input_file);

/* infer dlog(dN/dz)/dz from dN/dz */
ptr->nz_evo_dlog_nz[0] =
(ptr->nz_evo_dlog_nz[1]-ptr->nz_evo_dlog_nz[0])
/(ptr->nz_evo_z[1]-ptr->nz_evo_z[0]);
for (row=1; row<ptr->nz_evo_size-1; row++){
ptr->nz_evo_dlog_nz[row] =
(ptr->nz_evo_dlog_nz[row+1]-ptr->nz_evo_dlog_nz[row-1])
/(ptr->nz_evo_z[row+1]-ptr->nz_evo_z[row-1]);
}
ptr->nz_evo_dlog_nz[ptr->nz_evo_size-1] =
(ptr->nz_evo_dlog_nz[ptr->nz_evo_size-1]-ptr->nz_evo_dlog_nz[ptr->nz_evo_size-2])
/(ptr->nz_evo_z[ptr->nz_evo_size-1]-ptr->nz_evo_z[ptr->nz_evo_size-2]);

/* Call spline interpolation: */
class_call(array_spline_table_lines(ptr->nz_evo_z,
ptr->nz_evo_size,
ptr->nz_evo_dlog_nz,
1,
ptr->nz_evo_dd_dlog_nz,
_SPLINE_EST_DERIV_,
ptr->error_message),
ptr->error_message,
ptr->error_message);
}

/** (a.3.) workspace, allocated in a parallel zone since in openmp
version there is one workspace per thread */

Expand Down Expand Up @@ -408,6 +509,19 @@ int transfer_free(

}

if (ptr->nz_size > 0) {
free(ptr->nz_z);
free(ptr->nz_nz);
free(ptr->nz_ddnz);
}

if (ptr->nz_evo_size > 0) {
free(ptr->nz_evo_z);
free(ptr->nz_evo_nz);
free(ptr->nz_evo_dlog_nz);
free(ptr->nz_evo_dd_dlog_nz);
}

return _SUCCESS_;

}
Expand Down Expand Up @@ -2148,6 +2262,7 @@ int transfer_sources(
double f_evo = 0.;

/* when the selection function is multiplied by a function dNdz */
double z;
double dNdz;
double dln_dNdz_dz;

Expand Down Expand Up @@ -2369,14 +2484,35 @@ int transfer_sources(
f_evo = 2./pvecback[pba->index_bg_H]/pvecback[pba->index_bg_a]/tau0_minus_tau[index_tau]
+ pvecback[pba->index_bg_H_prime]/pvecback[pba->index_bg_H]/pvecback[pba->index_bg_H]/pvecback[pba->index_bg_a];

z = pba->a_today/pvecback[pba->index_bg_a]-1.;

if (ptr->has_nz_evo_file ==_TRUE_) {
class_stop(ptr->error_message,"Selection function in a file: not yet implemented");

class_test((z<ptr->nz_evo_z[0]) || (z>ptr->nz_evo_z[ptr->nz_evo_size-1]),
ptr->error_message,
"Your input file for the selection function only covers the redhsift range [%f : %f]. However, your input for the selection function requires z=%f",
ptr->nz_evo_z[0],
ptr->nz_evo_z[ptr->nz_evo_size-1],
z);

class_call(array_interpolate_spline(
ptr->nz_evo_z,
ptr->nz_evo_size,
ptr->nz_evo_dlog_nz,
ptr->nz_evo_dd_dlog_nz,
1,
z,
&last_index,
&dln_dNdz_dz,
1,
ptr->error_message),
ptr->error_message,
ptr->error_message);
}
else {

class_call(transfer_dNdz_analytic(ptr,
pba->a_today/pvecback[pba->index_bg_a]-1.,
z,
&dNdz,
&dln_dNdz_dz),
ptr->error_message,
Expand Down Expand Up @@ -2652,13 +2788,36 @@ int transfer_sources(
f_evo = 2./pvecback[pba->index_bg_H]/pvecback[pba->index_bg_a]/tau0_minus_tau[index_tau]
+ pvecback[pba->index_bg_H_prime]/pvecback[pba->index_bg_H]/pvecback[pba->index_bg_H]/pvecback[pba->index_bg_a];

z = pba->a_today/pvecback[pba->index_bg_a]-1.;

if (ptr->has_nz_evo_file == _TRUE_) {
class_stop(ptr->error_message,"Selection function in a file: not yet implemented");

class_test((z<ptr->nz_evo_z[0]) || (z>ptr->nz_evo_z[ptr->nz_evo_size-1]),
ptr->error_message,
"Your input file for the selection function only covers the redhsift range [%f : %f]. However, your input for the selection function requires z=%f",
ptr->nz_evo_z[0],
ptr->nz_evo_z[ptr->nz_evo_size-1],
z);

class_call(array_interpolate_spline(
ptr->nz_evo_z,
ptr->nz_evo_size,
ptr->nz_evo_dlog_nz,
ptr->nz_evo_dd_dlog_nz,
1,
z,
&last_index,
&dln_dNdz_dz,
1,
ptr->error_message),
ptr->error_message,
ptr->error_message);

}
else {

class_call(transfer_dNdz_analytic(ptr,
pba->a_today/pvecback[pba->index_bg_a]-1.,
z,
&dNdz,
&dln_dNdz_dz),
ptr->error_message,
Expand Down Expand Up @@ -2764,6 +2923,7 @@ int transfer_selection_function(
double x;
double dNdz;
double dln_dNdz_dz;
int last_index;

/* trivial dirac case */
if (ppt->selection==dirac) {
Expand All @@ -2789,7 +2949,27 @@ int transfer_selection_function(
if ((ptr->has_nz_file == _TRUE_) || (ptr->has_nz_analytic == _TRUE_)) {

if (ptr->has_nz_file == _TRUE_) {
class_stop(ptr->error_message,"Selection function in a file: not yet implemented");

class_test((z<ptr->nz_z[0]) || (z>ptr->nz_z[ptr->nz_size-1]),
ptr->error_message,
"Your input file for the selection function only covers the redhsift range [%f : %f]. However, your input for the selection function requires z=%f",
ptr->nz_z[0],
ptr->nz_z[ptr->nz_size-1],
z);

class_call(array_interpolate_spline(
ptr->nz_z,
ptr->nz_size,
ptr->nz_nz,
ptr->nz_ddnz,
1,
z,
&last_index,
&dNdz,
1,
ptr->error_message),
ptr->error_message,
ptr->error_message);
}
else {

Expand Down Expand Up @@ -2829,7 +3009,20 @@ int transfer_selection_function(
if ((ptr->has_nz_file == _TRUE_) || (ptr->has_nz_analytic == _TRUE_)) {

if (ptr->has_nz_file == _TRUE_) {
class_stop(ptr->error_message,"Selection function in a file: not yet implemented");

class_call(array_interpolate_spline(
ptr->nz_z,
ptr->nz_size,
ptr->nz_nz,
ptr->nz_ddnz,
1,
z,
&last_index,
&dNdz,
1,
ptr->error_message),
ptr->error_message,
ptr->error_message);
}
else {

Expand Down

0 comments on commit 300a50c

Please sign in to comment.