Loading iteration_functions.h 0 → 100644 +242 −0 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 Loading
iteration_functions.h 0 → 100644 +242 −0 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