Example2:modify input.c new param

Modify input.c: new input parameter

Based on https://lesgourg.github.io/class-tour/Tokyo2014/lecture6_input_module.pdf

We will modify input.c so that hi_class accepts the baryon fraction (eta_b) as input parameter, apart from Omega_b (the baryon energy density ratio) and omega_b (omega_b * h^2 = Omega_b), which are the variables already defined.

In the referenced slides, they show that:

Omega_b * h^2 = 1.81*10^6 * eta_b * (T_gamma_0/K)^3

where T_gamma_0 is the CMB temperature.

Different procedures could be attemped. Our recommendation, which works without knowing which functions are defined in each module, is using grep. As we want to write an alternative to Omega_b, let's search where it is filled. Supposing you are in the source folder:

grep -n Omega_b input.c

(the -n is to say grep to print the line number) which outputs:

663:  class_call(parser_read_double(pfc,"Omega_b",&param1,&flag1,errmsg),
671:             "In input file, you can only enter one of Omega_b or omega_b, choose one");

Now we know where Omega_b is filled (line 663), using our favorite editor we open input.c:663 and see the interesting part of code is:

  /* Omega_0_b (baryons) */
  class_call(parser_read_double(pfc,"Omega_b",&param1,&flag1,errmsg),
             errmsg,
             errmsg);
  class_call(parser_read_double(pfc,"omega_b",&param2,&flag2,errmsg),
             errmsg,
             errmsg);
  class_test(((flag1 == _TRUE_) && (flag2 == _TRUE_)),
             errmsg,
             "In input file, you can only enter one of Omega_b or omega_b, choose one");
  if (flag1 == _TRUE_)
    pba->Omega0_b = param1;
  if (flag2 == _TRUE_)
    pba->Omega0_b = param2/pba->h/pba->h;

  Omega_tot += pba->Omega0_b;

To keep it clean, we will add our new parameter after the class_call that fills omega_b. Just copy and paste that line and modify for eta_b:

  class_call(parser_read_double(pfc,"eta_b",&param3,&flag3,errmsg),
             errmsg,
             errmsg);

Note: We can use param3 and flag3 because they are already initialized. In case a fourth variable were desired, we would have to initialized them at the beginning of the function (input_read_parameters) with the others (lines 521 and 522).

Then, we have to modify the class_test conditions, so that only one of the three variables is given as input:

  class_test(((flag1 == _TRUE_) && (flag2 == _TRUE_) && (flag3 == _TRUE_)),
             errmsg,
             "In input file, you can only enter one of Omega_b, omega_b or eta_b, choose one");

Finally, as CLASS works internally with pba->Omega0_b, we need to add a line to convert eta_b to Omega_b. Let's do it before the Omega_tot += pba->Omega0_b line:

  if (flag3 == _TRUE_)
    pba->Omega0_b = 1.81e6 * param3 * pow(pba->T_cmb ,3) / pba->h / pba->h;

Putting all together:

  /* Omega_0_b (baryons) */
  class_call(parser_read_double(pfc,"Omega_b",&param1,&flag1,errmsg),
             errmsg,
             errmsg);
  class_call(parser_read_double(pfc,"omega_b",&param2,&flag2,errmsg),
             errmsg,
             errmsg);
  class_call(parser_read_double(pfc,"eta_b",&param3,&flag3,errmsg),
             errmsg,
             errmsg);
  class_test(((flag1 == _TRUE_) && (flag2 == _TRUE_) && (flag3 == _TRUE_)),
             errmsg,
             "In input file, you can only enter one of Omega_b, omega_b or eta_b, choose one");
  if (flag1 == _TRUE_)
    pba->Omega0_b = param1;
  if (flag2 == _TRUE_)
    pba->Omega0_b = param2/pba->h/pba->h;
  if (flag3 == _TRUE_)
    pba->Omega0_b = 1.81e6 * param3 * pow(pba->T_cmb ,3) / pba->h / pba->h;

  Omega_tot += pba->Omega0_b;