Wrote more detailed README.txt

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

        "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");

    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");

  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;

  }

    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);

  }

        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")

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;

}