For both the Montecarlo and numerical method, made a change so that the P(u) distribution function

is printed normalized to its maximum value

compute_results.c

@@ -80,6 +80,11 @@ void compute_results(int method,
   /*Print I(u) normalized at the value u=1*/
   /*===========================================================================*/
 
+  fprintf(dat, "!Col1: cos(theta)\n");
+  fprintf(dat, "!Col2: Pdeg\n");
+  fprintf(dat, "!Col3: I(u)/I(u=1)\n");
+  
+  
   if (disktau < 4)
   {	  
     for (jj = 0; jj < Nstep_mu; jj++)
@@ -89,6 +94,16 @@ void compute_results(int method,
     }
   }
   else
+	  
+	  fprintf(dat, "!Pdeg band 1\n");
+	    fprintf(dat, "!Pdeg band 2\n");
+	    fprintf(dat, "!Pdeg band 3\n");
+	    fprintf(dat, "!Pdeg band 4\n");
+	      fprintf(dat, "!Pdeg band 5\n");
+	      
+	     
+	   
+	  
   {
     for (jj = 0; jj < Nstep_mu; jj++)
     {
@@ -125,6 +140,7 @@ void compute_results(int method,
 
   double* diffusion = malloc(k_iter * sizeof(double));
   double sum_k;
+  double Fval_max = -10;
 
   for (kk = 0; kk < k_iter; kk++)
   {
@@ -144,7 +160,14 @@ void compute_results(int method,
     }
 
     diffusion[kk] = sum_k;
-    fprintf(dat, "%d  %5.4e\n", kk, diffusion[kk]);
+
+    if (diffusion[kk] > Fval_max)
+      Fval_max = diffusion[kk];
+  }
+
+  for (kk = 0; kk < k_iter; kk++)
+  {
+    fprintf(dat, "%d  %5.4f\n", kk, diffusion[kk] / Fval_max);
   }
 
   fclose(dat);

main_program.c

@@ -282,10 +282,10 @@ 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");
+    diffusion = set_filename("diffusion_mc_tau", tau_char, seed_char, albedo_char, "mc");
+    outspec = set_filename("spectrum_mc_tau", tau_char, seed_char, albedo_char, "mc");
+    integral = set_filename("integralpolar_mc_tau", tau_char, seed_char, albedo_char, "mc");
+    polarfile = set_filename("energypolar_mc_tau", tau_char, seed_char, albedo_char, "mc");
 
     status = slab_mc(nph, seed);
 

mc_slab.c

@@ -61,6 +61,8 @@ int slab_mc(int nphot, int seed)
   int FlagLocation;
   int n_sc;
   int nph_esc;
+  int Fval_max;
+  int nsc_max;
   int photon_diffusion[NSC_MAX];
   int status;
   double E_lab;
@@ -351,8 +353,6 @@ int slab_mc(int nphot, int seed)
       } // End of confition  if (tau > eps)... else...
 
     } // end of while (FlagLocation)
-
-    
   }
 
   /*==================================================================*/
@@ -388,12 +388,23 @@ int slab_mc(int nphot, int seed)
 
     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);
 
+    Fval_max=-10;
     for (ii = 0; ii < 300; ii++)
     {
-      fprintf(dat_diffusion, " %d  %d\n", ii, photon_diffusion[ii]);
+      if (photon_diffusion[ii] > Fval_max)
+      {
+        Fval_max = photon_diffusion[ii];
+      }
+    }
+
+    for (ii = 0; ii < 300; ii++)
+    {
+      fprintf(dat_diffusion, " %d  %3.2f\n", ii, (double)photon_diffusion[ii]/Fval_max);
     }
 
     fclose(dat_diffusion);
@@ -545,6 +556,8 @@ void comptonization(double E_lab,
   }
 }
 
+/*======================================================================*/
+
 double distance_top(double z_pos, double kz, double H)
 {
   double t;
@@ -554,6 +567,8 @@ double distance_top(double z_pos, double kz, double H)
   return t;
 }
 
+/*======================================================================*/
+
 double distance_bottom(double z_pos, double kz)
 {
   double t;

run_ode.sh

+array_tau=(0.2 0.5 1 2 5 8)
+array_seed=(2 3)
+array_albedo=(1)
+
+ii="0"
+jj="0"
+kk="0"
+
+#Loop over tau
+
+while [ $ii -lt ${#array_tau[@]}  ]
+do
+
+    #Loop over seed distribution
+    
+     jj="0"
+     while [ $jj -lt ${#array_seed[@]}  ]
+     do
+         
+         #Loop over albedo (numerical solution does not depend on it actually)
+	 
+         kk="0"
+         while [ $kk -lt ${#array_albedo[@]}  ]
+         do
+
+	     if (( $(echo "${array_tau[ii]} <= 2" | bc -l) )); then
+
+		 kiter=20
+
+	       else
+	     
+	     
+		   kiter=$(echo "scale=2; ${array_tau[ii]}^2 *3" | bc)
+
+	      fi
+
+	     command="start_iteration method=ode disktau=${array_tau[ii]}  seed=${array_seed[jj]} polardisk=n albedobase=${array_albedo[kk]} kiter=$kiter"
+
+	     echo $command
+
+	     #$command
+
+	           kk=$[$kk+1]
+         done
+
+	           jj=$[$jj+1]
+     done
+
+     	           ii=$[$ii+1]
+        done
+
+	 
+
+