From ec06623c543235f4c84ae0c95b60107bf89ba66d Mon Sep 17 00:00:00 2001
From: Ruben Farinelli <ruben.farinelli@inaf.it>
Date: Tue, 28 Nov 2023 16:24:14 +0100
Subject: [PATCH] Wrote more detailed README.txt
---
README.txt | 169 +++++++++++++++++++++++++++++++++++++++++++++++++
main_program.c | 16 ++---
make2tar.py | 8 +--
ode_example.c | 27 --------
4 files changed, 177 insertions(+), 43 deletions(-)
mode change 100644 => 100755 make2tar.py
delete mode 100644 ode_example.c
diff --git a/README.txt b/README.txt
index 35855ba..092b891 100644
--- a/README.txt
+++ b/README.txt
@@ -1,3 +1,172 @@
The code computes the polarization degree of radiation from a plane-parallel slab in the diffusion approximation for a pure scattering atmosphere.
+You may use either numerical solution of the coupled system of radiative transfer
+equations as defined in Pomraning (1973) or Monte Carlo method.
+
+1) Library requirements
+
+The following third-party libraries are required for compilation
+
+1) Gnu Scientific Library
+2) CFITIO library
+3) MPI library
+
+
+As these packages may be installed in non-default folders, it is necessary
+to define some enviroment variables pointing to library headers and shared object.
+These enviroment variables, read by the Makefile at compilation step, are:
+
+LIB_CFITSIO
+HEADERS_CFITSIO
+
+LIB_GSL
+HEADERS_GSL
+
+LIB_MPI
+HEADERS_MPI
+
+2) Creation of the executable
+
+From the local folder where the Makefile is located just type
+
+$make
+
+If everything worked well, you will see the executable which is called "start_iteration"
+To remove object files
+
+$make clean
+
+3) Running the code
+
+The executable is launched from the shell prompt with a list of parameters, with synopsis like
+
+$start_iteration par1=value par2=value etc etc
+
+======================================================================
+3.1) Numerical solution
+======================================================================
+
+The parameter selecting this choice of the algorithm is "method=ode"
+
+Three locations of the seed photons are possible, and are defined by the parameter "seed=N" where
+
+N=1 photons at the center of the slab
+N=2 photons uniformly distributed across the slab
+N=3 photons at the base of the slab
+
+The Thomson optical depth of the slab is defined with parameter "disktau=value"
+
+The number of iteration is defined by the parameter "kiter=N", where N must be an integer greater than zero.
+To obtain reliable results it is important to consider the optical depth of the slab.
+For instance, if disktau=1 then kiter=20 could be enough, on the other hand for e.g. tau=5, kiter=100 at least
+is recommended.
+
+At the end of the simulation the following files are produced
+
+intensity_ode_tauA_seedB.qdp
+diffusion_ode_tauA_seedB.qdp
+
+where A and B states for the value defined for parameters disktau and seed.
+The file format for intensity_ode_tauA_seedB.qdp is the following:
+
+Col1: u, the cosine with respect to the slab normal
+Col2: Pdeg, namely Ir(u)-Il(u)/Ir(u)+Il(u)/
+Col3: I(u)/I(u=1), specific intensity normalized at its value for u=1
+
+Output file diffusion_ode_tauA_seedB.qdp has two columns instead:
+
+Col1: k, namely number of scattering
+Col2: F(k), namely the total flux emerging from the top of the slab for each order of scattering
+
+
+======================================================================
+3.2) Monte Carlo method
+======================================================================
+
+Here you must select from command line "method=mc"
+While the slab optical depth must be defined as well, for Monte Carlo method the following parameters are required:
+
+polardisk=y|n
+
+which defines if seed photons are initilly polarized (default is no).
+
+albedobase=A
+
+where A must be defined between 0 and 1. The code assumes an energy-independent albedo at the base of the slab.
+To perform tight comparison of the results with numerical method, you must set A=0.
+
+nph=N
+
+the number of photons for each core.
+There are two additional optional parameters to be passed at the prompt, namely the BB seed photons and electron
+temperature, respectively.
+They can be passed simply as
+
+ktbb=value kte=value
+
+and their default values is 1 for both cases.
+
+To run simulation with single core, you must define at shell prompt
+
+$start_simulation (plus parameter list)
+
+Tu run multi-core simulation use instead
+
+$mpirun -n Ncores start_simulation (plus parameter list)
+
+The output files are the following:
+
+a) spectrum_mc_tauN_seedM_AQ.qdp
+
+Energy spectra at the top of the slab with format:
+
+Col1: E
+Col2: Delta_E
+Col3: I(E)
+Col4-Col25: angle-dependent I(E)
+
+b) integralpolar_mc_tauN_seedM_AQ.qdp
+
+Total, energy-integrated, polarimetric results with format:
+
+Col1: u
+Col2: Pdeg (u)
+Col3: Q(u)
+Col4: U(u)
+Col5: Itot(u)/I(1)*1/u
+
+c) energypolar_mc_tauN_seedM_AQ.qdp
+
+Energy-dependent polarimetric results with column description defined in the file
+
+d) diffusion_mc_tauN_seedM_AQ.qdp
+
+Col1: number of scatterings k
+Col2: P(k) normalized photon distribution over number of scatterings
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/main_program.c b/main_program.c
index 41549a5..20705e8 100644
--- a/main_program.c
+++ b/main_program.c
@@ -50,7 +50,7 @@ int main(int argc, char* argv[])
"distribution at the top of a plane-parallel scattering atmosphere\n");
printf("Synopsis\n");
- printf("start_interation disktau=disktau seed=[1,2,3] method=ode|st85|mc\n\n");
+ printf("start_interation disktau=disktau seed=[1,2,3] method=ode|mc\n\n");
printf("Shared parameters\n\n");
printf("disktau=vaue: slab optical depth\n\n");
@@ -63,10 +63,6 @@ int main(int argc, char* argv[])
printf("Parameters for iteration algorithms\n");
printf("method=ode : solution of the ODE system for Il and Ir\n");
- printf(
- "method=st85 : use ST85 algorithm (valid only for seed photon at the center or uniformly "
- "distributed)\n");
-
printf("kiter=N: number of iterations\n\n");
printf("Parameters for MC simulation\n");
@@ -169,7 +165,7 @@ int main(int argc, char* argv[])
if (method == NULL || method == 0x0)
{
- printf("Please select the method to compute results with method=ode|st85|mc\n");
+ printf("Please select the method to compute results with method=ode|mc\n");
return 1;
}
@@ -282,18 +278,14 @@ int main(int argc, char* argv[])
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);
-
- simul_sync();
-
+ status = slab_mc(nph, seed);
- // check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);
}
else
{
printf(
- "\nPlease select method=ode (solve ODE), method=st85 (ST85) or method=mc (Montecarlo)\n");
+ "\nPlease select method=ode (solve ODE) or method=mc (Montecarlo)\n");
return (1);
}
diff --git a/make2tar.py b/make2tar.py
old mode 100644
new mode 100755
index 46bc6ed..accbb7f
--- a/make2tar.py
+++ b/make2tar.py
@@ -43,14 +43,14 @@ def Create_Makefile(FileName=None):
file.write("#Trova tutti i file sorgenti in tutte le directory\n")
file.write("SOURCES := $(foreach dir, $(dirs),$(call find_sources, $(dir)))\n\n")
- file.write("# Genera i file oggetto corrispondenti ai file sorgenti\n")
+ file.write("#Genera i file oggetto corrispondenti ai file sorgenti\n")
file.write("OBJECTS = $(patsubst %.c, %.o, $(SOURCES))\n\n")
- file.write("#Nome dell'eseguibile\n")
- file.write("EXECUTABLE = $(BINDIR)/my_program\n\n\n")
+ file.write("#Executable name\n")
+ file.write("EXECUTABLE = $(BINDIR)/start_iteration\n\n\n")
- file.write("# Regola di default: compila l'eseguibile\n")
+ file.write("#Default rule: compile executable\n")
file.write("all: $(EXECUTABLE)\n\n")
file.write("$(EXECUTABLE): $(OBJECTS)\n")
diff --git a/ode_example.c b/ode_example.c
deleted file mode 100644
index 273c820..0000000
--- a/ode_example.c
+++ /dev/null
@@ -1,27 +0,0 @@
-int main(void)
-{
- double mu = 10;
- gsl_odeiv2_system sys = {func, jac, 2, &mu};
-
- gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, 1e-3, 1e-8, 1e-8);
-
- double t = 0.0;
- double y[2] = {1.0, 0.0};
- int i, s;
-
- for (i = 0; i < 100; i++)
- {
- s = gsl_odeiv2_driver_apply_fixed_step(d, &t, 1e-3, 1000, y);
-
- if (s != GSL_SUCCESS)
- {
- printf("error: driver returned %d\n", s);
- break;
- }
-
- printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
- }
-
- gsl_odeiv2_driver_free(d);
- return s;
-}
--
GitLab