Implementation of ROMBERG integration in background.c

[Physics-16]

Hi,

Actually I want to modify the file background.c through the implementation of ROMBERG integration for the continuity integral. The code for ROMBERG integration which I added in the file background.c; is the following (NOTE:- This ROMBERG code has been taken from this CLASS issue "https://github.com/lesgourg/class_public/pull/170"):-- ** integral_fld = 0.0;

        int  integration_successful = _FALSE_;
        int  max_steps = 33;
        int  i, j;
        double acc = 1.0e-5;
    double  R1[max_steps], R2[max_steps];           //buffers
        double  *Rp = &R1[0], *Rc = &R2[0];            //Rp is previous row, Rc is current row
        double  h_rom = (1.0 - a);                    //step size

        Rp[0] = 0.0;
        class_call(background_w_fld(pba,1.0,&w_fld,&dw_over_da_fld,&integral_fld), pba->error_message, pba->error_message);
        Rp[0] = 3.0 * (1.0 + w_fld) / 1.0 * h_rom *0.5;
        class_call(background_w_fld(pba,a,&w_fld,&dw_over_da_fld,&integral_fld), pba->error_message, pba->error_message);
        Rp[0] += 3.0 * (1.0 + w_fld) / a * h_rom *0.5;

      for(i = 1; i < max_steps; ++i){
     h_rom /= 2.0;
     double c = 0.0;
     size_t ep = 1 << (i-1);                //2^(n-1)
     for(j = 1; j <= ep; ++j){
        class_call(background_w_fld(pba,a+(2*j-1)*h_rom,&w_fld,&dw_over_da_fld,&integral_fld), pba->error_message,   pba->error_message);
        c += 3.0 * (1.0 + w_fld) / (a+(2*j-1)*h_rom);
     }
     Rc[0] = h_rom*c + 0.5*Rp[0];           //R(i,0)
     for(j = 1; j <= i; ++j){
       double n_k = pow(4, j);
       Rc[j] = (n_k*Rc[j-1] - Rp[j-1])/(n_k-1);     //compute R(i,j)
     }
     if(i > 1 && fabs(Rp[i-1]-Rc[i]) < acc){
        integral_fld = Rc[i-1];
        integration_successful = _TRUE_;
        break;
     }

  //swap Rn and Rc as we only need the last row
     double *rt = Rp;
     Rp = Rc;
     Rc = rt;
  }**  

Whenever I tried to run the new code with MontePython, I got the following error:--

Error in Class: thermodynamics_init(L:404) : error in thermodynamics_solve(ppr,pba,pth,ptw,pvecback) ;

=>thermodynamics_solve(L:1647) : error in thermodynamics_reionization_evolve_with_tau(&tpaw, interval_limit[index_interval], interval_limit[index_interval+1], mz_output, pth->tt_size) ;

=>thermodynamics_reionization_evolve_with_tau(L:2314) : condition (tau_sup < pth->tau_reio) is true; parameters are such that reionization cannot start after z_start_max ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ I think that the problem arises due to the inaccurate result produced by the ROMBERG code. But I can't found any error in the ROMBERG code. OR, It may be something else. Can anybody help me to reduce the error such that I can run MontePython with the presence of ROMBERG integration in CLASS ??