Added the albedo at the base of the slab as a parameter for the Montecarlo method

.clang-format

+Language: Cpp
+AccessModifierOffset: -4  
+AlignEscapedNewlinesLeft: true  
+AlignTrailingComments: true  
+AllowAllParametersOfDeclarationOnNextLine: true  
+AllowShortBlocksOnASingleLine: true  
+AllowShortFunctionsOnASingleLine: None  
+AllowShortIfStatementsOnASingleLine: false  
+AllowShortLoopsOnASingleLine: false  
+AlwaysBreakBeforeMultilineStrings: true  
+AlwaysBreakTemplateDeclarations: true  
+BinPackParameters: false  
+BreakBeforeBinaryOperators: false  
+BreakBeforeBraces: Allman  
+BreakBeforeTernaryOperators: true  
+BreakConstructorInitializersBeforeComma: false  
+ColumnLimit: 100
+CommentPragmas: ''  
+ConstructorInitializerAllOnOneLineOrOnePerLine: true  
+ConstructorInitializerIndentWidth: 4  
+ContinuationIndentWidth: 4  
+DerivePointerAlignment: false  
+DisableFormat: false  
+ExperimentalAutoDetectBinPacking: false  
+IndentCaseLabels: false  
+IndentWidth: 2  
+IndentWrappedFunctionNames: false  
+IndentFunctionDeclarationAfterType: false  
+MaxEmptyLinesToKeep: 1  
+KeepEmptyLinesAtTheStartOfBlocks: false  
+NamespaceIndentation: None  
+PenaltyBreakBeforeFirstCallParameter: 1  
+PenaltyBreakComment: 300  
+PenaltyBreakString: 1000  
+PenaltyBreakFirstLessLess: 120  
+PenaltyExcessCharacter: 1  
+PenaltyReturnTypeOnItsOwnLine: 1000  
+PointerAlignment: Left  
+SpaceBeforeAssignmentOperators: true  
+SpaceBeforeParens: ControlStatements  
+SpaceInEmptyParentheses: false  
+SpacesBeforeTrailingComments: 1  
+SpacesInAngles: false  
+SpacesInCStyleCastParentheses: false  
+SpacesInContainerLiterals: true  
+SpacesInParentheses: false  
+Cpp11BracedListStyle: true  
+Standard: Cpp11  
+TabWidth: 4  
+UseTab: Never 

compute_stokes.c

@@ -91,6 +91,12 @@ void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_para
     double* Pdeg = (double*)malloc(nstepangles * sizeof(double));
     dat_integral_polar = fopen(integral, "w");
 
+    fprintf(dat_integral_polar, "!Col1: cos(theta)\n");
+    fprintf(dat_integral_polar, "!Col2: Pdeg\n");
+    fprintf(dat_integral_polar, "!Col3: Q\n");
+    fprintf(dat_integral_polar, "!Col4: U\n");
+    fprintf(dat_integral_polar, "!Col5: Itot/I(u=1)*1/u\n\n");
+
     for (jj = 0; jj < nstepangles; jj++)
     {
       u_center = (array_ubounds[jj] + array_ubounds[jj + 1]) / 2;

main_program.c

@@ -2,24 +2,29 @@
 #include <mpi.h>
 #include <stdio.h>
 
+char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method);
+
 int seed;
 double kte;
 double ktbb;
 double tau_c;
 double disktau;
+double albedobase;
+double gmin;
+double gmax;
+double p1;
+
 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;
+int Albedoflag;
 
 char* outspec;
 char* polarfile;
@@ -27,8 +32,6 @@ 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;
@@ -57,7 +60,7 @@ int main(int argc, char* argv[])
     printf("Shared parameters\n\n");
     printf("disktau=vaue: slab optical depth\n\n");
 
-    printf("\nSeed photons distributions: three options are possible\n");
+    printf("\nSeed photons distributions: three options are possible with parameter seed=N\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");
@@ -85,6 +88,7 @@ int main(int argc, char* argv[])
 
   char* tau_char;
   char* seed_char;
+  char* albedo_char;
 
   for (t = 1; t < argc; t++)
   {
@@ -155,6 +159,13 @@ int main(int argc, char* argv[])
     {
       ktbb = numpar(argv[t]);
     }
+
+    if (strstr(argv[t], "albedobase=") != NULL)
+    {
+      albedobase = numpar(argv[t]);
+      albedo_char = stringpar(argv[t]);
+      Albedoflag = 1;
+    }
   }
 
   /*=========================================================*/
@@ -167,26 +178,6 @@ int main(int argc, char* argv[])
     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;
@@ -213,8 +204,8 @@ int main(int argc, char* argv[])
 
   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);
+    diffusion_file = set_filename("diffusion_ode_tau", tau_char, seed_char, albedo_char, "ode");
+    polarization_file = set_filename("intensity_ode_tau", tau_char, seed_char, albedo_char, "ode");
 
     if (kiter == 0)
     {
@@ -226,8 +217,8 @@ int main(int argc, char* argv[])
   }
   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);
