Project initialization

Makefile

+LDFLAGS       =  -L ${LIB_CFITSIO} -L ${LIB_GSL}  -L ${LIB_MPI}  -lcfitsio   -lm   -lgsl   -lgslcblas -lmpi
+INCLUDE       =  -I ${PWD}  -I ${HEADERS_CFITSIO}  -I ${HEADERS_GSL}  -I ${HEADERS_MPI} -I ${MC_INSTDIR}/../../include/
+
+CFLAGS   = -g  -Ofast -Wall   ${INCLUDE}
+
+
+LOCAL_DIR               = ${PWD}
+MC_SLAB                 = ${MC_INSTDIR}/../../gr_code/slab_rte
+MC_POLARIZATION   	= ${MC_INSTDIR}/polarization
+MC_ELE_DISTRIBUTION     = ${MC_INSTDIR}/electron_distributions
+MC_GENERAL 		= ${MC_INSTDIR}/general_routines/
+MC_COORD   		= ${MC_INSTDIR}/coordinates_systems
+MC_SEED    		= ${MC_INSTDIR}/seed_photons/
+MC_HYDRO                = ${MC_INSTDIR}/../../gr_code/hydrodynamics
+
+
+#========================================================================
+
+
+LOCAL_DIR_OBJ=  ${LOCAL_DIR}/main_program.o   ${LOCAL_DIR}/check_results.o  ${LOCAL_DIR}/compute_results.o ${LOCAL_DIR}/function_integrals.o ${LOCAL_DIR}/legendre_integration.o ${LOCAL_DIR}/set_deltatau.o ${LOCAL_DIR}/slab_polarization.o ${LOCAL_DIR}/solve_rte.o ${LOCAL_DIR}/mc_slab.o ${LOCAL_DIR}/compute_stokes.o ${LOCAL_DIR}/collect_photons.o  
+
+MC_POLARIZATION_OBJ=${MC_POLARIZATION}/ortonormal_axes.o ${MC_POLARIZATION}/unpolarized_seedphot.o ${MC_POLARIZATION}/polarization_probability.o ${MC_POLARIZATION}/depolarization_probability.o ${MC_POLARIZATION}/electfield_as_new.o ${MC_POLARIZATION}/sample_azim_angle.o ${MC_POLARIZATION}/allocate_stokes_pointer.o ${MC_POLARIZATION}/disk_polarizationelectricfield.o ${MC_POLARIZATION}/pdeg_atmosphere.o ${MC_POLARIZATION}/write_polardata.o ${MC_POLARIZATION}/select_quadrant.o $ ${MC_POLARIZATION}/detector_axis.o
+
+MC_ELE_DISTRIBUTION_OBJ=${MC_ELE_DISTRIBUTION}/hybrid/eledist_hybrid.o ${MC_ELE_DISTRIBUTION}/hybrid/integrate_maxwgamma.o ${MC_ELE_DISTRIBUTION}/sample_electron.o ${MC_ELE_DISTRIBUTION}/powerlaw/eledist_nth.o ${MC_ELE_DISTRIBUTION}/maxwellian/subrel_maxwellian.o
+
+MC_COORD_OBJ=${MC_COORD}/versor_boost.o ${MC_COORD}/elmag_lorentztransform.o 
+
+MC_GENERAL_OBJ=${MC_INSTDIR}/../../gr_code/clock_random.o ${MC_INSTDIR}/../../gr_code/dmatrix.o ${MC_INSTDIR}/../../gr_code/round_to_digits.o ${MC_INSTDIR}/general_routines/gauleg.o ${MC_INSTDIR}/general_routines/numpar.o ${MC_INSTDIR}/general_routines/stringpar.o ${MC_INSTDIR}/general_routines/strrev.o ${MC_INSTDIR}/general_routines/stringconcat.o ${MC_INSTDIR}/general_routines/interpolate_kn.o ${MC_INSTDIR}/../../gr_code/dot_prod.o  ${MC_GENERAL}/read_array.o ${MC_INSTDIR}/../../gr_code/vect_prod.o ${MC_GENERAL}/read_matrix.o ${MC_GENERAL}/energy_lorentz_boost.o ${MC_GENERAL}/new_energy_erf.o ${MC_GENERAL}/make_array_angles.o ${MC_INSTDIR}/../../gr_code/init_mpi.o ${MC_INSTDIR}/../../gr_code/abs_val.o ${MC_GENERAL}/make_array_energy.o
+
+MC_SEED_OBJ=${MC_SEED}/bb_spec.o
+
+MC_HYDRO_OBJ=${MC_HYDRO}/reverse_array.o
+
+
+MC_SLAB_OBJ=${MC_SLAB}/electricfield_fromfile.o ${MC_SLAB}/read_resultsfile.o  ${MC_SLAB}/setup_gsl_objects.o ${MC_SLAB}/compute_cumulative_angdist.o 
+
+
+all:  start_iteration
+
+%.o: %.c $(DEPS)
+	$(CC) -c  -o  $@ $< $(CFLAGS)
+
+start_iteration:  ${LOCAL_DIR_OBJ} ${MC_GENERAL_OBJ} ${MC_POLARIZATION_OBJ} ${MC_ELE_DISTRIBUTION_OBJ} ${MC_COORD_OBJ} ${MC_SEED_OBJ}   ${MC_HYDRO_OBJ} ${MC_SLAB_OBJ}
+	${CC} $^ -o $@ ${LDFLAGS} 
+
+
+clean:
+	rm -f main_program.o ${LOCAL_DIR_OBJ} ${MC_GENERAL_OBJ} ${MC_POLARIZATION_OBJ} ${MC_ELE_DISTRIBUTION_OBJ} ${MC_COORD_OBJ} ${MC_SEED_OBJ}   ${MC_HYDRO_OBJ} ${MC_SLAB_OBJ}	
+
+
+ciao:
+	echo ${MC_SLAB_OBJ}
+
+

README.txt

+The code computes the polarization degree of radiation from a plane-parallel slab in the diffusion approximation for a pure scattering atmosphere.
+
+

check_results.c

+#include "iteration_functions.h"
+
+void check_results(gsl_spline* Spline_sample_Iu, gsl_spline* Spline_eval_pdeg, gsl_spline* Spline_eval_limbdark)
+{
+  double u_test = 0, P = 0, Itot = 0;
+  /*double Ir, Il;*/
+  double csi;
+  int ii, status;
+
+  for (ii = 0; ii < 2500; ii++)
+  {
+    csi = clock_random();
+    status = gsl_spline_eval_e(Spline_sample_Iu, csi, xacc_1d, &u_test);
+
+    if (status)
+    {
+      printf("File %s\n", __FILE__);
+      printf("csi=%lf \n", csi);
+      printf("Error at line %d: %s\n", __LINE__, gsl_strerror(status));
+      exit(1);
+    }
+
+    status = gsl_spline_eval_e(Spline_eval_pdeg, u_test, xacc_1d, &P);
+    if (status)
+    {
+      printf("File %s\n", __FILE__);
+      printf("csi=%lf u_test=%lf\n", csi, u_test);
+      printf("Error at line %d: %s\n", __LINE__, gsl_strerror(status));
+      exit(1);
+    }
+
+    status = gsl_spline_eval_e(Spline_eval_limbdark, u_test, xacc_1d, &Itot);
+    if (status)
+    {
+      printf("File %s\n", __FILE__);
+      printf("Error at line %d: %s\n", __LINE__, gsl_strerror(status));
+      exit(1);
+    }
+
+    if (P > 1)
+    {
+      printf("%lf %lf %5.4e\n", u_test, P, Itot);
+    }
+
+    /*  Ir = 1 / 2. * Itot * (1 + P);
+      Il = 1 / 2. * Itot * (1 - P);
+  */
+  }
+}
\ No newline at end of file

collect_photons.c

+#include "iteration_functions.h"
+
+#define size 100
+
+int good_j, good_k, n;
+
+/*==================================================================================*/
+
+void collect_photons(double ene, double* k, double* array_ene, double** array_counts, int flag)
+{
+  double center_ene;
+  double delta_ene;
+
+  double total_counts;
+  double counts;
+
+  int i, j;
+  double kz;
+
+  if (k[2] < 0)
+  {
+    kz = -k[2];
+  }
+  else
+  {
+    kz = k[2];
+  }
+
+  // =====================================================================#
+  // Define the array of the total counts
+  // =====================================================================#
+
+  static int array_total_counts[NSTEP_ENE];
+
+  // =====================================================================#
+  // Loop on energies
+  // =====================================================================#
+
+  for (i = 0; i < nstepangles; i++)
+  {
+    if (array_ubounds[i] < -1 || array_ubounds[i] > 1)
+    {
+      printf("Error bounds %d %lf \n", i, array_ubounds[i]);
+    }
+  }
+
+  for (i = 0; i < NSTEP_ENE; i++)
+  {
+    if ((ene >= array_ene[i]) && (ene < array_ene[i + 1]))
+    {
+      // =====================================================================#
+      // If the parameter nstepangles is not set, then the total
+      // spectrum is collected and exit from the loop over energy
+      // Otherwise make additional loop over the polar angle defined
+      // by array_u
+      // =====================================================================#
+
+      if (nstepangles < 2)
+      {
+        break;
+      }
+      else
+      {
+        for (j = 0; j < nstepangles; j++)
+        {
+          if (kz <= array_ubounds[j] && kz > array_ubounds[j + 1])
+          {
+            good_j = j;
+
+            array_counts[i][good_j]++;
+
+            break;
+          }
+        }
+      }
+
+      // Else condition for photon energy within the boundaries defined by array_ene
+    }
+    else
+    {
+      continue;
+    }
+
+    // End of loop over energy
+  }
+
+  if (flag == 1)
+  {
+    dat_spectrum = fopen(outspec, "w");
+
+    if (dat_spectrum == NULL)
+    {
+      printf("Error opening file %s\n", outspec);
+      exit(1);
+    }
+
+    // ===========================================================================
+    /*Compute the total spectrum summed over all angles*/
+    // ===========================================================================
+
+    for (i = 0; i < NSTEP_ENE; i++)
+    {
+      for (j = 0; j < nstepangles; j++)
+      {
+        array_total_counts[i] = array_total_counts[i] + array_counts[i][j];
+      }
+    }
+
+    // ====================================================================================
+    // Print spectra values, both total and, if set, also angle-dependent
+    // ====================================================================================
+
+    for (i = 0; i < NSTEP_ENE - 1; i++)
+    {
+      center_ene = 0.5 * (array_ene[i + 1] + array_ene[i]);
+      delta_ene = 0.5 * (array_ene[i + 1] - array_ene[i]);
+
+      // ====================================================================================
+      // Total counts per keV
+      // ====================================================================================
+
+      total_counts = array_total_counts[i] / (array_ene[i + 1] - array_ene[i]);
+
+      fprintf(dat_spectrum, "%5.3e  %5.3e  %5.3e ", center_ene, delta_ene, total_counts * center_ene);
+
+      for (j = 0; j < nstepangles; j++)
+      {
+        nph = 0;
+
+        // ====================================================================================##
+        // Print the columns with angle-dependent spectra spectrum in the file
+        // ====================================================================================##
+
+        counts = array_counts[i][j] / (array_ene[i + 1] - array_ene[i]);
+        nph = nph + array_counts[i][j];
+
+        if (j < nstepangles - 1)
+        {
+          fprintf(dat_spectrum, " %5.4e ", counts * center_ene);
+        }
+        else
+        {
+          fprintf(dat_spectrum, " %5.4e\n", counts * center_ene);
+        }
+      }
+    }
+
+    fclose(dat_spectrum);
+
+  } // End of flag condition for actually producing the qdp file
+
+  // End of program
+}

compute_results.c

