Corrected bug in the MC routine for high values of the optical depth

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}
+CFLAGS   = -g  -O3 -Wall ${INCLUDE}
 
 
 LOCAL_DIR               = ${PWD}
@@ -20,7 +20,7 @@ MC_SYNC                 = ${MC_INSTDIR}/../../gr_code/synchrotron
 
 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  ${LOCAL_DIR}/simul_sync.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_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
 

main_program.c

-#include "iteration_functions.h"
 #include "globals.h"
+#include "iteration_functions.h"
 #include <mpi.h>
 #include <stdio.h>
 
@@ -9,11 +9,9 @@ int seed;
 double tau_c;
 double albedobase;
 
-
 char* polarization_file;
 char* diffusion_file;
 
-
 // int nph;
 // int rank;
 
@@ -152,6 +150,16 @@ int main(int argc, char* argv[])
       ktbb = numpar(argv[t]);
     }
 
+    if (strstr(argv[t], "emin=") != NULL)
+    {
+      emin = numpar(argv[t]);
+    }
+
+    if (strstr(argv[t], "emax=") != NULL)
+    {
+      emax = numpar(argv[t]);
+    }
+
     if (strstr(argv[t], "albedobase=") != NULL)
     {
       albedobase = numpar(argv[t]);
@@ -180,6 +188,16 @@ int main(int argc, char* argv[])
     ktbb = 1;
   }
 
+  if (!emin)
+  {
+    emin = 0.1;
+  }
+
+  if (!emax)
+  {
+    emax = 50.;
+  }
+
   /*=========================================================*/
 
   if (!disktau || disktau == 0)
@@ -274,23 +292,37 @@ int main(int argc, char* argv[])
       exit(1);
     }
 
+    if (diffusion == 0x0)
+    {
       diffusion = set_filename("diffusion_mc_tau", tau_char, seed_char, albedo_char, "mc");
+    }
+
+    if (outspec == 0x0)
+    {
       outspec = set_filename("spectrum_mc_tau", tau_char, seed_char, albedo_char, "mc");
+    }
+
+    if (integral == 0x0)
+    {
       integral = set_filename("integralpolar_mc_tau", tau_char, seed_char, albedo_char, "mc");
+    }
+
+    if (polarfile == 0x0)
+    {
       polarfile = set_filename("energypolar_mc_tau", tau_char, seed_char, albedo_char, "mc");
+    }
 
     status = slab_mc(nph, seed);
-    
   }
 
   else
   {
-    printf(
-        "\nPlease select method=ode (solve ODE) or method=mc (Montecarlo)\n");
+    printf("\nPlease select method=ode (solve ODE) or method=mc (Montecarlo)\n");
     return (1);
   }
 
-  if (strcmp(method, "mc") == 0) MPI_Finalize();
+  if (strcmp(method, "mc") == 0)
+    MPI_Finalize();
 
   return (0);
 }

mc_slab.c

@@ -17,7 +17,6 @@ void comptonization(double E_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);
 
@@ -28,14 +27,12 @@ double distance_top(double z0, double kz, double H);
 // double obsmindeg;
 // double obsmaxdeg;
 
-
 // double* array_angles;
 // double* array_ubounds;
 
 // double** array_gb;
 // double** ptr_2colfile;
 
-
 const gsl_root_fsolver_type* T_phi;
 gsl_root_fsolver* s_phi;
 
@@ -59,7 +56,6 @@ int slab_mc(int nphot, int seed)
   int photon_diffusion[NSC_MAX];
   int status;
 
-
   double E_lab;
   double E_new_lab;
   double k_lab[3];
@@ -105,7 +101,7 @@ int slab_mc(int nphot, int seed)
   array_gb = read_matrix(matrix_gb);
 
   double array_ene[NSTEP_ENE];
-  make_array_energy(array_ene, 0.01, 20, NSTEP_ENE);
+  make_array_energy(array_ene, emin, emax, NSTEP_ENE);
 
   double** array_counts;
   array_counts = dmatrix(0, NSTEP_ENE - 1, 0, nstepangles - 1);
@@ -139,7 +135,6 @@ int slab_mc(int nphot, int seed)
     setup_gsl_objects(radiation_field, MC_SLAB);
   }
 
-
   Pdeg = 0;
   nph_esc = 0;
   cos_theta = 0;
