From 4add7111cc892bed47a3c5f46ad06f45c4621e27 Mon Sep 17 00:00:00 2001
From: Ruben Farinelli <ruben@MacBook-Air-di-Ruben.local>
Date: Mon, 25 Sep 2023 17:19:37 +0200
Subject: [PATCH] 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 | 27 ++++++++++++++++++++++--
compute_stokes.c | 2 +-
main_program.c | 8 +++----
mc_slab.c | 21 +++++++++++++++---
run_ode.sh | 54 +++++++++++++++++++++++++++++++++++++++++++++++
5 files changed, 102 insertions(+), 10 deletions(-)
create mode 100755 run_ode.sh
diff --git a/compute_results.c b/compute_results.c
index 97c496f..2c899a5 100644
--- a/compute_results.c
+++ b/compute_results.c
@@ -80,8 +80,13 @@ 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++)
{
fprintf(dat, "%lf %5.4e %5.4e\n", array_mu[jj], Pdeg[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);
diff --git a/compute_stokes.c b/compute_stokes.c
index 4c46861..7dbc198 100644
--- a/compute_stokes.c
+++ b/compute_stokes.c
@@ -140,7 +140,7 @@ void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_para
}
fprintf(dat_integral_polar, "%lf %lf %lf %lf %5.4f\n", u_center, Pdeg[jj], Qtot / Itot,
- Utot / Itot, 1 / u_center * (Itot / Itot_u1));
+ Utot / Itot, 1/u_center*(Itot / Itot_u1));
}
fclose(dat_integral_polar);
diff --git a/main_program.c b/main_program.c
index d29571b..7519636 100644
--- a/main_program.c
+++ b/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);
diff --git a/mc_slab.c b/mc_slab.c
index 5cd29c7..70204ea 100644
--- a/mc_slab.c
+++ b/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++)
+ {
+ if (photon_diffusion[ii] > Fval_max)
+ {
+ Fval_max = photon_diffusion[ii];
+ }
+ }
+
for (ii = 0; ii < 300; ii++)
{
- fprintf(dat_diffusion, " %d %d\n", ii, photon_diffusion[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;
diff --git a/run_ode.sh b/run_ode.sh
new file mode 100755
index 0000000..94f4176
--- /dev/null
+++ b/run_ode.sh
@@ -0,0 +1,54 @@
+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
+
+
+
+
--
GitLab