+#include "iteration_functions.h"
+
+double sum_kiter(int k_min, int k_max, int jj, Iklr_intensity* ptr_Iklr, Ikl_intensity* ptr_Ikl, Ikr_intensity* ptr_Ikr, int method);
+
+void compute_results(int method,
+                     int k_iter,
+                     int Nstep_tau,
+                     int Nstep_mu,
+                     double* array_mu,
+                     double* weights_u,
+                     Iklr_intensity* ptr_Iklr,
+                     Ikl_intensity* ptr_Ikl,
+                     Ikr_intensity* ptr_Ikr)
+{
+  int jj, kk;
+  double* Pdeg = malloc(Nstep_mu * sizeof(double));
+
+  double* Pdeg_Band1 = malloc(Nstep_mu * sizeof(double));
+  double* Pdeg_Band2 = malloc(Nstep_mu * sizeof(double));
+  double* Pdeg_Band3 = malloc(Nstep_mu * sizeof(double));
+  double* Pdeg_Band4 = malloc(Nstep_mu * sizeof(double));
+  double* Pdeg_Band5 = malloc(Nstep_mu * sizeof(double));
+
+  double* Itot_upstream = malloc(Nstep_mu * sizeof(double));
+  double* Itot_downstream = malloc(Nstep_mu * sizeof(double));
+
+  double sum_Il, sum_Ir;
+
+  double sum_Il_downstream, sum_Ir_downstream;
+
+  FILE* dat;
+
+  dat = fopen(polarization_file, "w");
+
+  for (jj = 0; jj < Nstep_mu; jj++)
+  {
+    sum_Il = 0;
+    sum_Ir = 0;
+    sum_Il_downstream = 0;
+    sum_Ir_downstream = 0;
+
+    for (kk = 0; kk < k_iter; kk++)
+    {
+      if (method == 1)
+      {
+        sum_Il = sum_Il + ptr_Iklr[kk].Il_matrix_upstream[Nstep_tau - 1][jj];
+        sum_Ir = sum_Ir + ptr_Iklr[kk].Ir_matrix_upstream[Nstep_tau - 1][jj];
+
+        sum_Il_downstream = sum_Il_downstream + ptr_Iklr[kk].Il_matrix_downstream[Nstep_tau - 1][jj];
+        sum_Ir_downstream = sum_Ir_downstream + ptr_Iklr[kk].Ir_matrix_downstream[Nstep_tau - 1][jj];
+      }
+      else
+      {
+        sum_Il = sum_Il + ptr_Ikl[kk].matrix[0][jj];
+        sum_Ir = sum_Ir + ptr_Ikr[kk].matrix[0][jj];
+      }
+    }
+
+    if (disktau >= 4 && k_iter > 40)
+    {
+      Pdeg_Band1[jj] = sum_kiter(0, 10, jj, ptr_Iklr, ptr_Ikl, ptr_Ikr, method);
+      Pdeg_Band2[jj] = sum_kiter(10, 20, jj, ptr_Iklr, ptr_Ikl, ptr_Ikr, method);
+      Pdeg_Band3[jj] = sum_kiter(20, 30, jj, ptr_Iklr, ptr_Ikl, ptr_Ikr, method);
+      Pdeg_Band4[jj] = sum_kiter(30, 40, jj, ptr_Iklr, ptr_Ikl, ptr_Ikr, method);
+      Pdeg_Band5[jj] = sum_kiter(40, k_iter, jj, ptr_Iklr, ptr_Ikl, ptr_Ikr, method);
+    }
+
+    Pdeg[jj] = (sum_Ir - sum_Il) / (sum_Ir + sum_Il);
+
+    if (isnan(Pdeg[jj]))
+    {
+      printf("Error: Pdeg is NaN at u=%lf  sum_Ir=%4.3e sum_Il=%4.3e\n", array_mu[jj], sum_Ir, sum_Il);
+    }
+
+    Itot_upstream[jj] = (sum_Ir + sum_Il);
+    Itot_downstream[jj] = sum_Ir_downstream + sum_Il_downstream;
+  }
+
+  /*===========================================================================*/
+  /*Print I(u) normalized at the value u=1*/
+  /*===========================================================================*/
+
+  if (disktau < 4)
+  {
+    for (jj = 0; jj < Nstep_mu; jj++)
+    {
+      fprintf(dat, "%lf  %5.4e  %5.4e\n", array_mu[jj], Pdeg[jj],
+              Itot_upstream[jj] / Itot_upstream[Nstep_mu - 1]);
+    }
+  }
+  else
+  {
+    for (jj = 0; jj < Nstep_mu; jj++)
+    {
+      fprintf(dat, "%lf  %5.4e  %5.4e %5.4e  %5.4e  %5.4e  %5.4e  %5.4e \n", array_mu[jj], Pdeg[jj],
+              Itot_upstream[jj] / Itot_upstream[Nstep_mu - 1], Pdeg_Band1[jj], Pdeg_Band2[jj],
+              Pdeg_Band3[jj], Pdeg_Band4[jj], Pdeg_Band5[jj]);
+    }
+  }
+
+  fclose(dat);
+
+  /*===========================================================================*/
+
+  double TotalFlux = 0;
+  double TotalFlux_downstream = 0;
+
+  for (jj = 0; jj < Nstep_mu; jj++)
+  {
+    TotalFlux = TotalFlux + Itot_upstream[jj] * array_mu[jj] * weights_u[jj];
+
+    TotalFlux_downstream = TotalFlux_downstream + Itot_downstream[jj] * array_mu[jj] * weights_u[jj];
+  }
+
+  printf("Total flux at the top    (u>0)  %lf\n", TotalFlux);
+  printf("Total flux at the bottom (u<0)  %lf\n", TotalFlux_downstream);
+  printf("Percenatge of flux at the top %5.4f\n", TotalFlux / (TotalFlux + TotalFlux_downstream));
+
+  /*===========================================================================*/
+
+  dat = fopen(diffusion_file, "w");
+  fprintf(dat, "!Total flux at the top    (u>0)  %lf\n", TotalFlux);
+  fprintf(dat, "!Total flux at the bottom (u<0)  %lf\n", TotalFlux_downstream);
+  fprintf(dat, "!Percenatge of flux at the top %5.4f\n", TotalFlux / (TotalFlux + TotalFlux_downstream));
+
+  double* diffusion = malloc(k_iter * sizeof(double));
+  double sum_k;
+
+  for (kk = 0; kk < k_iter; kk++)
+  {
+    sum_k = 0;
+
+    for (jj = 0; jj < Nstep_mu; jj++)
+    {
+      if (method == 1)
+      {
+        sum_k = sum_k + ptr_Iklr[kk].Il_matrix_upstream[Nstep_tau - 1][jj] * weights_u[jj] * array_mu[jj] +
+                ptr_Iklr[kk].Ir_matrix_upstream[Nstep_tau - 1][jj] * weights_u[jj] * array_mu[jj];
+      }
+      else
+      {
+        sum_k = sum_k + ptr_Ikl[kk].matrix[0][jj] + ptr_Ikr[kk].matrix[0][jj];
+      }
+    }
+
+    diffusion[kk] = sum_k;
+    fprintf(dat, "%d  %5.4e\n", kk, diffusion[kk]);
+  }
+
+  fclose(dat);
+
+  /*===================================================================================*/
+
+  radiation_field_t* radiation_field = (radiation_field_t*)malloc(sizeof(radiation_field_t));
+  radiation_field->u_min = array_mu[0];
+  radiation_field->u_max = array_mu[Nstep_mu - 1];
+
+  read_resultsfile(polarization_file, radiation_field, RTE);
+
+  setup_gsl_objects(radiation_field, RTE);
+
+  check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);
+}
+
+/*========================================================================*/
+
+double sum_kiter(int k_min, int k_max, int jj, Iklr_intensity* ptr_Iklr, Ikl_intensity* ptr_Ikl, Ikr_intensity* ptr_Ikr, int method)
+{
+  int kk;
+  double Pdeg;
+  double sum_Il;
+  double sum_Ir;
+  double sum_Il_downstream;
+  double sum_Ir_downstream;
+
+  sum_Il = 0;
+  sum_Ir = 0;
+  sum_Il_downstream = 0;
+  sum_Ir_downstream = 0;
+
+  for (kk = k_min; kk < k_max; kk++)
+  {
+    if (method == 1)
+    {
+      sum_Il = sum_Il + ptr_Iklr[kk].Il_matrix_upstream[Nstep_tau - 1][jj];
+      sum_Ir = sum_Ir + ptr_Iklr[kk].Ir_matrix_upstream[Nstep_tau - 1][jj];
+
+      sum_Il_downstream = sum_Il_downstream + ptr_Iklr[kk].Il_matrix_downstream[Nstep_tau - 1][jj];
+      sum_Ir_downstream = sum_Ir_downstream + ptr_Iklr[kk].Ir_matrix_downstream[Nstep_tau - 1][jj];
+    }
+    else
+    {
+      sum_Il = sum_Il + ptr_Ikl[kk].matrix[0][jj];
+      sum_Ir = sum_Ir + ptr_Ikr[kk].matrix[0][jj];
+    }
+  }
+
+  Pdeg = (sum_Ir - sum_Il) / (sum_Ir + sum_Il);
+
+  return Pdeg;
+}

compute_stokes.c

+#include "iteration_functions.h"
+
+void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_parameters* ptr_stokes, double* array_ubounds, int flag)
+
+{
+  int FlagFound;
+  int ii, jj;
+  double Qs;
+  double Us;
+  double cos_theta;
+
+  static double Qp_cart[3];
+  static double Qm_cart[3];
+  static double Up_cart[3];
+  static double Um_cart[3];
+
+  cos_theta = k_lab[2];
+
+  /* =====================================================================# */
+  /* Loop on energies */
+  /* =====================================================================# */
+
+  FlagFound = FALSE;
+
+  if (flag == 0)
+  {
+    for (ii = 0; ii < POLAR_NBOUNDS - 1; ii++)
+    {
+      if ((ene >= ene_bound[ii]) && (ene <= ene_bound[ii + 1]))
+      {
+        /* =====================================================================# */
+        /* Loop on u array, the cosine of the polar angle. Indeed, k[2]=u */
+        /* =====================================================================# */
+
+        for (jj = 0; jj < nstepangles; jj++)
+        {
+          if ((cos_theta <= array_ubounds[jj]) && (cos_theta >= array_ubounds[jj + 1]))
+          {
+            /*	if (jj==nstepangles-1)
+                {
+                    printf("Angles [%5.4f - %5.4f]\n", array_ubounds[jj], array_ubounds[jj+1]);
+                }
+           */
+
+            detector_axis(k_lab, Qp_cart, Qm_cart, Up_cart, Um_cart);
+
+            Qs = pow(dot_prod(Qp_cart, polvect_cart), 2) - pow(dot_prod(Qm_cart, polvect_cart), 2);
+            Us = pow(dot_prod(Up_cart, polvect_cart), 2) - pow(dot_prod(Um_cart, polvect_cart), 2);
+
+            if (isnan(Qs))
+            {
+              printf("Error: Qs is Nan\n");
+              exit(1);
+            }
+
+            if (isnan(Us))
+            {
+              printf("Error: Us is Nan\n");
+              exit(1);
+            }
+
+            ptr_stokes[ii].array_Is[jj] = ptr_stokes[ii].array_Is[jj] + 1;
+            ptr_stokes[ii].array_Qs[jj] = ptr_stokes[ii].array_Qs[jj] + Qs;
+            ptr_stokes[ii].array_Us[jj] = ptr_stokes[ii].array_Us[jj] + Us;
+            ptr_stokes[ii].counter[jj]++;
+
+            FlagFound = TRUE;
+
+            break;
+          }
+        }
+      }
+    }
+
+    if (FlagFound == FALSE)
+    {
+      double zio = 1;
+      // printf("Warning: no correspondence for Ene=%5.4f k_z=%5.4f\n", ene, k_lab[2]);
+      // exit(1);
+    }
+  }
+  else
+  {
+    double u_center;
+    double Qtot;
+    double Utot;
+    double Itot_u1;
+    double Itot;
+    int Itot_allbands = 0;
+
+    double* Pdeg = (double*)malloc(nstepangles * sizeof(double));
+    dat_integral_polar = fopen(integral, "w");
+
+    for (jj = 0; jj < nstepangles; jj++)
+    {
+      u_center = (array_ubounds[jj] + array_ubounds[jj + 1]) / 2;
+
+      Qtot = 0;
+      Utot = 0;
+      Itot = 0;
+
+      for (ii = 0; ii < POLAR_NBOUNDS - 1; ii++)
+      {
+        Itot_allbands = Itot_allbands + ptr_stokes[ii].array_Is[jj];
+
+        Qtot = Qtot + ptr_stokes[ii].array_Qs[jj];
+        Utot = Utot + ptr_stokes[ii].array_Us[jj];
+        Itot = Itot + ptr_stokes[ii].array_Is[jj];
+
+        if (isnan(Qtot))
+        {
+          printf("Qtot is Nan at jj=%d ii=%d\n", jj, ii);
+          exit(1);
+        }
+
+        if (isnan(Utot))
+        {
+          printf("Utot is Nan at jj=%d ii=%d\n", jj, ii);
+          exit(1);
+        }
+
+        if (isnan((double)Itot))
+        {
+          printf("Itot is Nan at jj=%d ii=%d\n", jj, ii);
+          exit(1);
+        }
+      }
+
+      Pdeg[jj] = sqrt(Qtot * Qtot + Utot * Utot) / Itot;
+
+      if (jj == 0)
+      {
+        Itot_u1 = Itot;
+      }
+
+      fprintf(dat_integral_polar, "%lf %lf %lf %lf %5.4f\n", u_center, Pdeg[jj], Qtot / Itot,
+              Utot / Itot, 1 / u_center * (Itot / Itot_u1));
+    }
+
+    fclose(dat_integral_polar);
+
+    printf("Total counts all bands %d\n", Itot_allbands);
+  }
+}

dmatrix.c

+#include "iteration_functions.h"
+
+#define NR_END 1
+
+double** dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+  long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
+  double** m;
+
+  /* allocate pointers to rows */
+  m = (double**)malloc((size_t)((nrow + NR_END) * sizeof(double*)));
+
+  if (m == NULL)
+  {
+    printf("Error while allocating memory for 2D matrix, file %s\n", __FILE__);
+    exit(1);
+  }
+
+  m += NR_END;
+  m -= nrl;
+
+  /* allocate rows and set pointers to them */
+  m[nrl] = (double*)malloc((size_t)((nrow * ncol + NR_END) * sizeof(double)));
+  if (!m[nrl])
+  {
+    printf("Allocation failure 2 in matrix()");
+    exit(1);
+  }
+
+  m[nrl] += NR_END;
+  m[nrl] -= ncl;
+
+  for (i = nrl + 1; i <= nrh; i++)
+    m[i] = m[i - 1] + ncol;
+
+  /* return pointer to array of pointers to rows */
+  return m;
+}

function_integrals.c

