Set additional seed photon distribution with exponential law (c86fa72a) · Commits · polarization / slab

README.txt

+74 −54

Original line number	Original line	Diff line number	Diff line
	The code computes the polarization degree of radiation from a plane-parallel slab in the diffusion approximation for a pure scattering atmosphere.		The code computes the spectropolarimetric properties of radiation undergoing electron scattering in a plane-parallel slab.
	You may use either numerical solution of the coupled system of radiative transfer		One may use either a numerical solution of the coupled system of radiative transfer equations, as defined in Pomraning (1973), or a Monte Carlo method. The latter approach is better suited for computations at high electron temperatures, for which the Klein–Nishina cross-section provides more accurate results.
	equations as defined in Pomraning (1973) or Monte Carlo method.

			======================================================================
	1) Library requirements		1) Library requirements
			======================================================================

	The following third-party libraries are required for compilation		The following third-party libraries are required for compilation

	@@ -10,10 +11,12 @@ The following third-party libraries are required for compilation
	2) CFITIO library		2) CFITIO library
	3) MPI library		3) MPI library

			As these packages may be installed in non-default folders, it is necessary to define some environment variables. These variables tell the compiler and linker where to find the library files and their headers needed to build the code.
			Library headers are files ending in .h that contain definitions needed to use the library in the code.

	As these packages may be installed in non-default folders, it is necessary		Shared objects (.so) or static libraries (.a) are files containing the compiled library code used by the program at runtime (shared) or at compile-time (static).
	to define some enviroment variables pointing to library headers and shared object.
	These enviroment variables, read by the Makefile at compilation step, are:		These environment variables are read by the Makefile at the compilation step. The required variables are:

	LIB_CFITSIO		LIB_CFITSIO
	HEADERS_CFITSIO		HEADERS_CFITSIO
	@@ -24,7 +27,19 @@ HEADERS_GSL
	LIB_MPI		LIB_MPI
	HEADERS_MPI		HEADERS_MPI

			Each of them must point to the associated library and to the folder with header files, respectively.

			While you can set these variables in each terminal session, it is recommended to define them once and for all in your shell initialization files:

			For Bash: add the export commands to your ~/.bashrc or ~/.profile

			For C shell (csh/tcsh): add the setenv commands to your ~/.cshrc or ~/.tcshrc

			This way, the variables are automatically set every time you open a new terminal, and you won’t need to define them manually each time.

			======================================================================
	2) Creation of the executable		2) Creation of the executable
			======================================================================

	From the local folder where the Makefile is located just type		From the local folder where the Makefile is located just type

	@@ -42,16 +57,17 @@ The executable is launched from the shell prompt with a list of parameters, with
	$start_iteration par1=value par2=value etc etc		$start_iteration par1=value par2=value etc etc

	======================================================================		======================================================================
	3.1) Numerical solution		3) Iterative method
	======================================================================		======================================================================

	The parameter selecting this choice of the algorithm is "method=ode"		The parameter selecting this choice of the iterative algorithm is "method=ode"

	Three locations of the seed photons are possible, and are defined by the parameter "seed=N" where		Four 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=1 photons at the center
	N=2 photons uniformly distributed across the slab		N=2 photons uniformly distributed
	N=3 photons at the base of the slab		N=3 photons at the base
			N=4 photons with exponential distribution

	The Thomson optical depth of the slab is defined with parameter "disktau=value"		The Thomson optical depth of the slab is defined with parameter "disktau=value"

	@@ -79,71 +95,75 @@ Col2: F(k), namely the total flux emerging from the top of the slab for each ord


	======================================================================		======================================================================
	3.2) Monte Carlo method		4) Monte Carlo method
	======================================================================		======================================================================

	Here you must select from command line "method=mc"		To run the code using the Monte Carlo method, select at at the shell prompt:
	While the slab optical depth must be defined as well, for Monte Carlo method the following parameters are required:
			method=mc

			The slab optical depth must also be specified. The following Monte Carlo–specific parameters are required.

			Parameters

	polardisk =y\|n		polardisk =y\|n

	which defines if seed photons are initilly polarized (default is no).		Specifies whether the seed photons are initially polarized. This is valid only for the
			configuration with seed photons at the base of the slab (seed=3).

	albedobase = A		albedobase = A

	where A must be defined between 0 and 1. The code assumes an energy-independent albedo at the base of the slab.		Energy-independent albedo at the base of the slab.
	To perform tight comparison of the results with numerical method, you must set A=0.		A must be defined between 0 and 1.
			To perform a direct comparison with the iterative method, set A=0.

	nph = N		nph = N

	the number of photons for each core.		Number of photons per MPI process.
	There are two additional optional parameters to be passed at the prompt, namely the BB seed photons and electron
	temperature, respectively.		The following optional parameters control the BB seed photon and electron temperature, both expressed in keV:
	They can be passed simply as

	ktbb=value kte=value		ktbb=value kte=value

	and their default values is 1 for both cases.		Default values are ktbb=1 and kte=1.

	To run simulation with single core, you must define at shell prompt		By default, the results are computed for ten viewing angles. To change this value, use the parameter

	$start_simulation (plus parameter list)		nstepangles=value

	Tu run multi-core simulation use instead		Running the simulation

	$mpirun -n Ncores start_simulation (plus parameter list)		To run a single-core simulation, use:

	The output files are the following:		start_simulation [parameter list]

	a) spectrum_mc_tauN_seedM_AQ.qdp		To run a multi-core (MPI) simulation, use:

	Energy spectra at the top of the slab with format:		mpirun -n Ncores start_simulation [parameter list]

	Col1: E		Output files
	Col2: Delta_E
	Col3: I(E)
	Col4-Col25: angle-dependent I(E)

	b) integralpolar_mc_tauN_seedM_AQ.qdp		The code produces the following output files:

	Total, energy-integrated, polarimetric results with format:		spectrum_mc_tauN_seedM_AQ.qdp

	Col1: u		Energy spectrum at the top of the slab.
	Col2: Pdeg (u)
	Col3: Q(u)
	Col4: U(u)
	Col5: Itot(u)/I(1)*1/u

	c) energypolar_mc_tauN_seedM_AQ.qdp		Col 1: photon energy E
			Col 2: energy bin width Delta E
			Col 3: total intensity I(E)
			Col 4–nstepangles : angle-dependent intensity I(E), with u=cos(i) progressively decreasing

	Energy-dependent polarimetric results with column description defined in the file		integralpolar_mc_tauN_seedM_AQ.qdp

	d) diffusion_mc_tauN_seedM_AQ.qdp		Energy-integrated spectropolarimetric quantities, with more detailed column description provided at the beginning of the file.

	Col1: number of scatterings k		diffusion_mc_tauN_seedM_AQ.qdp
	Col2: P(k) normalized photon distribution over number of scatterings

			Photon diffusion statistics.

			Col1: number of scatterings k
			Col2: normalized photon distribution function P(k)

