From f82f66781b328da03a514f8fca01b65d774381e9 Mon Sep 17 00:00:00 2001 From: Ruben Farinelli <ruben.farinelli@inaf.it> Date: Wed, 27 Sep 2023 09:30:09 +0200 Subject: [PATCH] Moved the header file iteration_functions.h to the local folder --- iteration_functions.h | 242 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 242 insertions(+) create mode 100644 iteration_functions.h diff --git a/iteration_functions.h b/iteration_functions.h new file mode 100644 index 0000000..a6ea2c8 --- /dev/null +++ b/iteration_functions.h @@ -0,0 +1,242 @@ +#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 -- GitLab