Commit 0c63145f authored by Ruben Farinelli's avatar Ruben Farinelli
Browse files

Latest working version of the code fully consistent with the numerical iteration method

parent 08250fe8

compute_stokes.c

+2 −6
Original line number Diff line number Diff line
@@ -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);

    }

  }

iteration_functions.h

deleted100644 → 0
+0 −242
Original line number Diff line number Diff line 
#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

mc_slab.c

+46 −58
Original line number Diff line number Diff line
@@ -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,6 +145,7 @@ int slab_mc(int nphot, int seed) 
    setup_gsl_objects(radiation_field, MC_SLAB);

  }





  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,30 +265,27 @@ 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];



        Magfield_lab[0] = Magfield_lab_new[0];

        Magfield_lab[1] = Magfield_lab_new[1];

        Magfield_lab[2] = Magfield_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];

        }



      }

      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,7 +502,8 @@ 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)
@@ -523,7 +511,7 @@ void comptonization(double E_lab, 
    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);

solve_rte.c

+1 −1
Original line number Diff line number Diff line
@@ -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);

  }