legendre_integration.c

#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;
}