@@ -190,7 +185,10 @@ int slab_mc(int nphot, int seed)
       cos_theta = 2 * clock_random() - 1;
     }
     else if (seed == 2)
+    {if (!emin)
   {
+	emin=0.1;
+  }
       pos[2] = clock_random() * Hphot;
       cos_theta = 2 * clock_random() - 1;
     }
@@ -224,7 +222,6 @@ int slab_mc(int nphot, int seed)
     k_lab[1] = sin_theta * sin(phi);
     k_lab[2] = cos_theta;
 
-
     if (DiskPolarizationFlag == FROMFILE)
     {
       electricfield_fromfile(k_lab, Efield_lab, Magfield_lab, &Pdeg);
@@ -280,7 +277,6 @@ int slab_mc(int nphot, int seed)
           Efield_lab[jj] = Efield_lab_new[jj];
           Magfield_lab[jj] = Magfield_lab_new[jj];
         }
-
       }
       else
 
@@ -290,10 +286,16 @@ int slab_mc(int nphot, int seed)
         {
           // printf("Photon %d escaping after %d scattering\n", ii, n_sc);
 
-          if (n_sc==0) n_zero_sc++;
+          if (n_sc == 0)
+            n_zero_sc++;
 
           compute_stokes(E_lab, k_lab, Efield_lab, struct_stokes, array_ubounds, 0);
+
+          if (n_sc < NSC_MAX)
+          {
             photon_diffusion[n_sc]++;
+          }
+
           collect_photons(E_new_lab, k_lab, array_ene, array_counts, 0);
 
           nph_esc++;
@@ -366,9 +368,10 @@ int slab_mc(int nphot, int seed)
 
   MPI_Comm_size(MPI_COMM_WORLD, &nproc);
 
-
   /*==================================================================*/
 
+  // printf("Rank %d si trova qui %d\n", rank, __LINE__);
+
   if (rank == 0)
   {
     compute_stokes(1, k_lab, Efield_lab, struct_stokes_average, array_ubounds, 1);
@@ -396,24 +399,31 @@ int slab_mc(int nphot, int seed)
     }
 
     fclose(dat_diffusion);
-  }
-
-  collect_photons(0, k_lab, array_ene, reduce_array_counts, 1);
 
     dat_energy_polar = fopen(polarfile, "w");
+    if (dat_energy_polar == NULL)
+    {
+      fprintf(stderr, "Errore: fopen fallita su %s\n", polarfile);
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
 
     write_polardata(struct_stokes_average, ene_bound, dat_energy_polar);
-
     fclose(dat_energy_polar);
     
+    collect_photons(0, k_lab, array_ene, reduce_array_counts, 1);
+  }
+
+  
+
+  if (s_phi != NULL)
     gsl_root_fsolver_free(s_phi);
+  if (s_disk != NULL)
     gsl_root_fsolver_free(s_disk);
 
   return 0;
 
 } // End of function
 
-
 void comptonization(double E_lab,
                     double* k_lab,
                     double* Efield_lab,
@@ -424,7 +434,6 @@ void comptonization(double E_lab,
                     double* Magfield_lab_new)
 
 {
-
   int jj;
   int status;
   int PolarFlag;
@@ -444,7 +453,6 @@ void comptonization(double E_lab,
   double ElectField_old_erf[3];
   double MagField_Erf[3];
 
-
   PolarFlag = 0;
   Pdeg = 0;
 
@@ -462,7 +470,6 @@ void comptonization(double E_lab,
 
   E_erf = energy_lorentz_boost(E_lab, el_vel_vect, k_lab, 1);
 
-
   csi = clock_random();
 
   costheta_erf = (-1 + pow(-2 + 4 * csi + sqrt(5 - 16 * csi + 16 * pow(csi, 2)), 2 / 3.)) /
@@ -476,7 +483,6 @@ void comptonization(double E_lab,
 
   E_new_erf = new_energy_erf(E_erf, costheta_erf);
 
-
   phik = sample_azim_angle(E_erf, E_new_erf, costheta_erf, 1);
 
   //================================================================== */
@@ -499,7 +505,6 @@ void comptonization(double E_lab,
   PolarFlag = 1;
   status = depolarization_probability(E_erf, E_new_erf, costheta_erf, phik, &Pdeg, &PolarFlag);
 
-
   if (status == 1)
   {
     exit(1);