+#include "iteration_functions.h"
+
+/*==================================================================*/
+
+double Ak_integral(double x, void* params)
+{
+  gsl_set_error_handler_off();
+
+  int seed_location = ((AngularFunctionParams*)params)->seed_location;
+  int k = ((AngularFunctionParams*)params)->k;
+  int status;
+  double tau = ((AngularFunctionParams*)params)->tau;
+  double tau_s = ((AngularFunctionParams*)params)->tau_s;
+
+  double tau_round_high;
+
+  double Ikl_a, Ikl_b, value;
+
+  if (k == 1)
+  {
+    if (seed_location == SEED_CENTER)
+    {
+      Ikl_a = I0_center(tau0, tau, x);
+      Ikl_b = I0_center(tau0, 2 * tau0 - tau, x);
+      // printf("Center: Tau %lf 2tau0-tau %5.4f I(tau)=%5.4e I(2tau0-tau)=%5.4e  (u=%lf) \n", tau,
+      //                        2 * tau0 - tau, Ikl_a, Ikl_b, x);
+    }
+    else if (seed_location == SEED_UNIFORM)
+    {
+      Ikl_a = I0_uniform(tau0, tau, x);
+      Ikl_b = I0_uniform(tau0, 2 * tau0 - tau, x);
+    }
+    else if (seed_location == SEED_DELTA)
+    {
+      Ikl_a = I0_delta(tau, tau_s, x);
+      Ikl_b = I0_delta(2 * tau0 - tau, 2 * tau0 - tau_s, x);
+
+      printf("Delta: Tau %lf 2tau0-tau %5.4f I(tau)=%5.4e I(2tau0-tau)=%5.4e  (u=%lf) \n", tau,
+             2 * tau0 - tau, Ikl_a, Ikl_b, x);
+    }
+    else if (seed_location == SEED_BASE)
+    {
+      Ikl_a = I0_base(tau, tau_s, x, U_POSITIVE);
+      Ikl_b = I0_base(2 * tau0 - tau, tau_s, x, U_NEGATIVE);
+
+      printf("Base: Tau %lf 2tau0-tau %5.4f I(tau)=%5.4e I(2tau0-tau)=%5.4e  (u=%lf) \n", tau,
+             2 * tau0 - tau, Ikl_a, Ikl_b, x);
+    }
+    else if (seed_location == SEED_EXP)
+    {
+      Ikl_a = I0_exp(tau, tau0, x);
+      Ikl_b = I0_exp(2 * tau0 - tau, tau0, x);
+    }
+
+    else
+    {
+      Ikl_a = 1;
+      Ikl_b = 1;
+    }
+  }
+  else
+  {
+    if (2 * tau0 - tau == 2 * tau0)
+    {
+      tau_round_high = 2 * tau0 - tau - 1e-7;
+    }
+    else
+    {
+      tau_round_high = 2 * tau0 - tau;
+    }
+
+    status = gsl_spline2d_eval_e(Spline_I2d_l, tau, x, xacc_2d, yacc_2d, &Ikl_a);
+
+    if (status)
+    {
+      printf("Source file %s\n", __FILE__);
+      printf("u=%lf tau=%lf\n", x, tau);
+      printf("Line %d Error: %s\n", __LINE__, gsl_strerror(status));
+      exit(1);
+    }
+
+    status = gsl_spline2d_eval_e(Spline_I2d_l, tau_round_high, x, xacc_2d, yacc_2d, &Ikl_b);
+
+    if (status)
+    {
+      printf("error: %s\n", gsl_strerror(status));
+    }
+  }
+
+  value = 3. / 4 * (1 - x * x) * (Ikl_a + Ikl_b);
+
+  return value;
+}
+
+/*==================================================================*/
+
+double Bk_integral(double x, void* params)
+{
+  int status;
+  int seed_location = ((AngularFunctionParams*)params)->seed_location;
+  int k = ((AngularFunctionParams*)params)->k;
+  double tau = ((AngularFunctionParams*)params)->tau;
+  double tau_s = ((AngularFunctionParams*)params)->tau_s;
+  double tau_round_high;
+
+  double Ikl_a, Ikl_b, value;
+  if (k == 1)
+  {
+    if (seed_location == SEED_CENTER)
+    {
+      Ikl_a = I0_center(tau0, tau, x);
+      Ikl_b = I0_center(tau0, 2 * tau0 - tau, x);
+    }
+    else if (seed_location == SEED_UNIFORM)
+    {
+      Ikl_a = I0_uniform(tau0, tau, x);
+      Ikl_b = I0_uniform(tau0, 2 * tau0 - tau, x);
+    }
+    else if (seed_location == SEED_DELTA)
+    {
+      Ikl_a = I0_delta(tau, tau_s, x);
+      Ikl_b = I0_delta(2 * tau0 - tau, 2 * tau0 - tau_s, x);
+    }
+    else if (seed_location == SEED_BASE)
+    {
+      Ikl_a = I0_base(tau, tau_s, x, U_POSITIVE);
+      Ikl_b = I0_base(2 * tau0 - tau, tau_s, x, U_NEGATIVE);
+    }
+    else if (seed_location == SEED_EXP)
+    {
+      Ikl_a = I0_exp(tau, tau0, x);
+      Ikl_b = I0_exp(2 * tau0 - tau, tau0, x);
+    }
+
+    else
+    {
+      Ikl_a = 1;
+      Ikl_b = 1;
+    }
+  }
+  else
+  {
+    if (2 * tau0 - tau == 2 * tau0)
+    {
+      tau_round_high = 2 * tau0 - tau - 1e-7;
+    }
+    else
+    {
+      tau_round_high = 2 * tau0 - tau;
+    }
+
+    status = gsl_spline2d_eval_e(Spline_I2d_l, tau, x, xacc_2d, yacc_2d, &Ikl_a);
+
+    if (status)
+    {
+      printf("error: %s\n", gsl_strerror(status));
+    }
+
+    status = gsl_spline2d_eval_e(Spline_I2d_l, tau_round_high, x, xacc_2d, yacc_2d, &Ikl_b);
+
+    if (status)
+    {
+      printf("error: %s\n", gsl_strerror(status));
+    }
+  }
+
+  value = 3. / 8 * pow(x, 2) * (Ikl_a + Ikl_b);
+
+  return value;
+}
+
+double Ck_integral(double x, void* params)
+{
+  int status;
+  int seed_location = ((AngularFunctionParams*)params)->seed_location;
+  int k = ((AngularFunctionParams*)params)->k;
+  double tau = ((AngularFunctionParams*)params)->tau;
+  double tau_s = ((AngularFunctionParams*)params)->tau_s;
+  double tau_round_high;
+
+  double Ikr_a, Ikr_b, value;
+
+  if (k == 1)
+  {
+    if (seed_location == SEED_CENTER)
+    {
+      Ikr_a = I0_center(tau0, tau, x);
+      Ikr_b = I0_center(tau0, 2 * tau0 - tau, x);
+    }
+    else if (seed_location == SEED_UNIFORM)
+    {
+      Ikr_a = I0_uniform(tau0, tau, x);
+      Ikr_b = I0_uniform(tau0, 2 * tau0 - tau, x);
+    }
+    else if (seed_location == SEED_DELTA)
+    {
+      Ikr_a = I0_delta(tau, tau_s, x);
+      Ikr_b = I0_delta(2 * tau0 - tau, 2 * tau0 - tau_s, x);
+    }
+    else if (seed_location == SEED_BASE)
+    {
+      Ikr_a = I0_base(tau, tau_s, x, U_POSITIVE);
+      Ikr_b = I0_base(2 * tau0 - tau, 2 * tau0 - tau_s, x, U_NEGATIVE);
+    }
+    else if (seed_location == SEED_EXP)
+    {
+      Ikr_a = I0_exp(tau, tau_s, x);
+      Ikr_b = I0_exp(2 * tau0 - tau, tau0, x);
+    }
+
+    else
+    {
+      Ikr_a = 1;
+      Ikr_b = 1;
+    }
+  }
+  else
+  {
+    if (2 * tau0 - tau == 2 * tau0)
+    {
+      tau_round_high = 2 * tau0 - tau - 1e-8;
+    }
+    else
+    {
+      tau_round_high = 2 * tau0 - tau;
+    }
+
+    status = gsl_spline2d_eval_e(Spline_I2d_r, tau, x, xacc_2d, yacc_2d, &Ikr_a);
+
+    if (status)
+    {
+      printf("error: %s\n", gsl_strerror(status));
+    }
+
+    status = gsl_spline2d_eval_e(Spline_I2d_r, tau_round_high, x, xacc_2d, yacc_2d, &Ikr_b);
+
+    if (status)
+    {
+      printf("error: %s\n", gsl_strerror(status));
+    }
+  }
+
+  value = 3. / 8 * (Ikr_a + Ikr_b);
+
+  return value;
+}
+
+/*==================================================================*/
+double Il_integral(double x, void* params)
+{
+  double Ak_interp, Bk_interp, Ck_interp;
+  double value, tau, mu;
+
+  tau = ((double*)params)[0];
+  mu = ((double*)params)[1];
+
+  Ak_interp = gsl_spline_eval(Spline_tau_Ak, x, xacc_1d);
+  Bk_interp = gsl_spline_eval(Spline_tau_Bk, x, xacc_1d);
+  Ck_interp = gsl_spline_eval(Spline_tau_Ck, x, xacc_1d);
+
+  value = exp(-(x - tau) / mu) * ((1 - mu * mu) * Ak_interp + mu * mu * (Bk_interp + Ck_interp)) * 1 / mu;
+
+  return value;
+}
+
+/*==================================================================*/
+double Ir_integral(double x, void* params)
+{
+  double Bk_interp, Ck_interp;
+  double value, tau, mu;
+
+  tau = ((double*)params)[0];
+  mu = ((double*)params)[1];
+
+  Bk_interp = gsl_spline_eval(Spline_tau_Bk, x, xacc_1d);
+  Ck_interp = gsl_spline_eval(Spline_tau_Ck, x, xacc_1d);
+
+  value = exp(-(x - tau) / mu) * (Bk_interp + Ck_interp) * 1 / mu;
+
+  return value;
+}

legendre_integration.c

+
+#include "iteration_functions.h"
+
+/*
+Perform integration of the function I(tau, u) over the range u=[0-1] using
+Gauss-Legendre quadrature
+*/
+
+/*========================================================*/
+/*Equation (4.206)*/
+/*Compute the integral*/
+/* [2*(1-u^2)*(1-u'^2) + u^2 u'^2] I_l(u)  du' */
+/*========================================================*/
+
+double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u)
+{
+  int status;
+  int ii;
+  double integ_value;
+  double angular_term;
+  double Ik_l;
+
+  integ_value = 0;
+  Ik_l = 0;
+
+  for (ii = 0; ii < Nstep_mu; ii++)
+  {
+    angular_term = scattering_matrix(1, 1, u, u_prime[ii]);
+
+    /*===================================================================*/
+    /*Evaluate by 2D-interpolation the function I^{k-1}_l (tau, u')*/
+    /*====================================================================*/
+
+    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l);
+
+    if (Ik_l < 0)
+    {
+      printf("Warning: interpolated Ik_l < 0 for tau=%5.4f u_prime=%5.4f\n", tau, u_prime[ii]);
+      exit(1);
+    }
+
+    if (status)
+    {
+      printf("Error at line %d\n", __LINE__);
+      printf("File %s\n", __FILE__);
+      printf("tau=%5.4f u=%5.4f\n", tau, u_prime[ii]);
+
+      exit(1);
+    }
+
+    // printf("ii=%d w %5.4e\n", ii, weights_u[ii]);
+
+    integ_value = integ_value + angular_term * weights_u[ii] * Ik_l;
+  }
+
+  if (isnan(integ_value))
+  {
+    printf("Error: values is NaN\n");
+    printf("Line %d\n", __LINE__);
+    exit(1);
+  }
+
+  // printf("Integral (A) %lf\n", integ_value);
+  return integ_value;
+}
+
+/*========================================================*/
+/*Equation (4.206)*/
+/*Compute the integral*/
+/* u^2 Ir(u)*/
+/*========================================================*/
+
+double legendre_integration_B(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u)
+{
+  int status;
+  int ii;
+  double integ_value;
+  double angular_term;
+  double Ik_r;
+
+  integ_value = 0;
+  Ik_r = 0;
+
+  for (ii = 0; ii < Nstep_mu; ii++)
+  {
+    angular_term = scattering_matrix(1, 2, u, u_prime[ii]);
+
+    /*===================================================================*/
+    /*Evaluate by 2D-interpolation the function I^{k-1}_r (tau, u')*/
+    /*====================================================================*/
+
+    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r);
+
+    if (status)
+    {
+      printf("Error at line %d\n", __LINE__);
+      printf("File %s\n", __FILE__);
+      exit(1);
+    }
+
+    integ_value = integ_value + angular_term * weights_u[ii] * Ik_r;
+  }
+
+  if (isnan(integ_value))
+  {
+    printf("Error: values is NaN\n");
+    printf("Line %d\n", __LINE__);
+    exit(1);
+  }
+
+  return integ_value;
+}
+
+/*========================================================*/
+/*Equation (4.207)*/
+/*Compute the integral*/
+/* u'^2 Il(u)*/
+/*========================================================*/
+
+double legendre_integration_C(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u)
+{
+  int status;
+  int ii;
+  double integ_value;
+  double angular_term;
+  double Ik_l;
+
+  integ_value = 0;
+  Ik_l = 0;
+
+  for (ii = 0; ii < Nstep_mu; ii++)
+  {
+    angular_term = scattering_matrix(2, 1, u, u_prime[ii]);
+
+    /*===================================================================*/
+    /*Evaluate by 2D-interpolation the function I^{k-1}_l (tau, u')*/
+    /*====================================================================*/
+
+    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l);
+
+    if (status)
+    {
+      printf("Error at line %d\n", __LINE__);
+      printf("File %s\n", __FILE__);
+      exit(1);
+    }
+
+    integ_value = integ_value + angular_term * weights_u[ii] * Ik_l;
+  }
+  if (isnan(integ_value))
+  {
+    printf("Error: values is NaN\n");
+    printf("Line %d\n", __LINE__);
+    exit(1);
+  }
+
+  return integ_value;
+}
+
+/*========================================================*/
+/*Equation (4.207)*/
+/*Compute the integral*/
+/* Ir(u)*/
+/*========================================================*/
+
+double legendre_integration_D(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u)
+{
+  int status;
+  int ii;
+  double integ_value;
+  double angular_term;
+  double Ik_r;
+
+  integ_value = 0;
+  Ik_r = 0;
+
+  for (ii = 0; ii < Nstep_mu; ii++)
+  {
+    angular_term = scattering_matrix(2, 2, u, u_prime[ii]);
+
+    /*===================================================================*/
+    /*Evaluate by 2D-interpolation the function I^{k-1}_r (tau, u')*/
+    /*====================================================================*/
+
+    status = gsl_spline2d_eval_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r);
+
+    if (status)
+    {
+      printf("Error at line %d\n", __LINE__);
+      printf("File %s\n", __FILE__);
+      exit(1);
+    }
+
+    integ_value = integ_value + angular_term * weights_u[ii] * Ik_r;
+  }
+
+  if (isnan(integ_value))
+  {
+    printf("Error: values is NaN\n");
+    printf("Line %d\n", __LINE__);
+    exit(1);
+  }
+
+  return integ_value;
+}
+
+/*================================================================================*/
+
+double scattering_matrix(int row, int col, double u, double u_prime)
+{
+  double value;
+
+  if (row == 1 && col == 1)
+  {
+    value = 2 * (1 - u * u) * (1 - u_prime * u_prime) + u * u * u_prime * u_prime;
+  }
+  else if (row == 1 && col == 2)
+  {
+    value = u * u;
+  }
+  else if (row == 2 && col == 1)
+  {
+    value = u_prime * u_prime;
+  }
+  else if (row == 2 && col == 2)
+  {
+    value = 1;
+  }
+  else
+  {
+    printf("Uncorrect row or column value for scattering matrix\n");
+    printf("row=%d col=%d\n", row, col);
+    exit(1);
+  }
+
+  return value;
+}

main_program.c

