diff --git a/compute_stokes.c b/compute_stokes.c
index 7dbc1984b878614f6c50dec7bb19f40e550f7529..061a6e634d2850bc23b66d2fa2d598b2343d36b3 100644
--- a/compute_stokes.c
+++ b/compute_stokes.c
@@ -36,17 +36,13 @@ void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_para
{
if ((cos_theta <= array_ubounds[jj]) && (cos_theta >= array_ubounds[jj + 1]))
{
- /* if (jj==nstepangles-1)
- {
- printf("Angles [%5.4f - %5.4f]\n", array_ubounds[jj], array_ubounds[jj+1]);
- }
- */
detector_axis(k_lab, Qp_cart, Qm_cart, Up_cart, Um_cart);
Qs = pow(dot_prod(Qp_cart, polvect_cart), 2) - pow(dot_prod(Qm_cart, polvect_cart), 2);
Us = pow(dot_prod(Up_cart, polvect_cart), 2) - pow(dot_prod(Um_cart, polvect_cart), 2);
+
if (isnan(Qs))
{
printf("Error: Qs is Nan\n");
@@ -75,7 +71,7 @@ void compute_stokes(double ene, double* k_lab, double* polvect_cart, stokes_para
if (FlagFound == FALSE)
{
double zio = 1;
- // printf("Warning: no correspondence for Ene=%5.4f k_z=%5.4f\n", ene, k_lab[2]);
+ printf("Warning: no correspondence for Ene=%5.4f k_z=%5.4f\n", ene, k_lab[2]);
// exit(1);
}
}
diff --git a/iteration_functions.h b/iteration_functions.h
deleted file mode 100644
index a6ea2c8aedf50e1cd8dfcc607756699e057bfcec..0000000000000000000000000000000000000000
--- a/iteration_functions.h
+++ /dev/null
@@ -1,242 +0,0 @@
-#include "local_header.h"
-#include "globals.h"
-#include "common_header.h"
-#include "gsl/gsl_odeiv2.h"
-#include "define_values.h"
-
-
-#ifndef ITERATION_FUNCIONS_H_
-#define ITERATION_FUNCIONS_H_
-
-#define SEED_CENTER 1
-#define SEED_UNIFORM 2
-#define SEED_BASE 3
-#define SEED_DELTA 4
-#define SEED_EXP 5
-
-#define RTE 1
-#define MC_SLAB 2
-
-#define U_POSITIVE 1
-#define U_NEGATIVE -1
-
-extern FILE* dat_spectrum;
-extern FILE* dat_diffusion;
-extern FILE* dat_integral_polar;
-extern FILE* dat_energy_polar;
-
-extern int seed_location;
-extern int Nstep_mu;
-extern int Nstep_tau;
-
-extern double tau0;
-extern double tau_c;
-
-extern gsl_spline2d *Spline_I2d_l;
-extern gsl_spline2d *Spline_I2d_r;
-
-
-extern gsl_interp_accel* xacc_2d;
-extern gsl_interp_accel* yacc_2d;
-
-extern gsl_spline* Spline_sample_Iu;
-extern gsl_spline* Spline_eval_limbdark;
-extern gsl_spline* Spline_eval_pdeg;
-extern gsl_spline *Spline_tau_Ak;
-extern gsl_spline *Spline_tau_Bk;
-extern gsl_spline *Spline_tau_Ck;
-
-extern gsl_interp_accel* xacc_1d;
-extern gsl_interp_accel* yacc_1d;
-
-extern char* polarization_file;
-extern char* diffusion_file;
-extern char* integral;
-extern char* polarinput;
-
-typedef struct {
-
- double* array_mu;
- double* array_tau;
- double** matrix;
-
-
-
-} Ikl_intensity;
-
-typedef struct {
-
- double* array_mu;
- double* array_tau;
- double** matrix;
-
-} Ikr_intensity;
-
-
-typedef struct {
-
- int idx_tau;
- double * weights_u;
- double u;
- double* array_u;
- double* array_tau;
- double** Il_matrix_upstream;
- double** Ir_matrix_upstream;
-
- double** Il_matrix_downstream;
- double** Ir_matrix_downstream;
-
- gsl_spline2d *Spline_I2d_l_upstream;
- gsl_spline2d *Spline_I2d_l_downstream;
- gsl_spline2d *Spline_I2d_r_upstream;
- gsl_spline2d *Spline_I2d_r_downstream;
-
-
-} Iklr_intensity;
-
-
-typedef struct {
-
- int k;
- int seed_location;
- double tau_s;
- double tau0;
- double tau;
- double* array_mu;
- double* array_tau;
- double** matrix;
-
-} AngularFunctionParams;
-
-
-typedef struct {
-
- int nlines;
- double u_min;
- double u_max;
- double* u;
- double* pdeg;
- double* Q;
- double* U;
- double* limbdark;
- gsl_spline* Spline;
- gsl_interp_accel* acc;
-
-
-} radiation_field_t;
-
-typedef struct {
- double *array_rand;
- double *array_u;
- double *array_integ;
- double total_integral;
- int nlines;
- gsl_spline *spline;
- gsl_interp_accel *acc1;
- } cumulative_angdist_t;
-
-
-
-/*========================================================================*/
-/*Functions prototypes*/
-/*========================================================================*/
-
-
-double f_interp(double u, double tau);
-double I0_center(double tau0, double tau, double u);
-double I0_delta(double tau, double tau_s, double mu);
-double I0_uniform(double tau0, double tau, double mu);
-double I0_base(double tau0, double tau, double mu, int sign_mu);
-double I0_exp(double tau, double tau0, double mu);
-
-double theta_funct(double tau, double tau_c);
-void make_linarray(double a, double b, int nstep, double* my_array);
-double** dmatrix(long nrl, long nrh, long ncl, long nch);
-
-double Ak_integral (double x, void * params);
-double Bk_integral (double x, void * params);
-double Ck_integral (double x, void * params);
-
-double Il_integral (double x, void * params);
-double Ir_integral (double x, void * params);
-
-void Compute_Ak_array(AngularFunctionParams* ptr_data,
- double** Il_funct, double* za, double*array_tau, double*array_mu, int k,
- gsl_integration_workspace * w, double* Ak_array);
-
-void Compute_Bk_array(AngularFunctionParams* ptr_data,
- double** Il_funct, double* za, double*array_tau, double*array_mu, int k,
- gsl_integration_workspace * w, double* Bk_array);
-
-void Compute_Ck_array(AngularFunctionParams* ptr_data,
- double** Ir_funct, double* za, double*array_tau, double*array_mu, int k,
- gsl_integration_workspace * w, double* Ck_array);
-
-double round_to_digits(double value, int digits);
-
-void slab_polarization(int kiter, int seed_location);
-
-void read_resultsfile(char* filename, radiation_field_t* radiation_field, int Flag);
-void compute_cumulative_angdist(radiation_field_t* radiation_field, cumulative_angdist_t* cumulative_angdist);
-
-void check_results(gsl_spline* Spline_sample_Iu, gsl_spline* Spline_eval_pdeg, gsl_spline* Spline_eval_limbdark);
-
-/* ========================================================================= */
-/* Functions used to integrate the system of two ODE */
-/* ========================================================================= */
-
-
-extern gsl_spline* Spline_I2d_l_upstream;
-extern gsl_spline* Spline_I2d_l_downstream;
-extern gsl_spline* Spline_I2d_r_upstream;
-extern gsl_spline* Spline_I2d_l_downstream;
-
-double legendre_integration_A(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
-double legendre_integration_B(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
-double legendre_integration_C(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
-double legendre_integration_D(gsl_spline2d *Spline_I2d, double tau, double u, double* u_prime, double* weights_u);
-
-
-double scattering_matrix (int row, int col, double u, double u_prime);
-
-int solve_rte(int k_iter, int seed_distribution, int Nstep_tau);
-int RTE_Equations(double s, const double y[], double f[], void* params);
-
-
-void compute_results(int method,
- int k_iter,
- int Nstep_tau,
- int Nstep_mu,
- double* array_mu,
- double* weights_u,
- Iklr_intensity* ptr_Iklr,
- Ikl_intensity* ptr_Ikl,
- Ikr_intensity* ptr_Ikr);
-
-
-
-
-//==============================================================
-//Routines related to local MC
-//==============================================================
-
-int slab_mc(int nphot, int seed);
-void compute_stokes(double ene, double* k_lab, double* polvect_cart,
- stokes_parameters* ptr_stokes, double* array_ubounds, int flag);
-
-
-void collect_photons(double ene,
- double* k,
- double* array_ene,
- double** array_counts,
- int flag);
-
-void setup_gsl_objects(radiation_field_t* radiation_field, int Flag);
-void electricfield_fromfile (double* k_phot, double* elect_field_lab, double* mag_field_lab, double *Pdeg);
-
-double subrel_maxwellian (double kte);
-
-void optimize_taustep(int k_iter, int seed_distribution);
-
-
-#endif
diff --git a/mc_slab.c b/mc_slab.c
index 99291417d3a22b459c7a07609c102e167ad0e306..a4936af9b6cc1130237880d0c50ae08320f6d564 100644
--- a/mc_slab.c
+++ b/mc_slab.c
@@ -12,12 +12,12 @@ void comptonization(double E_lab,
double* k_lab,
double* Efield_lab,
double* Magfield_lab,
- int* PolarFlag,
double* E_new_lab,
double* k_new_lab,
double* Efield_lab_new,
double* Magfield_lab_new);
+
double distance_bottom(double z0, double kz);
double distance_top(double z0, double kz, double H);
@@ -57,14 +57,15 @@ int slab_mc(int nphot, int seed)
{
int ii;
int jj;
- int PolarFlag;
- int FlagLocation;
int n_sc;
+ int n_zero_sc;
int nph_esc;
int Fval_max;
- int nsc_max;
+ int FlagLocation;
int photon_diffusion[NSC_MAX];
int status;
+
+
double E_lab;
double E_new_lab;
double k_lab[3];
@@ -74,8 +75,6 @@ int slab_mc(int nphot, int seed)
double Magfield_lab[3];
double Magfield_lab_new[3];
double pos[3];
- double Pdeg;
-
double cos_theta;
double sin_theta;
double phi;
@@ -87,6 +86,7 @@ int slab_mc(int nphot, int seed)
double eps;
double l_int;
double csi;
+ double Pdeg;
T_disk = gsl_root_fsolver_bisection;
s_disk = gsl_root_fsolver_alloc(T_disk);
@@ -145,7 +145,8 @@ int slab_mc(int nphot, int seed)
setup_gsl_objects(radiation_field, MC_SLAB);
}
- Pdeg = 0;
+
+ Pdeg=0;
nph_esc = 0;
cos_theta = 0;
@@ -165,6 +166,7 @@ int slab_mc(int nphot, int seed)
stokes_parameters struct_stokes_average[POLAR_NBOUNDS - 1];
allocate_stokes_pointer(struct_stokes_average);
+ n_zero_sc=0;
/*=============================================================*/
/*Here starts loop over photon */
/*=============================================================*/
@@ -182,7 +184,6 @@ int slab_mc(int nphot, int seed)
}
n_sc = 0;
-
E_lab = bb_spec(ktbb);
E_new_lab = 0;
@@ -229,27 +230,14 @@ int slab_mc(int nphot, int seed)
k_lab[1] = sin_theta * sin(phi);
k_lab[2] = cos_theta;
+
if (DiskPolarizationFlag == FROMFILE)
{
electricfield_fromfile(k_lab, Efield_lab, Magfield_lab, &Pdeg);
-
- if (Pdeg < 0)
- Pdeg = -Pdeg;
-
- if (clock_random() < Pdeg)
- {
- PolarFlag = 1;
- }
- else
- {
- PolarFlag = 0;
- }
}
-
else
{
unpolarized_seedphot(k_lab, Efield_lab, Magfield_lab);
- PolarFlag = 0;
}
/*====================================================================*/
@@ -277,31 +265,28 @@ int slab_mc(int nphot, int seed)
if (tau > eps)
{
l_int = eps / (NeDisk * SIGMA_T * r_units);
- // printf("Photon %d will interact at dist %5.4f\n", ii, l_int);
pos[0] = pos[0] + k_lab[0] * l_int;
pos[1] = pos[1] + k_lab[1] * l_int;
pos[2] = pos[2] + k_lab[2] * l_int;
- comptonization(E_lab, k_lab, Efield_lab, Magfield_lab, &PolarFlag, &E_new_lab, k_new_lab,
- Efield_lab_new, Magfield_lab_new);
+ //printf("Photon %d before scattering %d Pdeg = %5.4f\n", ii, n_sc, Pdeg);
- n_sc++;
+ comptonization(E_lab, k_lab, Efield_lab, Magfield_lab, &E_new_lab, k_new_lab,
+ Efield_lab_new, Magfield_lab_new);
- k_lab[0] = k_new_lab[0];
- k_lab[1] = k_new_lab[1];
- k_lab[2] = k_new_lab[2];
+ //printf("After scattering Pdeg = %5.4f\n\n", Pdeg);
+ n_sc++;
E_lab = E_new_lab;
- Efield_lab[0] = Efield_lab_new[0];
- Efield_lab[1] = Efield_lab_new[1];
- Efield_lab[2] = Efield_lab_new[2];
+ for (jj=0; jj<3; jj++)
+ {
+ k_lab[jj] = k_new_lab[jj];
+ Efield_lab[jj] = Efield_lab_new[jj];
+ Magfield_lab[jj] = Magfield_lab_new[jj];
+ }
- Magfield_lab[0] = Magfield_lab_new[0];
- Magfield_lab[1] = Magfield_lab_new[1];
- Magfield_lab[2] = Magfield_lab_new[2];
-
}
else
@@ -311,6 +296,8 @@ int slab_mc(int nphot, int seed)
{
// printf("Photon %d escaping after %d scattering\n", ii, n_sc);
+ if (n_sc==0) n_zero_sc++;
+
compute_stokes(E_lab, k_lab, Efield_lab, struct_stokes, array_ubounds, 0);
photon_diffusion[n_sc]++;
collect_photons(E_new_lab, k_lab, array_ene, array_counts, 0);
@@ -340,7 +327,6 @@ int slab_mc(int nphot, int seed)
pos[2] = 0;
unpolarized_seedphot(k_lab, Efield_lab, Magfield_lab);
- PolarFlag = 0;
E_lab = bb_spec(ktbb);
}
@@ -348,6 +334,7 @@ int slab_mc(int nphot, int seed)
else
{
FlagLocation = 0;
+ break;
}
} // Close the else with condition k_lab[2] < 0
@@ -381,20 +368,26 @@ int slab_mc(int nphot, int seed)
}
}
+ // Get the number of processes
+
+ MPI_Comm_size(MPI_COMM_WORLD, &nproc);
+
+
/*==================================================================*/
if (rank == 0)
{
compute_stokes(1, k_lab, Efield_lab, struct_stokes_average, array_ubounds, 1);
- printf("Percentage of escaping photons %4.3f\n", 1. * nph_esc / nphot);
-
+ printf("Percentage of escaping photons from the top layer %4.3f\n", 1. * nph_esc / nphot);
+ printf("Percentage of escaping non scattered photons from the top layer %4.3f\n", (double) n_zero_sc/ 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)
@@ -426,19 +419,22 @@ int slab_mc(int nphot, int seed)
} // End of function
+
void comptonization(double E_lab,
double* k_lab,
double* Efield_lab,
double* Magfield_lab,
- int* PolarFlag,
double* E_new_lab,
double* k_new_lab,
double* Efield_lab_new,
double* Magfield_lab_new)
{
+
int jj;
int status;
+ int PolarFlag;
+ double Pdeg;
double phik;
double E_erf;
double E_new_erf;
@@ -450,12 +446,14 @@ void comptonization(double E_lab,
double k_erf[3];
double k_new_erf[3];
double el_vel_vect[3];
-
double ElectField_new_erf[3];
double ElectField_old_erf[3];
-
double MagField_Erf[3];
+
+ PolarFlag=0;
+ Pdeg=0;
+
//================================================================== */
// Find the K' components in the ERF */
// If also polarization is computed, transform also the electric */
@@ -466,15 +464,10 @@ void comptonization(double E_lab,
versor_boost(k_lab, el_vel_vect, k_erf, 1);
- double beta_new = subrel_maxwellian(kte);
-
- // printf("pluto %lf %lf\n", beta_el, beta_new);
-
elmag_lorentztransform(Efield_lab, Magfield_lab, el_vel_vect, ElectField_old_erf, MagField_Erf, 1);
E_erf = energy_lorentz_boost(E_lab, el_vel_vect, k_lab, 1);
- // printf("beta %lf %lf \n", beta_el, E_erf);
csi = clock_random();
@@ -489,14 +482,8 @@ void comptonization(double E_lab,
E_new_erf = new_energy_erf(E_erf, costheta_erf);
- if (*PolarFlag == 1)
- {
- phik = sample_azim_angle(E_erf, E_new_erf, costheta_erf);
- }
- else
- {
- phik = 2 * PI * clock_random();
- }
+
+ phik = sample_azim_angle(E_erf, E_new_erf, costheta_erf, 1);
//================================================================== */
/*Find the components of the photon versor after scattering*/
@@ -515,15 +502,16 @@ void comptonization(double E_lab,
// and compute the new electric field vector
//================================================================== */
- status = depolarization_probability(E_erf, E_new_erf, costheta_erf, phik, PolarFlag);
+ PolarFlag=1;
+ status = depolarization_probability(E_erf, E_new_erf, costheta_erf, phik, &Pdeg, &PolarFlag);
-
+
if (status == 1)
{
exit(1);
}
- electfield_as_new(k_new_erf, ElectField_old_erf, ElectField_new_erf, *PolarFlag);
+ electfield_as_new(k_new_erf, ElectField_old_erf, ElectField_new_erf, PolarFlag);
vect_prod(k_new_erf, ElectField_new_erf, MagField_Erf);
diff --git a/solve_rte.c b/solve_rte.c
index 9d7ee65aec981b62b8edfe33168fb844997ac5b1..20dd9df6b35fd7e89b7baec30b9e1846ee4f1766 100644
--- a/solve_rte.c
+++ b/solve_rte.c
@@ -405,7 +405,7 @@ double I0_center_updown(double tau, double u)
{
double value;
- if (tau >= tau0)
+ if (tau > tau0)
{
value = 1 / 2. * 1 / u * exp(-(tau - tau0) / u);
}