diff --git a/compute_stokes.c b/compute_stokes.c
index 061a6e634d2850bc23b66d2fa2d598b2343d36b3..d980dc827b08b353146bf18b0130e62925f8aae6 100644
--- a/compute_stokes.c
+++ b/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))
{
diff --git a/main_program.c b/main_program.c
index 20705e87e328103511e5776ea7608a2627760475..00308f304e6e42b20364d4d28549499a47104e26 100644
--- a/main_program.c
+++ b/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");
diff --git a/mc_slab.c b/mc_slab.c
index 2a658c6883c604c2844d3f356d9523e5a2968d79..8cdf0c6c22c518e1e5e64ade296a78dd9baf5cec 100644
--- a/mc_slab.c
+++ b/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,
diff --git a/plot_data.py b/plot_data.py
new file mode 100755
index 0000000000000000000000000000000000000000..20321ac7979fa390ed50efae456cab3adc606a76
--- /dev/null
+++ b/plot_data.py
@@ -0,0 +1,53 @@
+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)
diff --git a/solve_rte.c b/solve_rte.c
index 20dd9df6b35fd7e89b7baec30b9e1846ee4f1766..e474d1c9d9164b9553d1b323a5367b4fe3b622ff 100644
--- a/solve_rte.c
+++ b/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)*/