+#include "iteration_functions.h"
+#include <mpi.h>
+#include <stdio.h>
+
+int seed;
+double kte;
+double ktbb;
+double tau_c;
+double disktau;
+char* polarization_file;
+char* diffusion_file;
+char* polardisk;
+
+double gmin;
+double gmax;
+double p1;
+int FlagEleDist;
+int nph;
+int rank;
+//int nproc;
+int kiter;
+int DiskPolarizationFlag;
+
+char* outspec;
+char* polarfile;
+char* integral;
+char* diffusion;
+char* polarinput;
+
+char* set_filename(char* root, char* tau_char, char* seed_char);
+
+FILE* dat_spectrum;
+FILE* dat_diffusion;
+FILE* dat_integral_polar;
+FILE* dat_energy_polar;
+
+int main(int argc, char* argv[])
+{
+  int t;
+  int status;
+  char* method;
+  FlagEleDist = MAXW;
+
+  /*======================================================================*/
+  /*Initialize MPI*/
+  /*======================================================================*/
+
+  if (argc == 1)
+  {
+    printf(
+        "Program to compute with an iterative method the polarization degree and angular "
+        "distribution at the top of a plane-parallel scattering atmosphere\n");
+
+    printf("Synopsis\n");
+    printf("start_interation  disktau=disktau seed=[1,2,3] method=ode|st85|mc\n\n");
+
+    printf("Shared parameters\n\n");
+    printf("disktau=vaue: slab optical depth\n\n");
+
+    printf("\nSeed photons distributions: three options are possible\n");
+    printf("1: photons at the center of the slab\n");
+    printf("2: photons uniform distribution\n");
+    printf("3: photons at the base (valid only for ode and mc methods)\n\n");
+
+    printf("Parameters for iteration algorithms\n");
+
+    printf("method=ode  : solution of the ODE system for Il and Ir\n");
+    printf(
+        "method=st85 : use ST85 algorithm (valid only for seed photon at the center or uniformly "
+        "distributed)\n");
+
+    printf("kiter=N: number of iterations\n\n");
+
+    printf("Parameters for MC simulation\n");
+    printf("method=mc\n");
+    printf("nph=N       : number of simulated photons\n");
+    printf("kte=value   : electron temperature (default is 1 keV)\n");
+    printf("For seed photons at the bottom of the slab select if they are unpolarized\n");
+    printf(
+        "with option polardisk=n or the polarization and angular distribution is read from file\n");
+    printf("with polardisk=file\n");
+
+    return (1);
+  }
+
+  char* tau_char;
+  char* seed_char;
+
+  for (t = 1; t < argc; t++)
+  {
+    if (strstr(argv[t], "disktau=") != NULL)
+    {
+      disktau = numpar(argv[t]);
+      tau_char = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "kiter=") != NULL)
+      kiter = (int)numpar(argv[t]);
+
+    if (strstr(argv[t], "seed=") != NULL)
+    {
+      seed = (int)numpar(argv[t]);
+      seed_char = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "tau_c=") != NULL)
+      tau_c = numpar(argv[t]);
+
+    if (strstr(argv[t], "method=") != NULL)
+    {
+      method = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "nph=") != NULL)
+    {
+      nph = (int)numpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "polarfile=") != NULL)
+    {
+      polarfile = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "diffusion=") != NULL)
+    {
+      diffusion = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "outspec=") != NULL)
+    {
+      outspec = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "integral=") != NULL)
+    {
+      integral = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "polardisk=") != NULL)
+    {
+      polardisk = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "polarinput=") != NULL)
+    {
+      polarinput = stringpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "kte=") != NULL)
+    {
+      kte = numpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "ktbb=") != NULL)
+    {
+      ktbb = numpar(argv[t]);
+    }
+  }
+
+  /*=========================================================*/
+  /*Default values*/
+  /*=========================================================*/
+
+  if (method == NULL || method == 0x0)
+  {
+    printf("Please select the method to compute results with method=ode|st85|mc\n");
+    return 1;
+  }
+
+  if (diffusion == 0x0)
+  {
+    diffusion = "mc_diffusion.qdp";
+  }
+
+  if (integral == 0x0)
+  {
+    integral = "integral_polar_mc.qdp";
+  }
+
+  if (polarfile == 0x0)
+  {
+    polarfile = "energy_polar_mc.qdp";
+  }
+
+  if (outspec == 0x0)
+  {
+    outspec = "espectrum.qdp";
+  }
+
+  if (!kte)
+  {
+    kte = 1;
+  }
+
+  if (!ktbb)
+  {
+    ktbb = 1;
+  }
+
+  /*=========================================================*/
+
+  if (!disktau || disktau == 0)
+  {
+    printf("Please select the slab optical depth with parameter disktau=value\n");
+    return (0);
+  }
+
+  if (!seed)
+  {
+    printf("Please select the seed photon distribution with parameter seed=1,2,3\n\n");
+    return (1);
+  }
+
+  if (strcmp(method, "ode") == 0)
+  {
+    diffusion_file = set_filename("diffusion_ode_tau", tau_char, seed_char);
+    polarization_file = set_filename("intensity_ode_tau", tau_char, seed_char);
+
+    if (kiter == 0)
+    {
+      printf("Please select the number of iterations with parameter kiter=kiter\n");
+      exit(1);
+    }
+
+    solve_rte(kiter, seed);
+  }
+  else if (strcmp(method, "st85") == 0)
+  {
+    polarization_file = set_filename("intensity_st85_tau", tau_char, seed_char);
+    diffusion_file = set_filename("diffusion_st85_tau", tau_char, seed_char);
+
+    if (kiter == 0)
+    {
+      printf("Please select the number of iterations with parameter kiter=kiter\n");
+      exit(1);
+    }
+
+    slab_polarization(kiter, seed);
+  }
+  else if (strcmp(method, "mc") == 0)
+  {
+    MPI_Init(&argc, &argv); 
+    init_mpi(); 
+
+    if (polardisk == 0x0 && seed == 3 && rank == 0)
+    {
+      printf("Set polardisk=n or polardisk=file\n");
+      exit(1);
+    }
+
+    if (seed == 1 || seed == 2)
+    {
+      polardisk = "n";
+    }
+
+    if (nph == 0)
+    {
+      printf("Please select the number of photons for MC simulation with parametere nph=nph\n");
+      exit(1);
+    }
+
+    if (strcmp(polardisk, "n") == 0)
+    {
+      DiskPolarizationFlag = FALSE;
+    }
+    else if (strcmp(polardisk, "file") == 0)
+    {
+      DiskPolarizationFlag = FROMFILE;
+
+      if (polarinput == 0x0)
+      {
+        printf("Select the pre-tabulated reference file with parameter polarinput=filename\n");
+        exit(1);
+      }
+    }
+    else
+    {
+      printf("Unknonw option for parameter polardisk\n");
+      exit(1);
+    }
+
+    status = slab_mc(nph, seed);
+    // check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);
+  }
+
+  else
+  {
+    printf(
+        "\nPlease select method=ode (solve ODE), method=st85 (ST85) or method=mc (Montecarlo)\n");
+    return (1);
+  }
+
+  MPI_Finalize();
+
+  return (0);
+}
+
+char* set_filename(char* root, char* tau_char, char* seed_char)
+{
+  // root="intensity_ode";
+
+  char* filename;
+  char* tmp_a;
+  char* tmp_b;
+  char* tmp_c;
+
+  tmp_a = stringconcat(root, tau_char);
+
+  // printf("tmp_a = %s\n", tmp_a);
+
+  tmp_b = stringconcat("_seed", seed_char);
+
+  // printf("tmp_b = %s\n", tmp_b);
+
+  tmp_c = stringconcat(tmp_a, tmp_b);
+
+  // printf("tmp_c = %s\n", tmp_c);
+
+  filename = stringconcat(tmp_c, ".dat");
+
+  // printf("Final name %s\n\n", filename);
+
+  return filename;
+}

make2tar.py

+
+import os, sys
+import re
+import shutil
+import subprocess
+
+makefile='Makefile'
+Array_SubFolders=list()
+
+#======================================================================
+#Function to replace ${key} with os.environ[key]
+#======================================================================
+
+def replace_env_variables(value=None):
+
+    if 'MC_INSTDIR' in value:
+        modified_value=value.replace('${MC_INSTDIR}', os.environ['MC_INSTDIR'])
+    else:
+
+        modified_value=value
+
+    return modified_value
+
+#======================================================================
+def Create_Makefile(FileName=None):
+
+
+    with open(FileName, 'w') as file:
+        # Definizione delle variabili
+        file.write("CC = gcc\n")
+        file.write("CFLAGS = -Wall -Ofast -I${HEADERS_GSL} -I ./include -I ${HEADERS_MPI} -I ${HEADERS_CFITSIO}\n")
+        file.write("LDFLAGS = -L ${LIB_CFITSIO} -L ${LIB_GSL} -L ${LIB_MPI} -lcfitsio -lm -lgsl -lgslcblas -lmpi\n\n\n")
+
+        file.write("SRCDIR =./\n")
+        file.write("BINDIR =./\n")
+        file.write("OBJDIR =./obj\n")
+
+        file.write("# Funzione ricorsiva per cercare file sorgenti in tutte le directory\n")
+        file.write("find_sources = $(wildcard $1/*.c)\n")
+
+        file.write("dirs := $(wildcard $(SRCDIR)*/)\n\n")
+
+        file.write("#Trova tutti i file sorgenti in tutte le directory\n")
+        file.write("SOURCES := $(foreach dir, $(dirs),$(call find_sources, $(dir)))\n\n")
+
+        file.write("# Genera i file oggetto corrispondenti ai file sorgenti\n")
+
+        file.write("OBJECTS = $(patsubst %.c, %.o, $(SOURCES))\n\n")
+
+        file.write("#Nome dell'eseguibile\n")
+        file.write("EXECUTABLE = $(BINDIR)/my_program\n\n\n")
+
+        file.write("# Regola di default: compila l'eseguibile\n")
+        file.write("all: $(EXECUTABLE)\n\n")
+
+        file.write("$(EXECUTABLE): $(OBJECTS)\n")
+        file.write("\t${CC} $^ -o $@ ${LDFLAGS}\n\n")
+
+        file.write("clean:\n")
+        file.write("\trm -rf $(OBJDIR) $(BINDIR) $(OBJECTS)\n\n")
+
+        # Dichiarazione di phony targets
+        file.write(".PHONY: all\n")
+        file.write("\tclean\n\n")
+
+        file.write("%.o: %.c $(DEPS)\n")
+        file.write("\t$(CC) -c -o  $@ $< $(CFLAGS)\n\n")
+
+        # Crea le directory obj e bin se non esistono
+        file.write("$(shell mkdir -p $(OBJDIR) $(BINDIR))\n")
+
+        file.write("debug:\n")
+        file.write("\t@echo \"SOURCES: $(SOURCES)\"\n")
+        file.write("\t@echo \"OBJECTS: $(OBJECTS)\" \n\n")
+
+    file.close()
+
+#===========================================================
+#First check if the enviroment variable is defined
+#===========================================================
+
+if 'MC_INSTDIR' not in os.environ:
+    print('Warning: the enviroment variable MC_INST_DIR is not defined')
+    os.sys.exit()
+
+#===========================================================
+#Define the root folder for storing the subfolders
+#===========================================================
+
+MainFolder='tarball/'
+
+if not os.path.exists(MainFolder):
+    os.makedirs(MainFolder)
+
+#=====================================================================
+#Creates the subfolder with include files and copy to it
+#=====================================================================
+
+LocalHeaders_Folder = MainFolder + "include/"
+try:
+    os.makedirs(LocalHeaders_Folder)
+except:
+    print('Failed to create folder ', LocalHeaders_Folder)
+    #os.sys.exit()
+
+
+MC_Headers_Folder = os.environ['MC_INSTDIR']+'/../../include/'
+
+for filename in os.listdir(MC_Headers_Folder):
+
+    file_path = os.path.join(MC_Headers_Folder, filename)
+    try:
+        shutil.copy2(file_path, LocalHeaders_Folder)
+    except Exception as e:
+        print(f"Error: {e}")
+
+#===========================================================
+# Open the file for reading
+#===========================================================
+
+with open(makefile, 'r') as file:
+    # Read the lines of the file
+    all_lines = file.readlines()
+
+#==============================================================
+#First read the local enviroment variables from the Makefile
+#==============================================================
+
+Enviroment_Variables = dict()
+
+for line in all_lines:
+
+    if "=" in line:
+
+        string_value = line.split('=')
+        string_left = string_value[0]
+        string_right = string_value[1]
+
+        if ".o" not in string_right and len(string_right.split(' '))==2:
+
+            Enviroment_Variables[string_left.strip()]=string_right.strip()
+
+for key, value in Enviroment_Variables.items():
+    if 'PWD' in value:
+        Enviroment_Variables[key] = os.environ['PWD']
+
+Enviroment_Variables['MC_INSTDIR']=os.environ['MC_INSTDIR']
+
+
+for key, value in Enviroment_Variables.items():
+
+    modified_value=replace_env_variables(value=value)
+    Enviroment_Variables[key]=modified_value
+
+
+#=======================================
+#Select lines containing "OBJ=
+#==================================================
+
+pattern = r'\${(.*?)}'
+
+for line in all_lines:
+
+    if "OBJ=" in line:
+
+        string_value = line.split('=')
+        string_left = string_value[0]
+        string_right = string_value[1].strip()
+
+        #print('==================================================')
+
+        for key in Enviroment_Variables.keys():
+
+            if key in string_left:
+
+                #=======================================================
+                #Create the subfolder
+                #=======================================================
+
+                sub_folder = MainFolder + key.lower()
+                Array_SubFolders.append(key.lower())
+
+                if not os.path.exists(sub_folder):
+                    try:
+                        os.makedirs(sub_folder)
+                    except:
+                        print('Failed to create folder ', sub_folder)
+                        os.sys.exit()
+                #=======================================================
+
+                #print('Files corresponding to Enviroment Var ', key)
+                object_files = string_right.split(' ')
+
+                #=====================================================================
+                #Substitute the .o extension with .c
+                #=====================================================================
+
+                #print('Destination folder ', sub_folder)
+                for file_obj in object_files:
+
+                    if len(file_obj) > 0:
+
+                        source_file_fullpath=file_obj.replace(".o", ".c")
+
+                        for key in Enviroment_Variables.keys():
+
+                            if key in source_file_fullpath:
+
+                                source_file_fullpath=source_file_fullpath.replace('${'+key+'}', Enviroment_Variables[key])
+
+                                try:
+                                    # Copy the file to the destination folder
+                                    shutil.copy(source_file_fullpath, sub_folder)
+                                    #print(f"File '{explicit_file_path}' copied to '{sub_folder}' successfully.")
+                                except FileNotFoundError:
+                                    print(f"Error: Source file '{source_file_fullpath}' not found.")
+                                    os.sys.exit()
+                                except PermissionError:
+                                    print("Error: Permission denied. Make sure you have the necessary permissions.")
+                                except Exception as e:
+                                    print(f"An error occurred: {e}")
+
+
+Create_Makefile(FileName=MainFolder+'Makefile')
+
+print('\nThe source tree has been created in subfolder %s' % (MainFolder))
+print('Type command tar cfv my_archive.tar  %s to create the tarball' % (MainFolder))
+
+
+
+
+
+
+
+

make_header_polar.py

