#include "iteration_functions.h" /* Perform integration of the function I(tau, u) over the range u=[0-1] using Gauss-Legendre quadrature */ /*========================================================*/ /*Equation (4.206)*/ /*Compute the integral*/ /* [2*(1-u^2)*(1-u'^2) + u^2 u'^2] I_l(u) du' */ /*========================================================*/ double legendre_integration_A(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u) { int status; int ii; double integ_value; double angular_term; double Ik_l; double tolerance; integ_value = 0; Ik_l = 0; tolerance=1e-7; for (ii = 0; ii < Nstep_mu; ii++) { angular_term = scattering_matrix(1, 1, u, u_prime[ii]); /*===================================================================*/ /*Evaluate by 2D-interpolation the function I^{k-1}_l (tau, u')*/ /*====================================================================*/ status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l); //printf("u %lf tau %lf Ik_l %5.4e\n", u_prime[ii], tau, Ik_l); if (Ik_l < -tolerance) { //("Warning: interpolated Ik_l < 0 for tau=%5.4f u_prime=%5.4f\n", tau, u_prime[ii]); //exit(1); } if (status) { printf("Error at line %d\n", __LINE__); printf("File %s\n", __FILE__); printf("tau=%5.4f u=%5.4f I^{k}(l)=%5.4e\n", tau, u_prime[ii], Ik_l); exit(1); } // printf("ii=%d w %5.4e\n", ii, weights_u[ii]); integ_value = integ_value + angular_term * weights_u[ii] * Ik_l; } if (isnan(integ_value)) { printf("Error: values is NaN\n"); printf("Line %d\n", __LINE__); exit(1); } // printf("Integral (A) %lf\n", integ_value); return integ_value; } /*========================================================*/ /*Equation (4.206)*/ /*Compute the integral*/ /* u^2 Ir(u)*/ /*========================================================*/ double legendre_integration_B(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u) { int status; int ii; double integ_value; double angular_term; double Ik_r; integ_value = 0; Ik_r = 0; for (ii = 0; ii < Nstep_mu; ii++) { angular_term = scattering_matrix(1, 2, u, u_prime[ii]); /*===================================================================*/ /*Evaluate by 2D-interpolation the function I^{k-1}_r (tau, u')*/ /*====================================================================*/ status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r); if (status) { printf("Error at line %d\n", __LINE__); printf("File %s\n", __FILE__); exit(1); } integ_value = integ_value + angular_term * weights_u[ii] * Ik_r; } if (isnan(integ_value)) { printf("Error: values is NaN\n"); printf("Line %d\n", __LINE__); exit(1); } return integ_value; } /*========================================================*/ /*Equation (4.207)*/ /*Compute the integral*/ /* u'^2 Il(u)*/ /*========================================================*/ double legendre_integration_C(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u) { int status; int ii; double integ_value; double angular_term; double Ik_l; integ_value = 0; Ik_l = 0; for (ii = 0; ii < Nstep_mu; ii++) { angular_term = scattering_matrix(2, 1, u, u_prime[ii]); /*===================================================================*/ /*Evaluate by 2D-interpolation the function I^{k-1}_l (tau, u')*/ /*====================================================================*/ status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_l); if (status) { printf("Error at line %d\n", __LINE__); printf("File %s\n", __FILE__); exit(1); } integ_value = integ_value + angular_term * weights_u[ii] * Ik_l; } if (isnan(integ_value)) { printf("Error: values is NaN\n"); printf("Line %d\n", __LINE__); exit(1); } return integ_value; } /*========================================================*/ /*Equation (4.207)*/ /*Compute the integral*/ /* Ir(u)*/ /*========================================================*/ double legendre_integration_D(gsl_spline2d* Spline_I2d, double tau, double u, double* u_prime, double* weights_u) { int status; int ii; double integ_value; double angular_term; double Ik_r; integ_value = 0; Ik_r = 0; for (ii = 0; ii < Nstep_mu; ii++) { angular_term = scattering_matrix(2, 2, u, u_prime[ii]); /*===================================================================*/ /*Evaluate by 2D-interpolation the function I^{k-1}_r (tau, u')*/ /*====================================================================*/ status = gsl_spline2d_eval_extrap_e(Spline_I2d, tau, u_prime[ii], xacc_2d, yacc_2d, &Ik_r); if (status) { printf("Error at line %d\n", __LINE__); printf("File %s\n", __FILE__); exit(1); } integ_value = integ_value + angular_term * weights_u[ii] * Ik_r; } if (isnan(integ_value)) { printf("Error: values is NaN\n"); printf("Line %d\n", __LINE__); exit(1); } return integ_value; } /*================================================================================*/ double scattering_matrix(int row, int col, double u, double u_prime) { double value; if (row == 1 && col == 1) { value = 2 * (1 - u * u) * (1 - u_prime * u_prime) + u * u * u_prime * u_prime; } else if (row == 1 && col == 2) { value = u * u; } else if (row == 2 && col == 1) { value = u_prime * u_prime; } else if (row == 2 && col == 2) { value = 1; } else { printf("Uncorrect row or column value for scattering matrix\n"); printf("row=%d col=%d\n", row, col); exit(1); } return value; }