Added exponential distribution of seed photons for the ODE method

compute_stokes.c

@@ -55,9 +55,9 @@ void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_para
               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].array_I[jj] = ptr_stokes[ii].array_I[jj] + 1;
+            ptr_stokes[ii].array_Q[jj] = ptr_stokes[ii].array_Q[jj] + Qs;
+            ptr_stokes[ii].array_U[jj] = ptr_stokes[ii].array_U[jj] + Us;
             ptr_stokes[ii].counter[jj]++;
 
             FlagFound = TRUE;
@@ -103,11 +103,11 @@ void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_para
 
       for (ii = 0; ii < POLAR_NBOUNDS - 1; ii++)
       {
-        Itot_allbands = Itot_allbands + ptr_stokes[ii].array_Is[jj];
+        Itot_allbands = Itot_allbands + ptr_stokes[ii].array_I[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];
+        Qtot = Qtot + ptr_stokes[ii].array_Q[jj];
+        Utot = Utot + ptr_stokes[ii].array_U[jj];
+        Itot = Itot + ptr_stokes[ii].array_I[jj];
 
         if (isnan(Qtot))
         {

main_program.c

@@ -58,7 +58,8 @@ int main(int argc, char* argv[])
     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");
+    printf("3: photons at the base (valid only for ode and mc methods)\n");
+    printf("4: photons with exponential distribution (valid only for ode and mc methods)\n\n");
 
     printf("Parameters for iteration algorithms\n");
 

mc_slab.c

@@ -341,13 +341,13 @@ int slab_mc(int nphot, int seed)
 
   for (ii = 0; ii < POLAR_NBOUNDS - 1; ii++)
   {
-    MPI_Reduce(struct_stokes[ii].array_Is, struct_stokes_average[ii].array_Is, nstepangles,
+    MPI_Reduce(struct_stokes[ii].array_I, struct_stokes_average[ii].array_I, nstepangles,
                MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
 
-    MPI_Reduce(struct_stokes[ii].array_Qs, struct_stokes_average[ii].array_Qs, nstepangles,
+    MPI_Reduce(struct_stokes[ii].array_Q, struct_stokes_average[ii].array_Q, nstepangles,
                MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
 
-    MPI_Reduce(struct_stokes[ii].array_Us, struct_stokes_average[ii].array_Us, nstepangles,
+    MPI_Reduce(struct_stokes[ii].array_U, struct_stokes_average[ii].array_U, nstepangles,
                MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
 
     MPI_Reduce(struct_stokes[ii].counter, struct_stokes_average[ii].counter, nstepangles, MPI_INT,

plot_data.py

+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+def read_data(file_path):
+    """
+    Read the data from the input file and extract X, Y, and Z values.
+    Assumes the data format is X Y Z in each line.
+    """
+    X, Y, Z = np.loadtxt(file_path, unpack=True)
+    return X, Y, Z
+
+def plot_3d_histogram(X, Y, Z, bins_x, bins_y):
+    """
+    Plot a 3D histogram using bar3d with square bins.
+    """
+    H, xedges, yedges = np.histogram2d(X, Y, bins=[bins_x, bins_y])
+    X, Y = np.meshgrid(xedges[:-1], yedges[:-1])
+
+    # Set the widths and depths of the bars to make them square
+    dx = (xedges[1] - xedges[0])
+    dy = (yedges[1] - yedges[0])
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+
+    # Flatten the arrays
+    X = X.ravel()
+    Y = Y.ravel()
+    Z = np.zeros_like(X)
+    dx = np.full_like(X, dx)
+    dy = np.full_like(X, dy)
+
+    ax.bar3d(X, Y, Z, dx, dy, H.ravel(), shade=True)
+    ax.set_xlabel('X')
+    ax.set_ylabel('Y')
+    ax.set_zlabel('Frequency')
+    ax.set_title('3D Histogram of Z(x, y)')
+    plt.show()
+
+if __name__ == "__main__":
+    # Path to the input file (replace with your file path)
+    file_path = 'data.txt'
+
+    # Read data from the input file
+    X, Y, Z = read_data(file_path)
+
+    # Number of bins for X and Y axes
+    bins_x = 40  # Adjust as needed
+    bins_y = 40  # Adjust as needed
+
+    # Plot the 3D histogram
+    plot_3d_histogram(X, Y, Z, bins_x, bins_y)

solve_rte.c

@@ -4,6 +4,9 @@
 double I0_center_updown(double tau, double u);
 double I0_upward(double tau, double u);
 double I0_uniform_updown(double tau, double u);
+double I0_exp_upward(double tau, double u);
+double I0_exp_downward(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);
@@ -117,6 +120,19 @@ int solve_rte(int k_iter, int seed_distribution, int Nstep_tau)
         ptr_Iklr[0].Il_matrix_downstream[ii][jj] = 0;
         ptr_Iklr[0].Ir_matrix_downstream[ii][jj] = 0;
       }
+      
+      else if (seed_distribution == SEED_EXP)
+      {
+        ptr_Iklr[0].Il_matrix_upstream[ii][jj] = I0_exp_upward(array_tau[ii], array_mu[jj]);
+        ptr_Iklr[0].Ir_matrix_upstream[ii][jj] = I0_exp_upward(array_tau[ii], array_mu[jj]);
+
+        ptr_Iklr[0].Il_matrix_downstream[ii][jj] = I0_exp_downward(array_tau[ii], array_mu[jj]);
+        ptr_Iklr[0].Ir_matrix_downstream[ii][jj] = I0_exp_downward(array_tau[ii], array_mu[jj]);
+        
+        
+      }
+ 
+      
       else
       {
         printf(
@@ -361,6 +377,38 @@ int RTE_Equations(double s, const double y[], double f[], void* params)
   return GSL_SUCCESS;
 }
 
+
+/*====================================================================*/
+/*Value of the I^k(tau,u) function for k=0 and seed photons*/
+/*with exponential distribution for u > 0 and u < 0*/
+/*====================================================================*/
+
+double I0_exp_upward(double tau, double u)
+{
+	double value;
+	
+	value= ((-1 + exp(((-1 + u) * tau) / u))) / 
+           (2.0 * exp(tau) * (-1 + u));
+	
+	return value;
+	
+}
+
+
+double I0_exp_downward(double tau, double u)
+{
+	double value;
+	  
+   	value=(exp(-2 * tau0) * (exp(tau) - exp(-tau / u))) / 
+           (2.0 * (1 + u));
+    
+          
+	return value;
+	
+}
+
+
+
 /*====================================================================*/
 /*Value of the I^k(tau,u) function for k=0 and seed photons*/
 /*at the base of the slab (tau=0)*/