+#!/usr/bin/env python3
+import math
+import numpy as np
+import sys
+import os
+
+
+#==============================================================
+
+if len(sys.argv) < 4:
+    print("\nScript to create the polarization_structure.h of the MC ray-tracing code")
+    print("Synopsis: make_header_polar.py emin=emin emax=emax  nbin=nbin\n")
+    os.sys.exit()
+
+
+
+
+for inputval in sys.argv:
+
+    if inputval.startswith('emin') == True:
+        emin = float(inputval[inputval.rfind('=') + 1:])
+
+    if inputval.startswith('emax') == True:
+        emax = float(inputval[inputval.rfind('=') + 1:])
+
+    if inputval.startswith('nbin') == True:
+        nbin = int(inputval[inputval.rfind('=') + 1:])
+
+
+
+
+#==============================================================
+
+array_energy=np.zeros(nbin)
+ymin = math.log(emin / 511.)
+ymax = math.log(emax / 511.)
+
+hy = (ymax - ymin) / (nbin - 1)
+y = ymin
+
+c_array="static const double ene_bound[POLAR_NBOUNDS]={"
+for ii in range(nbin):
+
+    energy=math.exp(y)
+
+    array_energy[ii]=energy*511
+    if ii == nbin-1:
+        final_string="};"
+    else:
+        final_string=", "
+
+    c_array=c_array+str(round(array_energy[ii],3))+final_string
+
+    y = y + hy
+
+#===========================================================
+output="polarization_structure.h"
+f = open(output, "w")
+
+f.write("#ifndef POLARIZATION_STRUCTURE_H_ \n")
+f.write("#define POLARIZATION_STRUCTURE_H_ \n\n")
+
+f.write("#include <functions.h>\n\n")
+
+f.write("#define POLAR_NBOUNDS %d\n\n" % (nbin))
+
+f.write(c_array)
+f.write("\n\n\n\n")
+
+f.write("typedef struct {\n\n")
+
+f.write("double* array_Is;\n")
+f.write("double* array_Qs;\n")
+f.write("double* array_Us;\n")
+
+f.write("double* polar_degree;\n")
+f.write("double* csi_angle;\n")
+
+f.write("uint32_t *counter;\n\n")
+
+
+f.write("} stokes_parameters;\n\n\n")
+
+f.write("#endif\n")
+
+f.close()
+
+#===========================================================
+
+print('Written file', output)
\ No newline at end of file

mc_slab.c