+    polarization_file = set_filename("intensity_st85_tau", tau_char, seed_char, albedo_char, "st85");
+    diffusion_file = set_filename("diffusion_st85_tau", tau_char, seed_char, albedo_char, "st85");
 
     if (kiter == 0)
     {
@@ -239,8 +230,8 @@ int main(int argc, char* argv[])
   }
   else if (strcmp(method, "mc") == 0)
   {
-    MPI_Init(&argc, &argv); 
-    init_mpi(); 
+    MPI_Init(&argc, &argv);
+    init_mpi();
 
     if (polardisk == 0x0 && seed == 3 && rank == 0)
     {
@@ -253,6 +244,18 @@ int main(int argc, char* argv[])
       polardisk = "n";
     }
 
+    if (Albedoflag != 1)
+    {
+      printf("Please select the albedo at the base of the slab with albedobase=value\n");
+      exit(1);
+    }
+
+    if (albedobase < 0 || albedobase > 1)
+    {
+      printf("Error: albedo must be in the range 0-1\n");
+      exit(1);
+    }
+
     if (nph == 0)
     {
       printf("Please select the number of photons for MC simulation with parametere nph=nph\n");
@@ -279,7 +282,13 @@ int main(int argc, char* argv[])
       exit(1);
     }
 
+    diffusion = set_filename("mc_diffusion_tau", tau_char, seed_char, albedo_char, "mc");
+    outspec = set_filename("mc_spectrum_tau", tau_char, seed_char, albedo_char, "mc");
+    integral = set_filename("mc_integralpolar_tau", tau_char, seed_char, albedo_char, "mc");
+    polarfile = set_filename("mc_energypolar_tau", tau_char, seed_char, albedo_char, "mc");
+
     status = slab_mc(nph, seed);
+
     // check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);
   }
 
@@ -295,7 +304,7 @@ int main(int argc, char* argv[])
   return (0);
 }
 
-char* set_filename(char* root, char* tau_char, char* seed_char)
+char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method)
 {
   // root="intensity_ode";
 
@@ -303,6 +312,7 @@ char* set_filename(char* root, char* tau_char, char* seed_char)
   char* tmp_a;
   char* tmp_b;
   char* tmp_c;
+  char* tmp_d;
 
   tmp_a = stringconcat(root, tau_char);
 
@@ -316,9 +326,20 @@ char* set_filename(char* root, char* tau_char, char* seed_char)
 
   // printf("tmp_c = %s\n", tmp_c);
 
-  filename = stringconcat(tmp_c, ".dat");
+  if (strcmp(method, "mc") == 0)
+  {
+    tmp_d = stringconcat(tmp_c, stringconcat("_A", albedo));
+
+    filename = stringconcat(tmp_d, ".qdp");
+  }
+  else
+  {
+    filename = stringconcat(tmp_c, ".qdp");
+  }
 
   // printf("Final name %s\n\n", filename);
 
   return filename;
 }
+
+/*==============================================================================*/

mc_slab.c

@@ -2,6 +2,8 @@
 
 #define NSC_MAX 1000
 
+extern double albedobase;
+
 /*==================================================================*/
 /*Function prototypes*/
 /*==================================================================*/
@@ -40,17 +42,14 @@ 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;
-
+gsl_root_fsolver* s_maxw;
 
 /*==================================================================*/
 
@@ -92,13 +91,12 @@ int slab_mc(int nphot, int seed)
 
   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.;
@@ -303,6 +301,8 @@ int slab_mc(int nphot, int seed)
         Magfield_lab[2] = Magfield_lab_new[2];
       }
       else
+
+      /*Photon potentially escaping from the slab	  */
       {
         if (k_lab[2] > 0)
         {
@@ -313,18 +313,46 @@ int slab_mc(int nphot, int seed)
           collect_photons(E_new_lab, k_lab, array_ene, array_counts, 0);
 
           nph_esc++;
+          FlagLocation = 0;
         }
         else
         {
-          // printf("Photon %d escaping\n", ii);
-        }
+          csi = clock_random();
 
-        FlagLocation = 0;
-      }
+          /*Photon absorbed and newly re-emitted with probability A=albedobase*/
+          if (csi <= albedobase)
+          {
+            // printf("Photon %d reflected\n", ii);
+
+            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;
+
+            pos[0] = 0;
+            pos[1] = 0;
+            pos[2] = 0;
+
+            unpolarized_seedphot(k_lab, Efield_lab, Magfield_lab);
+            PolarFlag = 0;
+            E_lab = bb_spec(ktbb);
+          }
+
+          /*Photon absorbed*/
+          else
+          {
+            FlagLocation = 0;
+          }
+        } // Close the else with condition k_lab[2] < 0
+
+      } // End of confition  if (tau > eps)... else...
 
     } // end of while (FlagLocation)
 
-    // newphot:;
+    
   }
 
   /*==================================================================*/
@@ -426,9 +454,9 @@ void comptonization(double E_lab,
 
   versor_boost(k_lab, el_vel_vect, k_erf, 1);
 
-  double beta_new=subrel_maxwellian (kte);
+  double beta_new = subrel_maxwellian(kte);
 
-  //printf("pluto %lf %lf\n", beta_el, beta_new);
+  // printf("pluto %lf %lf\n", beta_el, beta_new);
 
   elmag_lorentztransform(Efield_lab, Magfield_lab, el_vel_vect, ElectField_old_erf, MagField_Erf, 1);
 

run_mc.sh

+array_tau=(0.2 0.5 1 2 5 8)
+array_seed=(2 3)
+array_albedo=(0 0.5 1)
+
+ii="0"
+jj="0"
+kk="0"
+
+nph=10000
+
+#Loop over tau
+while [ $ii -lt ${#array_tau[@]}  ]
+do
+
+     #Loop over seed distribution
+     jj="0"
+     while [ $jj -lt ${#array_seed[@]}  ]
+     do
+         
+         #Loop over slabrmax
+         kk="0"
+         while [ $kk -lt ${#array_albedo[@]}  ]
+         do
+
+
+	     command="start_iteration method=mc disktau=${array_tau[ii]}  seed=${array_seed[jj]} polardisk=n nph=$nph albedobase=${array_albedo[kk]}"
+
+	     echo $command
+
+	     $command
+
+	           kk=$[$kk+1]
+         done
+
+	           jj=$[$jj+1]
+     done
+
+     	           ii=$[$ii+1]
+        done
+
+	 
+
+