main_program.c

+4 −4

Original line number	Original line	Diff line number	Diff line
	#include "globals.h"		#include "globals.h"
	#include "variables.h"		#include "variables.h"
	#include <mpi.h>
	#include <stdio.h>
	#include <common_header.h>		#include <common_header.h>
	#include <iteration_functions.h>		#include <iteration_functions.h>
			#include <mpi.h>
			#include <stdio.h>

	char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method);		char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method);
	char* set_mc_filename(const char* root, double ktbb, double kte, double tau, int seed, double albedo, const char* method);		char* set_mc_filename(const char* root, double ktbb, double kte, double tau, int seed, double albedo, const char* method);

mc_slab.c

+67 −51

Original line number	Original line	Diff line number	Diff line
	@@ -48,7 +48,7 @@ gsl_root_fsolver* s_maxw;

	/==================================================================/		/==================================================================/

	int slab_mc(int nphot, int seed)		int slab_mc(int nphot, int seed_distribution)
	{		{
	int ii;		int ii;
	int jj;		int jj;
	@@ -81,6 +81,7 @@ int slab_mc(int nphot, int seed)
	double csi;		double csi;
	double Pdeg;		double Pdeg;
	double kn_corr;		double kn_corr;
			double tau_sampled;

	int naxis1 = 0;		int naxis1 = 0;
	int naxis2 = 0;		int naxis2 = 0;
	@@ -102,29 +103,27 @@ int slab_mc(int nphot, int seed)

	/===================================================================================/		/===================================================================================/

			if (getenv("MC_SLAB") == NULL)
	if (getenv("MC_SLAB") == NULL) {		{
	printf("Please set the enviroment variable MC_SLAB pointing to the folder where it is defined the code executabe\n");		printf(
			"Please set the enviroment variable MC_SLAB pointing to the folder where it is defined "
			"the code executabe\n");
	exit(EXIT_FAILURE);		exit(EXIT_FAILURE);
	}		}

	ptr_gb = read_array(stringconcat(getenv("MC_SLAB"), "/temp_random.txt"));		ptr_gb = read_array(stringconcat(getenv("MC_SLAB"), "/temp_random.txt"));


	matrix_gb = stringconcat(getenv("MC_SLAB"), "/gammabeta_lab.fits");		matrix_gb = stringconcat(getenv("MC_SLAB"), "/gammabeta_lab.fits");
	array_gb = read_matrix(matrix_gb);		array_gb = read_matrix(matrix_gb);

	/===================================================================================/		/===================================================================================/

	char* knfile_maxw;		char* knfile_maxw;
	knfile_maxw =		knfile_maxw = stringconcat(getenv("MC_SLAB"), "/S00_lab.fits");
	stringconcat(getenv("MC_SLAB"), "/S00_lab.fits");


	matrix_kn = read_matrix(knfile_maxw);		matrix_kn = read_matrix(knfile_maxw);

	ptr_kn = read_array(		ptr_kn = read_array(stringconcat(getenv("MC_SLAB"), "/energy_temp.txt"));
	stringconcat(getenv("MC_SLAB"), "/energy_temp.txt"));

	/===================================================================================/		/===================================================================================/

	@@ -258,6 +257,11 @@ int slab_mc(int nphot, int seed)
	printf("Photon %d\n", ii);		printf("Photon %d\n", ii);
	}		}

			if (!emin)
			{
			emin = 0.1;
			}

	n_sc = 0;		n_sc = 0;
	E_lab = bb_spec(ktbb);		E_lab = bb_spec(ktbb);
	E_new_lab = 0;		E_new_lab = 0;
	@@ -265,21 +269,17 @@ int slab_mc(int nphot, int seed)
	pos[0] = 0;		pos[0] = 0;
	pos[1] = 0;		pos[1] = 0;

	if (seed == 1)		if (seed_distribution == SEED_CENTER)
	{		{
	pos[2] = Hphot / 2;		pos[2] = Hphot / 2;
	cos_theta = 2 * clock_random() - 1;		cos_theta = 2 * clock_random() - 1;
	}		}
	else if (seed == 2)		else if (seed_distribution == SEED_UNIFORM)
	{
	if (!emin)
	{		{
	emin = 0.1;
	}
	pos[2] = clock_random() * Hphot;		pos[2] = clock_random() * Hphot;
	cos_theta = 2 * clock_random() - 1;		cos_theta = 2 * clock_random() - 1;
	}		}
	else		else if (seed_distribution == SEED_BASE)
	{		{
	pos[2] = 0;		pos[2] = 0;

	@@ -301,6 +301,21 @@ int slab_mc(int nphot, int seed)
	cos_theta = clock_random();		cos_theta = clock_random();
	}		}
	}		}
			else if (seed_distribution == SEED_EXP)
			{
			csi = clock_random();

			tau_sampled = log(1 / (1 + csi * (-1 + exp(-disktau))));

			pos[2] = tau_sampled / disktau * Hphot;

			cos_theta = 2 * clock_random() - 1;
			}
			else
			{
			printf("Error: values for seed photon can be 1,2,3,4\n");
			exit(1);
			}

	sin_theta = sqrt(1 - cos_theta * cos_theta);		sin_theta = sqrt(1 - cos_theta * cos_theta);
	phi = 2 * PI * clock_random();		phi = 2 * PI * clock_random();
	@@ -653,9 +668,6 @@ void comptonization(double E_lab,
	printf("Error in photon energy in lab frame after scattering: value %5.4f\n", *E_new_lab);		printf("Error in photon energy in lab frame after scattering: value %5.4f\n", *E_new_lab);
	exit(1);		exit(1);
	}		}



	}		}

	/======================================================================/		/======================================================================/
	@@ -685,17 +697,22 @@ void write_kn_histo(const char* filename)
	{		{
	// Trova il minimo		// Trova il minimo
	int min_val = histo_kn[0];		int min_val = histo_kn[0];
	for (int i = 1; i < NBIN_KN; i++) {		for (int i = 1; i < NBIN_KN; i++)
	if (histo_kn[i] < min_val) min_val = histo_kn[i];		{
			if (histo_kn[i] < min_val)
			min_val = histo_kn[i];
	}		}

	if (min_val == 0) {		if (min_val == 0)
	fprintf(stderr, "Attenzione: il valore minimo dell'istogramma è 0, normalizzazione impossibile!\n");		{
			fprintf(stderr,
			"Attenzione: il valore minimo dell'istogramma è 0, normalizzazione impossibile!\n");
	exit(EXIT_FAILURE);		exit(EXIT_FAILURE);
	}		}

	FILE* fp = fopen(filename, "w");		FILE* fp = fopen(filename, "w");
	if (!fp) {		if (!fp)
			{
	perror("Errore apertura file");		perror("Errore apertura file");
	exit(EXIT_FAILURE);		exit(EXIT_FAILURE);
	}		}
	@@ -714,4 +731,3 @@ void write_kn_histo(const char* filename)

	fclose(fp);		fclose(fp);
	}		}

solve_rte.c

+1 −1

Original line number	Original line	Diff line number	Diff line
	@@ -131,7 +131,7 @@ int solve_rte(int k_iter, int seed_distribution, int Nstep_tau)
	else		else
	{		{
	printf(		printf(
	"Error: please select seed=4 (photons at the base) or seed=2 (uniform distribution)\n");		"Error: values for seed photon can be 1,2,3,4\n");
	exit(1);		exit(1);
	}		}
	}		}