+#include "iteration_functions.h"
+
+#define NSC_MAX 1000
+
+/*==================================================================*/
+/*Function prototypes*/
+/*==================================================================*/
+
+void comptonization(double E_lab,
+                    double* k_lab,
+                    double* Efield_lab,
+                    double* Magfield_lab,
+                    int* PolarFlag,
+                    double* E_new_lab,
+                    double* k_new_lab,
+                    double* Efield_lab_new,
+                    double* Magfield_lab_new);
+
+double distance_bottom(double z0, double kz);
+double distance_top(double z0, double kz, double H);
+
+/*==================================================================*/
+/*Datatypes defined in the header globals.h*/
+/*==================================================================*/
+
+double obsmindeg;
+double obsmaxdeg;
+
+double* col1;
+double* col2;
+double* array_angles;
+double* array_ubounds;
+
+double** array_gb;
+double** ptr_gb;
+double** ptr_2colfile;
+double** array_gb;
+
+int nstepangles;
+
+char* matrix_gb;
+
+
+const gsl_root_fsolver_type* T_phi;
+gsl_root_fsolver* s_phi;
+
+
+const gsl_root_fsolver_type* T_disk;
+gsl_root_fsolver* s_disk;
+
+const gsl_root_fsolver_type* T_maxw;
+gsl_root_fsolver *s_maxw;
+
+
+/*==================================================================*/
+
+int slab_mc(int nphot, int seed)
+{
+  int ii;
+  int jj;
+  int PolarFlag;
+  int FlagLocation;
+  int n_sc;
+  int nph_esc;
+  int photon_diffusion[NSC_MAX];
+  int status;
+  double E_lab;
+  double E_new_lab;
+  double k_lab[3];
+  double k_new_lab[3];
+  double Efield_lab[3];
+  double Efield_lab_new[3];
+  double Magfield_lab[3];
+  double Magfield_lab_new[3];
+  double pos[3];
+  double Pdeg;
+
+  double cos_theta;
+  double sin_theta;
+  double phi;
+  double NeDisk;
+  double r_units;
+  double Hphot;
+  double L_top;
+  double tau;
+  double eps;
+  double l_int;
+  double csi;
+
+  T_disk = gsl_root_fsolver_bisection;
+  s_disk = gsl_root_fsolver_alloc(T_disk);
+
+  T_phi = gsl_root_fsolver_bisection;
+  s_phi = gsl_root_fsolver_alloc(T_phi);
+  
+  /*This is for sampling beta from subrelativistics maxwellian*/
+
+  T_maxw = gsl_root_fsolver_bisection;
+  s_maxw = gsl_root_fsolver_alloc(T_maxw);
+
+
+  nstepangles = 40;
+  obsmindeg = 0.;
+  obsmaxdeg = 90.;
+
+  ptr_gb = read_array(
+      stringconcat(getenv("MC_INSTDIR"), "/electron_distributions/maxwellian/temp_random.txt"));
+
+  matrix_gb =
+      stringconcat(getenv("MC_INSTDIR"), "/electron_distributions/maxwellian/gammabeta_lab.fits");
+  array_gb = read_matrix(matrix_gb);
+
+  double array_ene[NSTEP_ENE];
+  make_array_energy(array_ene, 0.01, 20, NSTEP_ENE);
+
+  double** array_counts;
+  array_counts = dmatrix(0, NSTEP_ENE - 1, 0, nstepangles - 1);
+
+  double** reduce_array_counts;
+  reduce_array_counts = dmatrix(0, NSTEP_ENE - 1, 0, nstepangles - 1);
+
+  for (ii = 0; ii < NSTEP_ENE - 1; ii++)
+  {
+    for (jj = 0; jj < nstepangles; jj++)
+    {
+      array_counts[ii][jj] = 0;
+      reduce_array_counts[ii][jj] = 0;
+    }
+  }
+
+  /*=======================================================================*/
+
+  Hphot = 0.5;
+  r_units = NS_MASS * R_SCW_SUN;
+  NeDisk = disktau / (SIGMA_T * r_units * Hphot);
+
+  /*=======================================================================*/
+
+  if (DiskPolarizationFlag == FROMFILE)
+  {
+    radiation_field_t* radiation_field = (radiation_field_t*)malloc(sizeof(radiation_field_t));
+
+    read_resultsfile(polarinput, radiation_field, MC_SLAB);
+
+    setup_gsl_objects(radiation_field, MC_SLAB);
+  }
+
+  Pdeg = 0;
+  nph_esc = 0;
+  cos_theta = 0;
+
+  for (ii = 0; ii < NSC_MAX; ii++)
+  {
+    photon_diffusion[ii] = 0;
+  }
+  /*===================================================================*/
+
+  array_angles = (double*)malloc((nstepangles) * sizeof(double));
+  array_ubounds = (double*)malloc((nstepangles + 1) * sizeof(double));
+  make_array_angles(array_angles, array_ubounds);
+
+  stokes_parameters struct_stokes[POLAR_NBOUNDS - 1];
+  allocate_stokes_pointer(struct_stokes);
+
+  stokes_parameters struct_stokes_average[POLAR_NBOUNDS - 1];
+  allocate_stokes_pointer(struct_stokes_average);
+
+  /*=============================================================*/
+  /*Here starts loop over photon */
+  /*=============================================================*/
+
+  for (ii = 0; ii < nphot; ii++)
+  {
+    //======================================================================== */
+    // New unpolarized photon */
+    // ===================================================================== */
+
+    if (rank == 0 && ii % 5000 == 0)
+    {
+      printf("======================\n");
+      printf("Photon %d\n", ii);
+    }
+
+    n_sc = 0;
+
+    E_lab = bb_spec(ktbb);
+    E_new_lab = 0;
+
+    pos[0] = 0;
+    pos[1] = 0;
+
+    if (seed == 1)
+    {
+      pos[2] = Hphot / 2;
+      cos_theta = 2 * clock_random() - 1;
+    }
+    else if (seed == 2)
+    {
+      pos[2] = clock_random() * Hphot;
+      cos_theta = 2 * clock_random() - 1;
+    }
+    else
+    {
+      pos[2] = 0;
+
+      if (DiskPolarizationFlag == FROMFILE)
+      {
+        csi = clock_random();
+        status = gsl_spline_eval_e(Spline_sample_Iu, csi, xacc_1d, &cos_theta);
+
+        if (status)
+        {
+          printf("Error\n");
+          printf("File %s\n", __FILE__);
+        }
+
+        printf("sampling %lf\n", cos_theta);
+      }
+      else
+      {
+        cos_theta = clock_random();
+      }
+    }
+
+    sin_theta = sqrt(1 - cos_theta * cos_theta);
+    phi = 2 * PI * clock_random();
+
+    k_lab[0] = sin_theta * cos(phi);
+    k_lab[1] = sin_theta * sin(phi);
+    k_lab[2] = cos_theta;
+
+    if (DiskPolarizationFlag == FROMFILE)
+    {
+      electricfield_fromfile(k_lab, Efield_lab, Magfield_lab, &Pdeg);
+
+      if (Pdeg < 0)
+        Pdeg = -Pdeg;
+
+      if (clock_random() < Pdeg)
+      {
+        PolarFlag = 1;
+      }
+      else
+      {
+        PolarFlag = 0;
+      }
+    }
+
+    else
+    {
+      unpolarized_seedphot(k_lab, Efield_lab, Magfield_lab);
+      PolarFlag = 0;
+    }
+
+    /*====================================================================*/
+
+    FlagLocation = 1;
+
+    while (FlagLocation)
+    {
+      if (k_lab[2] >= 0)
+      {
+        L_top = distance_top(pos[2], k_lab[2], Hphot);
+      }
+      else
+      {
+        L_top = distance_bottom(pos[2], k_lab[2]);
+      }
+
+      /*====================================================================*/
+      /*Compute the optical depth from top of bottom of the atmosphere*/
+      /*=====================================================================*/
+
+      tau = NeDisk * SIGMA_T * r_units * L_top;
+      eps = -log(clock_random());
+
+      if (tau > eps)
+      {
+        l_int = eps / (NeDisk * SIGMA_T * r_units);
+        // printf("Photon %d will interact at dist %5.4f\n", ii, l_int);
+
+        pos[0] = pos[0] + k_lab[0] * l_int;
+        pos[1] = pos[1] + k_lab[1] * l_int;
+        pos[2] = pos[2] + k_lab[2] * l_int;
+
+        comptonization(E_lab, k_lab, Efield_lab, Magfield_lab, &PolarFlag, &E_new_lab, k_new_lab,
+                       Efield_lab_new, Magfield_lab_new);
+
+        n_sc++;
+
+        k_lab[0] = k_new_lab[0];
+        k_lab[1] = k_new_lab[1];
+        k_lab[2] = k_new_lab[2];
+
+        E_lab = E_new_lab;
+
+        Efield_lab[0] = Efield_lab_new[0];
+        Efield_lab[1] = Efield_lab_new[1];
+        Efield_lab[2] = Efield_lab_new[2];
+
+        Magfield_lab[0] = Magfield_lab_new[0];
+        Magfield_lab[1] = Magfield_lab_new[1];
+        Magfield_lab[2] = Magfield_lab_new[2];
+      }
+      else
+      {
+        if (k_lab[2] > 0)
+        {
+          // printf("Photon %d escaping after %d scattering\n", ii, n_sc);
+
+          compute_stokes(E_lab, k_lab, Efield_lab, struct_stokes, array_ubounds, 0);
+          photon_diffusion[n_sc]++;
+          collect_photons(E_new_lab, k_lab, array_ene, array_counts, 0);
+
+          nph_esc++;
+        }
+        else
+        {
+          // printf("Photon %d escaping\n", ii);
+        }
+
+        FlagLocation = 0;
+      }
+
+    } // end of while (FlagLocation)
+
+    // newphot:;
+  }
+
+  /*==================================================================*/
+
+  for (ii = 0; ii < POLAR_NBOUNDS - 1; ii++)
+  {
+    MPI_Reduce(struct_stokes[ii].array_Is, struct_stokes_average[ii].array_Is, nstepangles,
+               MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+
+    MPI_Reduce(struct_stokes[ii].array_Qs, struct_stokes_average[ii].array_Qs, nstepangles,
+               MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+
+    MPI_Reduce(struct_stokes[ii].array_Us, struct_stokes_average[ii].array_Us, nstepangles,
+               MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+
+    MPI_Reduce(struct_stokes[ii].counter, struct_stokes_average[ii].counter, nstepangles, MPI_INT,
+               MPI_SUM, 0, MPI_COMM_WORLD);
+  }
+
+  for (ii = 0; ii < NSTEP_ENE; ii++)
+  {
+    for (jj = 0; jj < nstepangles; jj++)
+    {
+      MPI_Reduce(&array_counts[ii][jj], &reduce_array_counts[ii][jj], 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+    }
+  }
+
+  /*==================================================================*/
+
+  if (rank == 0)
+  {
+    compute_stokes(1, k_lab, Efield_lab, struct_stokes_average, array_ubounds, 1);
+
+    printf("Percentage of escaping photons %4.3f\n", 1. * nph_esc / nphot);
+
+    dat_diffusion = fopen(diffusion, "w");
+    fprintf(dat_diffusion, "!Percentage of escaping photons %4.3f\n", 1. * nph_esc / nphot);
+
+    for (ii = 0; ii < 300; ii++)
+    {
+      fprintf(dat_diffusion, " %d  %d\n", ii, photon_diffusion[ii]);
+    }
+
+    fclose(dat_diffusion);
+  }
+
+  collect_photons(0, k_lab, array_ene, reduce_array_counts, 1);
+
+  dat_energy_polar = fopen(polarfile, "w");
+
+  write_polardata(struct_stokes_average, ene_bound, dat_energy_polar);
+
+  fclose(dat_energy_polar);
+
+  gsl_root_fsolver_free(s_phi);
+  gsl_root_fsolver_free(s_disk);
+
+  return 0;
+
+} // End of function
+
+void comptonization(double E_lab,
+                    double* k_lab,
+                    double* Efield_lab,
+                    double* Magfield_lab,
+                    int* PolarFlag,
+                    double* E_new_lab,
+                    double* k_new_lab,
+                    double* Efield_lab_new,
+                    double* Magfield_lab_new)
+
+{
+  int jj;
+  int status;
+  double phik;
+  double E_erf;
+  double E_new_erf;
+  double csi;
+  double costheta_erf;
+  double sintheta_erf;
+  double beta_el;
+
+  double k_erf[3];
+  double k_new_erf[3];
+  double el_vel_vect[3];
+
+  double ElectField_new_erf[3];
+  double ElectField_old_erf[3];
+
+  double MagField_Erf[3];
+
+  //================================================================== */
+  // Find the K' components in the ERF */
+  // If also polarization is computed, transform also the electric */
+  // field vector to ERF */
+  //================================================================== */
+
+  beta_el = sample_electron(E_lab, k_lab, kte, ACCRETION_DISK, el_vel_vect);
+
+  versor_boost(k_lab, el_vel_vect, k_erf, 1);
+
+  double beta_new=subrel_maxwellian (kte);
+
+  //printf("pluto %lf %lf\n", beta_el, beta_new);
+
+  elmag_lorentztransform(Efield_lab, Magfield_lab, el_vel_vect, ElectField_old_erf, MagField_Erf, 1);
+
+  E_erf = energy_lorentz_boost(E_lab, el_vel_vect, k_lab, 1);
+
+  // printf("beta %lf %lf \n", beta_el, E_erf);
+
+  csi = clock_random();
+
+  costheta_erf = (-1 + pow(-2 + 4 * csi + sqrt(5 - 16 * csi + 16 * pow(csi, 2)), 2 / 3.)) /
+                 pow(-2 + 4 * csi + sqrt(5 - 16 * csi + 16 * pow(csi, 2)), 1 / 3.);
+
+  sintheta_erf = sqrt(1 - costheta_erf * costheta_erf);
+
+  // ==================================================================
+  // Find the new energy photon in the erf after scattering
+  //==================================================================
+
+  E_new_erf = new_energy_erf(E_erf, costheta_erf);
+
+  if (*PolarFlag == 1)
+  {
+    phik = sample_azim_angle(E_erf, E_new_erf, costheta_erf);
+  }
+  else
+  {
+    phik = 2 * PI * clock_random();
+  }
+
+  //================================================================== */
+  /*Find the components of the photon versor after scattering*/
+  /*Use Equation(2) from Peplow, 1999*/
+  /*Nuclear Science and Engineering, 131, 132*/
+  //================================================================== */
+
+  for (jj = 0; jj < 3; jj++)
+  {
+    k_new_erf[jj] = k_erf[jj] * costheta_erf + ElectField_old_erf[jj] * cos(phik) * sintheta_erf +
+                    MagField_Erf[jj] * sin(phik) * sintheta_erf;
+  }
+
+  //================================================================== */
+  // Check the probability for depolarization*/
+  // and compute the new electric field vector
+  //================================================================== */
+
+  status = depolarization_probability(E_erf, E_new_erf, costheta_erf, phik, PolarFlag);
+
+  if (status == 1)
+  {
+    exit(1);
+  }
+
+  electfield_as_new(k_new_erf, ElectField_old_erf, ElectField_new_erf, *PolarFlag);
+
+  vect_prod(k_new_erf, ElectField_new_erf, MagField_Erf);
+
+  // =========================================================================== */
+  // Back convert from (K1)' in ERF to lab*/
+  // =========================================================================== */
+
+  versor_boost(k_new_erf, el_vel_vect, k_new_lab, 2);
+
+  // ================================================================== */
+  // Back-transform electric field versor to the lab frame */
+  // ================================================================== */
+
+  elmag_lorentztransform(ElectField_new_erf, MagField_Erf, el_vel_vect, Efield_lab_new, Magfield_lab_new, 2);
+
+  // =========================================================================== */
+  // Back convert from (K1)' in ERF to lab*/
+  // =========================================================================== */
+
+  versor_boost(k_new_erf, el_vel_vect, k_new_lab, 2);
+
+  // ================================================================== */
+  // Back convert electron energy from ERF to lab */
+  // ================================================================== */
+
+  *E_new_lab = energy_lorentz_boost(E_new_erf, el_vel_vect, k_new_erf, 2);
+
+  if (isnan(*E_new_lab) || *E_new_lab < 0)
+  {
+    printf("Error in photon energy in lab frame after scattering: value %5.4f\n", *E_new_lab);
+    exit(1);
+  }
+}
+
+double distance_top(double z_pos, double kz, double H)
+{
+  double t;
+
+  t = (H - z_pos) / kz;
+
+  return t;
+}
+
+double distance_bottom(double z_pos, double kz)
+{
+  double t;
+
+  t = -z_pos / kz;
+
+  return t;
+}

ode_example.c

+int main(void)
+{
+  double mu = 10;
+  gsl_odeiv2_system sys = {func, jac, 2, &mu};
+
+  gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, 1e-3, 1e-8, 1e-8);
+
+  double t = 0.0;
+  double y[2] = {1.0, 0.0};
+  int i, s;
+
+  for (i = 0; i < 100; i++)
+  {
+    s = gsl_odeiv2_driver_apply_fixed_step(d, &t, 1e-3, 1000, y);
+
+    if (s != GSL_SUCCESS)
+    {
+      printf("error: driver returned %d\n", s);
+      break;
+    }
+
+    printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
+  }
+
+  gsl_odeiv2_driver_free(d);
+  return s;
+}

read_iterationresults.c

+#include "iteration_functions.h"
+
+void read_iterationresults(char* filename, radiation_field_t* radiation_field)
+{
+  FILE* txt;
+  char line[NDIM];
+  int32_t tot_lines, count = 0;
+
+  /* int ii; */
+
+  txt = fopen(filename, "r");
+
+  if (txt == NULL)
+  {
+    printf("\nWARNING\n");
+    printf("Routine read_phdensity in file %s\n ", __FILE__);
+    printf("File %s not found\n\n", filename);
+    exit(1);
+  }
+
+  tot_lines = 0;
+
+  while (!feof(txt))
+  {
+    if (fgets(line, NDIM, txt) == NULL)
+    {
+      break;
+    }
+    else if (strstr(line, "#") != NULL)
+    {
+      continue;
+    }
+    else if (strstr(line, "!") != NULL)
+    {
+      continue;
+    }
+    else
+    {
+      tot_lines++;
+    }
+  }
+
+  rewind(txt);
+
+  /*===================================================================*/
+  /*Set pointer to structure defined in header file functions.h*/
+  /*===================================================================*/
+
+  radiation_field->u = (double*)malloc(tot_lines * sizeof(double));
+  radiation_field->pdeg = (double*)malloc(tot_lines * sizeof(double));
+  radiation_field->limbdark = (double*)malloc(tot_lines * sizeof(double));
+
+  count = 0;
+
+  while (!feof(txt))
+  {
+    if (fgets(line, NDIM, txt) == NULL)
+
+      break;
+
+    else if (strstr(line, "#") != NULL)
+    {
+      continue;
+    }
+    else if (strstr(line, "!") != NULL)
+    {
+      continue;
+    }
+
+    else
+    {
+      sscanf(line, "%lg %lg %lg\n", &radiation_field->u[count], &radiation_field->pdeg[count],
+             &radiation_field->limbdark[count]);
+
+      printf("Angle %lf  Pdeg %lf I(u) %5.4e\n", radiation_field->u[count],
+             radiation_field->pdeg[count], radiation_field->limbdark[count]);
+
+      count++;
+    }
+  }
+
+  radiation_field->nlines = count;
+}

round_to_digits.c

+#include "iteration_functions.h"
+
+double round_to_digits(double value, int digits)
+{
+  if (value == 0.0) // otherwise it will return 'nan' due to the log10() of zero
+    return 0.0;
+
+  double factor = pow(10.0, digits - ceil(log10(fabs(value))));
+  return round(value * factor) / factor;
+}

slab_polarization.c

+#include "iteration_functions.h"
+
+gsl_spline2d* Spline_I2d_l;
+gsl_spline2d* Spline_I2d_r;
+gsl_interp_accel* xacc_2d;
+gsl_interp_accel* yacc_2d;
+
+gsl_spline* Spline_sample_Iu;
+gsl_spline* Spline_eval_limbdark;
+gsl_spline* Spline_eval_pdeg;
+
+gsl_spline* Spline_tau_Ak;
+gsl_spline* Spline_tau_Bk;
+gsl_spline* Spline_tau_Ck;
+
+gsl_interp_accel* xacc_1d;
+gsl_interp_accel* yacc_1d;
+
+int Nstep_tau;
+int Nstep_mu;
+double eps;
+double u_min;
+double u_max;
+double tau_s;
+
+/*========================================================================*/
+
+void slab_polarization(int kiter, int seed_location)
+{
+  Nstep_mu = 20;
+  Nstep_tau = 100;
+
+  tau0 = disktau / 2.;
+
+  int ii, jj, kk;
+
+  double* array_mu;
+  double* array_tau;
+  double* array_weights;
+  double* Ak_array;
+  double* Bk_array;
+  double* Ck_array;
+
+  double error;
+
+  printf("tau0=%4.3f\n", tau0);
+
+  array_mu = malloc(Nstep_mu * sizeof(double));
+  array_weights = malloc(Nstep_mu * sizeof(double));
+  array_tau = malloc(Nstep_tau * sizeof(double));
+
+  Ak_array = malloc(Nstep_tau * sizeof(double));
+  Bk_array = malloc(Nstep_tau * sizeof(double));
+  Ck_array = malloc(Nstep_tau * sizeof(double));
+
+  eps = 1e-5;
+
+  tau_s = 2 * tau0 - tau_c;
+
+  /*=================================================================*/
+
+  float* roots_u = malloc((Nstep_mu + 1) * sizeof(float));
+  float* weights = malloc((Nstep_mu + 1) * sizeof(float));
+
+  gauleg(0, 1, roots_u, weights, Nstep_mu);
+
+  for (ii = 1; ii <= Nstep_mu; ii++)
+  {
+    array_mu[ii - 1] = roots_u[ii];
+    array_weights[ii - 1] = weights[ii];
+  }
+
+  u_min = array_mu[0];
+  u_max = array_mu[Nstep_mu - 1];
+
+  /*make_linarray(u_min, 1, Nstep_mu, array_mu);*/
+
+  /* exit(1);*/
+  /*=================================================================*/
+
+  make_linarray(0, 2. * tau0, Nstep_tau, array_tau);
+
+  Ikl_intensity* ptr_Ikl;
+  Ikr_intensity* ptr_Ikr;
+
+  ptr_Ikl = malloc(kiter * sizeof(Ikl_intensity));
+  ptr_Ikr = malloc(kiter * sizeof(Ikr_intensity));
+
+  for (ii = 0; ii < kiter; ii++)
+  {
+    ptr_Ikl[ii].matrix = dmatrix(0, Nstep_tau - 1, 0, Nstep_mu - 1);
+    ptr_Ikr[ii].matrix = dmatrix(0, Nstep_tau - 1, 0, Nstep_mu - 1);
+  }
+
+  /*=================================================================*/
+  /*Set the function for photons undegoing no scattering*/
+  /*=================================================================*/
+
+  for (ii = 0; ii < Nstep_tau; ii++)
+  {
+    for (jj = 0; jj < Nstep_mu; jj++)
+    {
+      if (seed_location == SEED_CENTER)
+      {
+        ptr_Ikl[0].matrix[ii][jj] = I0_center(tau0, array_tau[ii], array_mu[jj]);
+        ptr_Ikr[0].matrix[ii][jj] = I0_center(tau0, array_tau[ii], array_mu[jj]);
+      }
+      else if (seed_location == SEED_UNIFORM)
+      {
+        ptr_Ikl[0].matrix[ii][jj] = I0_uniform(tau0, array_tau[ii], array_mu[jj]);
+        ptr_Ikr[0].matrix[ii][jj] = I0_uniform(tau0, array_tau[ii], array_mu[jj]);
+      }
+      else if (seed_location == SEED_DELTA)
+      {
+        ptr_Ikl[0].matrix[ii][jj] = I0_delta(array_tau[ii], tau_s, array_mu[jj]);
+        ptr_Ikr[0].matrix[ii][jj] = I0_delta(array_tau[ii], tau_s, array_mu[jj]);
+        printf("Do not use this configuration for seed photons fot ST85 algorithm\n");
+        exit(1);
+      }
+      else if (seed_location == SEED_BASE)
+      {
+        ptr_Ikl[0].matrix[ii][jj] = I0_base(array_tau[ii], tau_s, array_mu[jj], U_POSITIVE);
+        ptr_Ikr[0].matrix[ii][jj] = I0_base(array_tau[ii], tau_s, array_mu[jj], U_POSITIVE);
+        printf("Do not use this configuration for seed photons fot ST85 algorithm\n");
+        exit(1);
+      }
+      else if (seed_location == SEED_EXP)
+      {
+        ptr_Ikl[0].matrix[ii][jj] = I0_exp(array_tau[ii], tau0, array_mu[jj]);
+        ptr_Ikr[0].matrix[ii][jj] = I0_exp(array_tau[ii], tau0, array_mu[jj]);
+        printf("Do not use this configuration for seed photons fot ST85 algorithm\n");
+        exit(1);
+      }
+
+      else
+      {
+        printf("Wait, do not use this!\n");
+        exit(1);
+      }
+    }
+  }
+
+  /*=================================================================*/
+  /*Setup for 2D interpolation*/
+  /*=================================================================*/
+
+  double* za = malloc(Nstep_tau * Nstep_mu * sizeof(double));
+  const gsl_interp2d_type* T_interp2d = gsl_interp2d_bilinear;
+
+  Spline_I2d_l = gsl_spline2d_alloc(T_interp2d, Nstep_tau, Nstep_mu);
+  Spline_I2d_r = gsl_spline2d_alloc(T_interp2d, Nstep_tau, Nstep_mu);
+
+  xacc_2d = gsl_interp_accel_alloc();
+  yacc_2d = gsl_interp_accel_alloc();
+
+  /*=================================================================*/
+  /*Initialize objects for 1D interpolation*/
+  /*=================================================================*/
+
+  Spline_tau_Ak = gsl_spline_alloc(gsl_interp_cspline, Nstep_tau);
+  Spline_tau_Bk = gsl_spline_alloc(gsl_interp_cspline, Nstep_tau);
+  Spline_tau_Ck = gsl_spline_alloc(gsl_interp_cspline, Nstep_tau);
+
+  xacc_1d = gsl_interp_accel_alloc();
+
+  /*=================================================================*/
+
+  AngularFunctionParams* ptr_data = malloc(sizeof(AngularFunctionParams));
+  ptr_data->seed_location = seed_location;
+  ptr_data->tau_s = tau_s;
+
+  /*=============================================================*/
+  /*Il and Ir functions*/
+  /*=============================================================*/
+
+  gsl_function I_l;
+  I_l.function = &Il_integral;
+
+  gsl_function I_r;
+  I_r.function = &Ir_integral;
+
+  /*=============================================================*/
+
+  gsl_integration_workspace* w = gsl_integration_workspace_alloc(1000);
+
+  double params_Ifunct[2];
+  double result_Ir;
+  double result_Il;
+
+  /*=============================================================*/
+  /*Here start the iteration over scattering*/
+  /*=============================================================*/
+
+  printf("Running iteration procedure using ST85 algorithm...\n");
+
+  for (kk = 1; kk < kiter; kk++)
+  {
+    if (kk % 5 == 0)
+      printf("iteration k=%d\n", kk);
+    /*=============================================*/
+    /*Loop over tau*/
+    /*=============================================*/
+
+    for (ii = 0; ii < Nstep_tau; ii++)
+    {
+      ptr_data->tau = array_tau[ii];
+      ptr_data->k = kk;
+
+      Compute_Ak_array(ptr_data, ptr_Ikl[kk - 1].matrix, za, array_tau, array_mu, kk, w, Ak_array);
+
+      gsl_spline_init(Spline_tau_Ak, array_tau, Ak_array, Nstep_tau);
+
+      Compute_Bk_array(ptr_data, ptr_Ikl[kk - 1].matrix, za, array_tau, array_mu, kk, w, Bk_array);
+
+      gsl_spline_init(Spline_tau_Bk, array_tau, Bk_array, Nstep_tau);
+
+      Compute_Ck_array(ptr_data, ptr_Ikr[kk - 1].matrix, za, array_tau, array_mu, kk, w, Ck_array);
+
+      gsl_spline_init(Spline_tau_Ck, array_tau, Ck_array, Nstep_tau);
+
+      /*=============================================*/
+      /*Loop over mu*/
+      /*=============================================*/
+
+      for (jj = 0; jj < Nstep_mu; jj++)
+      {
+        params_Ifunct[0] = array_tau[ii];
+        params_Ifunct[1] = array_mu[jj];
+
+        I_l.params = params_Ifunct;
+        I_r.params = params_Ifunct;
+
+        if (array_tau[ii] == 2 * tau0)
+        {
+          ptr_Ikl[kk].matrix[ii][jj] = 0;
+          ptr_Ikr[kk].matrix[ii][jj] = 0;
+        }
+        else
+        {
+          gsl_integration_qags(&I_l, array_tau[ii], 2 * tau0, eps, eps, 1000, w, &result_Il, &error);
+          gsl_integration_qags(&I_r, array_tau[ii], 2 * tau0, eps, eps, 1000, w, &result_Ir, &error);
+
+          if (isnan(result_Il))
+          {
+            printf("Warning: Il is nan for tau=%lf mu=%lf\n", array_tau[ii], array_mu[jj]);
+            result_Il = 0;
+          }
+
+          if (isnan(result_Ir))
+          {
+            printf("Warning: Ir is nan for tau=%lf mu=%lf\n", array_tau[ii], array_mu[jj]);
+            result_Ir = 0;
+          }
+
+          ptr_Ikl[kk].matrix[ii][jj] = result_Il;
+          ptr_Ikr[kk].matrix[ii][jj] = result_Ir;
+        }
+      }
+
+    } /*Loop over tau*/
+
+  } /*Loop over kiter*/
+
+  /*========================================================================*/
+  /*Now sum results*/
+  /*========================================================================*/
+
+  Iklr_intensity* ptr_void = 0;
+
+  compute_results(2, kiter, Nstep_tau, Nstep_mu, array_mu, array_weights, ptr_void, ptr_Ikl, ptr_Ikr);
+
+} // End of function
+
+/*==========================================================================*/
+void make_linarray(double a, double b, int nstep, double* my_array)
+{
+  int ii;
+  double hstep;
+
+  hstep = (b - a) / (nstep - 1);
+
+  for (ii = 0; ii < nstep; ii++)
+  {
+    my_array[ii] = a + hstep * ii;
+    /*printf("ii=%d value=%lf\n", ii, my_array[ii]);*/
+  }
+}
+
+/*========================================================================*/
+/*Distribution of photons with zero scattering*/
+/*========================================================================*/
+
+double theta_funct(double tau, double tau_c)
+{
+  double x;
+  double flag;
+  x = tau_c - tau;
+
+  if (x > 0)
+  {
+    flag = 1;
+  }
+  else
+  {
+    flag = 0;
+  }
+
+  // printf("Argument of theta : tau %lf tau_delta %lf flag= %d \n", tau, tau_c, flag);
+  return flag;
+}
+
+double I0_center(double tau0, double tau, double mu)
+{
+  double arg_exp, value;
+
+  arg_exp = (tau0 - tau) / mu;
+
+  value = 1 / mu * theta_funct(tau, tau0) * exp(-arg_exp);
+
+  return value;
+}
+
+double I0_delta(double tau, double tau_s, double mu)
+
+{
+  double value;
+
+  value = 1 / mu * theta_funct(tau, tau_s) * exp(-(tau_s - tau) / mu);
+
+  return value;
+}
+
+double I0_uniform(double tau0, double tau, double mu)
+{
+  double value;
+
+  value = 1 - exp(-(2 * tau0 - tau) / mu);
+
+  return value;
+}
+
+double I0_exp(double tau, double tau0, double mu)
+{
+  double value;
+  double tau_fold = 0.2;
+
+  value = tau_fold / (tau_fold - mu) * (exp(-(2 * tau0 - tau) / tau_fold) - exp(-(2 * tau0 - tau) / mu));
+
+  return value;
+}
+
+double I0_base(double tau, double tau_s, double mu, int sign_mu)
+{
+  double value;
+
+  // printf("tau_s=%4.3f tau_c %4.3f\n", tau_s, tau_c);
+
+  if (sign_mu > 0)
+  {
+    if (tau > tau_s)
+    {
+      value = 1 - exp(-(2 * tau0 - tau) / mu);
+    }
+    else
+    {
+      value = (-exp(-(2 * tau0 - tau) / mu) + exp(-(tau_s - tau) / mu));
+    }
+  }
+  else
+  {
+    if (tau < tau_c)
+    {
+      value = 1 - exp(-(tau_c - tau) / mu);
+    }
+    else
+    {
+      value = 0;
+    }
+  }
+
+  return value;
+}
+
+/*========================================================================*/
+/*Functions A^{k-1}, B^{k-1}, C^{k-1}*/
+/*========================================================================*/
+
+void Compute_Ak_array(AngularFunctionParams* ptr_data,
+                      double** Il_funct,
+                      double* za,
+                      double* array_tau,
+                      double* array_mu,
+                      int k,
+                      gsl_integration_workspace* w,
+                      double* Ak_array)
+{
+  int ii, tt, qq;
+  double result, error;
+
+  gsl_function F_Ak;
+  F_Ak.function = &Ak_integral;
+
+  /*==================================*/
+  /*Loop over tau*/
+  /*==================================*/
+
+  for (ii = 0; ii < Nstep_tau; ii++)
+  {
+    for (tt = 0; tt < Nstep_tau; tt++)
+    {
+      for (qq = 0; qq < Nstep_mu; qq++)
+      {
+        gsl_spline2d_set(Spline_I2d_l, za, tt, qq, Il_funct[tt][qq]);
+      }
+    }
+
+    gsl_spline2d_init(Spline_I2d_l, array_tau, array_mu, za, Nstep_tau, Nstep_mu);
+
+    ptr_data->tau = array_tau[ii];
+    ptr_data->k = k;
+
+    F_Ak.params = ptr_data;
+
+    // printf("Integration between %lf - %lf\n", u_min, u_max);
+
+    gsl_integration_qags(&F_Ak, u_min, u_max, eps, eps, 1000, w, &result, &error);
+    Ak_array[ii] = result;
+  }
+}
+
+/*=================================================================================*/
+
+void Compute_Bk_array(AngularFunctionParams* ptr_data,
+                      double** Il_funct,
+                      double* za,
+                      double* array_tau,
+                      double* array_mu,
+                      int k,
+                      gsl_integration_workspace* w,
+                      double* Bk_array)
+{
+  int ii, tt, qq;
+  double result, error;
+
+  gsl_function F_Bk;
+  F_Bk.function = &Bk_integral;
+
+  /*==================================*/
+  /*Loop over tau*/
+  /*==================================*/
+
+  for (ii = 0; ii < Nstep_tau; ii++)
+  {
+    for (tt = 0; tt < Nstep_tau; tt++)
+    {
+      for (qq = 0; qq < Nstep_mu; qq++)
+      {
+        gsl_spline2d_set(Spline_I2d_l, za, tt, qq, Il_funct[tt][qq]);
+      }
+    }
+
+    gsl_spline2d_init(Spline_I2d_l, array_tau, array_mu, za, Nstep_tau, Nstep_mu);
+
+    ptr_data->tau = array_tau[ii];
+    ptr_data->k = k;
+    F_Bk.params = ptr_data;
+
+    gsl_integration_qags(&F_Bk, u_min, u_max, eps, eps, 1000, w, &result, &error);
+    Bk_array[ii] = result;
+  }
+}
+
+/*=========================================================================*/
+
+void Compute_Ck_array(AngularFunctionParams* ptr_data,
+                      double** Ir_funct,
+                      double* za,
+                      double* array_tau,
+                      double* array_mu,
+                      int k,
+                      gsl_integration_workspace* w,
+                      double* Ck_array)
+{
+  int ii, tt, qq;
+  double result, error;
+
+  gsl_function F_Ck;
+  F_Ck.function = &Ck_integral;
+
+  /*==================================*/
+  /*Loop over tau*/
+  /*==================================*/
+
+  for (ii = 0; ii < Nstep_tau; ii++)
+  {
+    for (tt = 0; tt < Nstep_tau; tt++)
+    {
+      for (qq = 0; qq < Nstep_mu; qq++)
+      {
+        gsl_spline2d_set(Spline_I2d_r, za, tt, qq, Ir_funct[tt][qq]);
+      }
+    }
+
+    gsl_spline2d_init(Spline_I2d_r, array_tau, array_mu, za, Nstep_tau, Nstep_mu);
+
+    ptr_data->tau = array_tau[ii];
+    ptr_data->k = k;
+    F_Ck.params = ptr_data;
+
+    gsl_integration_qags(&F_Ck, u_min, u_max, eps, eps, 1000, w, &result, &error);
+    Ck_array[ii] = result;
+  }
+}

solve_rte.c

+#include "iteration_functions.h"
+
+double I0_center_updown(double tau, double u);
+double I0_upward(double tau, double u);
+double I0_uniform_updown(double tau, double u);
+
+void set_spline2d_obj(gsl_spline2d* Spline_I2d, double* za, double** I_funct);
+int set_array_idx(double tau_prev, double tau, double* array, int nstep);
+
+double tau0;
+
+/*================================================================================*/
+
+void solve_rte(int k_iter, int seed_distribution)
+{
+  int ii, jj, kk, status;
+
+  Nstep_mu = 30;
+  Nstep_tau = 100;
+
+  tau0 = disktau / 2.;
+
+  // printf("tau0=%4.3f\n", tau0);
+
+  static double y[4];
+  void* jac = NULL;
+  double step_tau;
+  double* array_tau;
+  double* array_mu;
+  double* array_weights;
+
+  gsl_spline2d* Spline_I2d_l_upstream;
+  gsl_spline2d* Spline_I2d_l_downstream;
+  gsl_spline2d* Spline_I2d_r_upstream;
+  gsl_spline2d* Spline_I2d_r_downstream;
+
+  /*==============================================================================================*/
+  float* roots = malloc((Nstep_mu + 1) * sizeof(float));
+  float* weights = malloc((Nstep_mu + 1) * sizeof(float));
+
+  array_mu = malloc(Nstep_mu * sizeof(double));
+  array_weights = malloc(Nstep_mu * sizeof(double));
+
+  gauleg(0, 1, roots, weights, Nstep_mu);
+
+  for (ii = 1; ii <= Nstep_mu; ii++)
+  {
+    array_mu[ii - 1] = roots[ii];
+    array_weights[ii - 1] = weights[ii];
+
+    // printf("%lf 1 \n", array_mu[ii - 1]);
+  }
+
+  array_tau = malloc(Nstep_tau * sizeof(double));
+  make_linarray(0, 2. * tau0, Nstep_tau, array_tau);
+
+  Iklr_intensity* ptr_Iklr = malloc(k_iter * sizeof(Iklr_intensity));
+
+  for (ii = 0; ii < k_iter; ii++)
+  {
+    ptr_Iklr[ii].Il_matrix_upstream = dmatrix(0, Nstep_tau - 1, 0, Nstep_mu - 1);
+    ptr_Iklr[ii].Ir_matrix_upstream = dmatrix(0, Nstep_tau - 1, 0, Nstep_mu - 1);
+
+    ptr_Iklr[ii].Il_matrix_downstream = dmatrix(0, Nstep_tau - 1, 0, Nstep_mu - 1);
+    ptr_Iklr[ii].Ir_matrix_downstream = dmatrix(0, Nstep_tau - 1, 0, Nstep_mu - 1);
+  }
+
+  /*==============================================================================*/
+
+  for (ii = 0; ii < Nstep_tau; ii++)
+  {
+    for (jj = 0; jj < Nstep_mu; jj++)
+    {
+      if (seed_distribution == SEED_CENTER)
+      {
+        ptr_Iklr[0].Il_matrix_upstream[ii][jj] = I0_center_updown(array_tau[ii], array_mu[jj]);
+        ptr_Iklr[0].Ir_matrix_upstream[ii][jj] = I0_center_updown(array_tau[ii], array_mu[jj]);
+
+        ptr_Iklr[0].Il_matrix_downstream[ii][jj] = I0_center_updown(array_tau[ii], array_mu[jj]);
+        ptr_Iklr[0].Ir_matrix_downstream[ii][jj] = I0_center_updown(array_tau[ii], array_mu[jj]);
+
+        /*if (jj==10)
+        {
+            printf("tau=%5.4f I=%5.4e\n", array_tau[ii], ptr_Iklr[0].Il_matrix_upstream[ii][jj]);
+        }*/
+      }
+      else if (seed_distribution == SEED_UNIFORM)
+      {
+        ptr_Iklr[0].Il_matrix_upstream[ii][jj] = I0_uniform_updown(array_tau[ii], array_mu[jj]);
+        ptr_Iklr[0].Ir_matrix_upstream[ii][jj] = I0_uniform_updown(array_tau[ii], array_mu[jj]);
+
+        ptr_Iklr[0].Il_matrix_downstream[ii][jj] = I0_uniform_updown(array_tau[ii], array_mu[jj]);
+        ptr_Iklr[0].Ir_matrix_downstream[ii][jj] = I0_uniform_updown(array_tau[ii], array_mu[jj]);
+      }
+
+      else if (seed_distribution == SEED_BASE)
+      {
+        ptr_Iklr[0].Il_matrix_upstream[ii][jj] = I0_upward(array_tau[ii], array_mu[jj]);
+        ptr_Iklr[0].Ir_matrix_upstream[ii][jj] = I0_upward(array_tau[ii], array_mu[jj]);
+
+        ptr_Iklr[0].Il_matrix_downstream[ii][jj] = 0;
+        ptr_Iklr[0].Ir_matrix_downstream[ii][jj] = 0;
+      }
+      else
+      {
+        printf(
+            "Error: please select seed=4 (photons at the base) or seed=2 (uniform distribution)\n");
+        exit(1);
+      }
+    }
+  }
+  /*=================================================================*/
+  /*Setup for 2D interpolation*/
+  /*=================================================================*/
+
+  double* Il_upstream_za = malloc(Nstep_tau * Nstep_mu * sizeof(double));
+  double* Il_downstream_za = malloc(Nstep_tau * Nstep_mu * sizeof(double));
+  double* Ir_upstream_za = malloc(Nstep_tau * Nstep_mu * sizeof(double));
+  double* Ir_downstream_za = malloc(Nstep_tau * Nstep_mu * sizeof(double));
+
+  const gsl_interp2d_type* T_interp2d = gsl_interp2d_bilinear;
+
+  Spline_I2d_l_upstream = gsl_spline2d_alloc(T_interp2d, Nstep_tau, Nstep_mu);
+  Spline_I2d_l_downstream = gsl_spline2d_alloc(T_interp2d, Nstep_tau, Nstep_mu);
+
+  Spline_I2d_r_upstream = gsl_spline2d_alloc(T_interp2d, Nstep_tau, Nstep_mu);
+  Spline_I2d_r_downstream = gsl_spline2d_alloc(T_interp2d, Nstep_tau, Nstep_mu);
+
+  xacc_2d = gsl_interp_accel_alloc();
+  yacc_2d = gsl_interp_accel_alloc();
+
+  step_tau = (2 * tau0) / (Nstep_tau - 1);
+
+  // printf("Step tau %5.4f\n", step_tau);
+  /*double check_val=0;*/
+
+  Iklr_intensity* ptr_data = malloc(sizeof(Iklr_intensity));
+  Iklr_intensity* params;
+
+  // const gsl_odeiv_step_type* T = gsl_odeiv_step_rk4;
+
+  /*================================================================================*/
+  /*Start the loop over k-iterations*/
+  /*================================================================================*/
+
+  double tau, tau_prev;
+  // step_tau = 0.01;
+
+  for (kk = 1; kk < k_iter; kk++)
+  {
+    printf("Iteration k=%d....\n\n", kk);
+
+    /*==========================================================================================*/
+    /*Initialize the spline for 2D interpolation for Il(tau,u)^{k-1} and Ir(tau,u)^{k-1}*/
+    /*both for upward and downward intensities*/
+    /*==========================================================================================*/
+
+    set_spline2d_obj(Spline_I2d_l_upstream, Il_upstream_za, ptr_Iklr[kk - 1].Il_matrix_upstream);
+    set_spline2d_obj(Spline_I2d_l_downstream, Il_downstream_za, ptr_Iklr[kk - 1].Il_matrix_downstream);
+
+    set_spline2d_obj(Spline_I2d_r_upstream, Ir_upstream_za, ptr_Iklr[kk - 1].Ir_matrix_upstream);
+    set_spline2d_obj(Spline_I2d_r_downstream, Ir_downstream_za, ptr_Iklr[kk - 1].Ir_matrix_downstream);
+
+    gsl_spline2d_init(Spline_I2d_l_upstream, array_tau, array_mu, Il_upstream_za, Nstep_tau, Nstep_mu);
+    gsl_spline2d_init(Spline_I2d_l_downstream, array_tau, array_mu, Il_downstream_za, Nstep_tau, Nstep_mu);
+
+    gsl_spline2d_init(Spline_I2d_r_upstream, array_tau, array_mu, Ir_upstream_za, Nstep_tau, Nstep_mu);
+    gsl_spline2d_init(Spline_I2d_r_downstream, array_tau, array_mu, Ir_downstream_za, Nstep_tau, Nstep_mu);
+
+    ptr_data->Spline_I2d_l_upstream = Spline_I2d_l_upstream;
+    ptr_data->Spline_I2d_l_downstream = Spline_I2d_l_downstream;
+
+    ptr_data->Spline_I2d_r_upstream = Spline_I2d_r_upstream;
+    ptr_data->Spline_I2d_r_downstream = Spline_I2d_r_downstream;
+
+    /*=======================================================================*/
+    /*Solve the RTE at the quadrature points*/
+    /*=======================================================================*/
+
+    status = GSL_SUCCESS + 1;
+
+    gsl_set_error_handler_off();
+
+    for (jj = 0; jj < Nstep_mu; jj++)
+    {
+      // printf("\nSolve system for u=%4.3f\n", array_mu[jj]);
+
+      ptr_data->u = array_mu[jj];
+      ptr_data->array_u = array_mu;
+      ptr_data->weights_u = array_weights;
+
+      params = ptr_data;
+
+      /*======================================*/
+      /*Initial boundary conditions*/
+      /*======================================*/
+
+      tau = 0;
+      tau_prev = 0;
+      int idx;
+
+      y[0] = 0;
+      y[1] = 0;
+      y[2] = 0;
+      y[3] = 0;
+
+      /*   h_step=1e-7;
+         gsl_odeiv_evolve * e_rte       = gsl_odeiv_evolve_alloc (4);
+         gsl_odeiv_control * c_rte      = gsl_odeiv_control_y_new (1e-8, 1e-8);
+         gsl_odeiv_step * s_rte         = gsl_odeiv_step_alloc (T, 4);
+   */
+      gsl_odeiv2_system sys = {RTE_Equations, jac, 4, params};
+      gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, step_tau, 1e-4, 1e-4);
+
+      while (tau < 2 * tau0)
+      {
+        status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
+
+        while (status != GSL_SUCCESS)
+        {
+          tau = 0;
+          step_tau = step_tau * 0.75;
+          printf("step_tau now is %lf\n", step_tau);
+
+          // gsl_odeiv2_driver_free(d);
+
+          gsl_odeiv2_driver* d =
+              gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, step_tau, 1e-3, 1e-3);
+
+          status = gsl_odeiv2_driver_apply_fixed_step(d, &tau, step_tau, 5, y);
+        }
+
+        idx = set_array_idx(tau_prev, tau, array_tau, Nstep_tau);
+
+        if (status)
+        {
+          printf("Error during integration: %s\n", gsl_strerror(status));
+          printf("Status %d\n", status);
+          exit(1);
+        }
+
+        // printf("tau_prev=%lf tau=%lf idx %d\n", tau_prev, tau, idx);
+
+        // status = gsl_odeiv_evolve_apply(e_rte, c_rte, s_rte, &sys, &tau, 2*tau0, &h_step, y);
+
+        // printf("tau=%lf  y[0]=%4.3e  y[1]=%4.3e y[2]=%4.3e y[3]=%4.3e idx %d\n", tau, y[0], y[1],
+        // y[2], y[3], idx);
+
+        if (idx > 0)
+        {
+          // printf("tau=%lf  y[0]=%4.3e  y[1]=%4.3e y[2]=%4.3e y[3]=%4.3e idx %d\n", tau, y[0], y[1],
+          //      y[2], y[3], idx);
+
+          ptr_Iklr[kk].Il_matrix_upstream[idx][jj] = y[0];
+          ptr_Iklr[kk].Ir_matrix_upstream[idx][jj] = y[1];
+
+          ptr_Iklr[kk].Il_matrix_downstream[idx][jj] = y[2];
+          ptr_Iklr[kk].Ir_matrix_downstream[idx][jj] = y[3];
+        }
+
+        tau_prev = tau;
+      }
+
+      gsl_odeiv2_driver_free(d);
+
+      // gsl_odeiv_evolve_free (e_rte);
+      // gsl_odeiv_control_free (c_rte);
+      // gsl_odeiv_step_free (s_rte);
+
+    } // End of loop over array_mu
+
+  } // End of loop over k_iter
+
+  Ikl_intensity* ptr_a = 0;
+  Ikr_intensity* ptr_b = 0;
+
+  compute_results(1, k_iter, Nstep_tau, Nstep_mu, array_mu, array_weights, ptr_Iklr, ptr_a, ptr_b);
+
+  exit(1);
+}
+
+/*================================================================================*/
+/*Solution of equations (4.206) and (4.207) of Pomraning 1973*/
+/*================================================================================*/
+
+int RTE_Equations(double s, const double y[], double f[], void* params)
+{
+  double* array_u = ((Iklr_intensity*)params)->array_u;
+  double* weights_u = (double*)((Iklr_intensity*)params)->weights_u;
+  double u = ((Iklr_intensity*)params)->u;
+
+  gsl_spline2d* Spline_I2d_l_upstream = ((Iklr_intensity*)params)->Spline_I2d_l_upstream;
+  gsl_spline2d* Spline_I2d_l_downstream = ((Iklr_intensity*)params)->Spline_I2d_l_downstream;
+
+  gsl_spline2d* Spline_I2d_r_upstream = ((Iklr_intensity*)params)->Spline_I2d_r_upstream;
+  gsl_spline2d* Spline_I2d_r_downstream = ((Iklr_intensity*)params)->Spline_I2d_r_downstream;
+
+  if (s > 2 * tau0 - 1e-5)
+  {
+    s = 2 * tau0;
+  }
+
+  /*==============================================*/
+  /*Equation for I_l  upstream*/
+  /*==============================================*/
+
+  f[0] = -1 / u * y[0] +
+         3 / (8 * u) *
+             (legendre_integration_A(Spline_I2d_l_upstream, s, u, array_u, weights_u) +
+              legendre_integration_A(Spline_I2d_l_downstream, 2 * tau0 - s, u, array_u, weights_u) +
+
+              legendre_integration_B(Spline_I2d_r_upstream, s, u, array_u, weights_u) +
+              legendre_integration_B(Spline_I2d_r_downstream, 2 * tau0 - s, u, array_u, weights_u));
+
+  /*==============================================*/
+  /*Equation for I_r upstream*/
+  /*==============================================*/
+
+  f[1] = -1 / u * y[1] +
+         3 / (8 * u) *
+             (legendre_integration_C(Spline_I2d_l_upstream, s, u, array_u, weights_u) +
+              legendre_integration_C(Spline_I2d_l_downstream, 2 * tau0 - s, u, array_u, weights_u) +
+              legendre_integration_D(Spline_I2d_r_upstream, s, u, array_u, weights_u) +
+              legendre_integration_D(Spline_I2d_r_downstream, 2 * tau0 - s, u, array_u, weights_u));
+
+  /*==============================================*/
+  /*Equation for I_l downstream*/
+  /*==============================================*/
+
+  f[2] = -1 / u * y[2] +
+         3 / (8 * u) *
+             (legendre_integration_A(Spline_I2d_l_downstream, s, u, array_u, weights_u) +
+              legendre_integration_A(Spline_I2d_l_upstream, 2 * tau0 - s, u, array_u, weights_u) +
+
+              legendre_integration_B(Spline_I2d_r_downstream, s, u, array_u, weights_u) +
+              legendre_integration_B(Spline_I2d_r_upstream, 2 * tau0 - s, u, array_u, weights_u));
+
+  /*==============================================*/
+  /*Equation for I_r downstream*/
+  /*==============================================*/
+
+  f[3] = -1 / u * y[3] +
+         3 / (8 * u) *
+             (legendre_integration_C(Spline_I2d_l_downstream, s, u, array_u, weights_u) +
+              legendre_integration_C(Spline_I2d_l_upstream, 2 * tau0 - s, u, array_u, weights_u) +
+
+              legendre_integration_D(Spline_I2d_r_downstream, s, u, array_u, weights_u) +
+              legendre_integration_D(Spline_I2d_r_upstream, 2 * tau0 - s, u, array_u, weights_u));
+
+  // printf("SLURM tau=%lf u=%5.4f f[0]=%5.4e f[1]=%5.4e f[2]=%5.4e f[3]=%5.4e\n", s, u, f[0], f[1], f[2], f[3]);
+
+  return GSL_SUCCESS;
+}
+
+/*====================================================================*/
+/*Value of the I^k(tau,u) function for k=0 and seed photons*/
+/*at the base of the slab (tau=0)*/
+/*====================================================================*/
+
+double I0_upward(double tau, double u)
+{
+  double value;
+
+  if (tau == 0)
+  {
+    value = 0;
+  }
+  else
+  {
+    value = 1 / 2. * 1 / u * exp(-tau / u);
+  }
+
+  return value;
+}
+
+/*====================================================================*/
+/*Value of the I^k(tau,u) function for k=0 and seed photons*/
+/*uniformly distributed*/
+/*====================================================================*/
+
+double I0_uniform_updown(double tau, double u)
+{
+  double value;
+
+  value = 1 / 2. * (1 - exp(-tau / u));
+
+  return value;
+}
+
+/*====================================================================*/
+/*Value of the I^k(tau,u) function for k=0 and seed photons*/
+/*uniformly distributed*/
+/*====================================================================*/
+
+double I0_center_updown(double tau, double u)
+{
+  double value;
+
+  if (tau >= tau0)
+  {
+    value = 1 / 2. * 1 / u * exp(-(tau - tau0) / u);
+  }
+  else
+  {
+    value = 0;
+  }
+
+  return value;
+}
+
+/*==============================================================================*/
+void set_spline2d_obj(gsl_spline2d* Spline_I2d, double* za, double** I_funct)
+{
+  int tt, qq;
+
+  for (tt = 0; tt < Nstep_tau; tt++)
+  {
+    for (qq = 0; qq < Nstep_mu; qq++)
+    {
+      /*printf("tt=%d qq=%d funct=%4.3e\n", tt, qq, I_funct[tt][qq]); */
+      gsl_spline2d_set(Spline_I2d, za, tt, qq, I_funct[tt][qq]);
+    }
+  }
+}
+
+int set_array_idx(double tau_prev, double tau, double* array, int nstep)
+{
+  int idx, ii;
+
+  idx = -1;
+
+  if (tau_prev == array[0])
+  {
+    idx = 0;
+  }
+  else if (tau == array[nstep - 1])
+  {
+    idx = nstep - 1;
+  }
+  else
+  {
+    for (ii = 1; ii < nstep; ii++)
+    {
+      if (tau_prev <= array[ii] && tau >= array[ii])
+      {
+        idx = ii;
+        // printf("Sono qui tau_prev %4.3f tau %4.3f array[ii] %4.3f idx %d\n", tau_prev, tau,
+        //		array[ii], idx);
+
+        break;
+      }
+    }
+  }
+
+  return idx;
+}