[Help-gsl] gsl integration error

[Silke Diedenhofen]

Hello!

I have a problem running my C programm. I have included the gsl libraray to
calculate the Airy function, the derivative of the airy function and a qagi
integral. When I run the programm, the following error occurs:

gsl: ../integration/qags.c:543: ERROR: number of iterations was insufficient
Default GSL error handler invoked.

abnormal program termination

I have compiled the source code with DEV-cpp. This is a part of the source
code:

void Delta_Function(double lambda_m, double n_im, double E_field, double
*delta_n, double *deltaalpha){
     double q, m1, m2, m1_stern, m2_stern, h_quer, dAi_1, dAi_2, m0 ;
     double theta_1, theta_2, zeta_1, zeta_2, tmp_2_1, tmp_2_2, dA1_1,
dA1_2, tmp_3_1, tmp_3_2, alpha;
     m0 = 9.1095e-31;
     q = 1.6 * pow(10, -19);
     m1_stern = 0.058 * m0;
     m2_stern = 0.037 * m0;
     m1 = 0.084 * m0;
     m2 = 0.45 * m0;
     h_quer = H/(2*PI);
     theta_1 = pow((pow(q, 2) * pow(E_field,2))/(2*m1_stern*h_quer) , 1/3);
     theta_2 = pow((pow(q, 2) * pow(E_field,2))/(2*m2_stern*h_quer) , 1/3);
     zeta_1 = ((H*C/0.00000087)-(H*C/lambda_m))/(h_quer * theta_1);
     zeta_2 = ((H*C/0.00000087)-(H*C/lambda_m))/(h_quer * theta_2);
     tmp_2_1 = zeta_1 * pow(gsl_sf_airy_Ai(zeta_1, GSL_PREC_SINGLE),2);
     tmp_2_2 = zeta_2 * pow(gsl_sf_airy_Ai(zeta_2, GSL_PREC_SINGLE),2);
     dA1_1 = pow(gsl_sf_airy_Ai_deriv(zeta_1, GSL_PREC_SINGLE), 2);
     dA1_2 = pow(gsl_sf_airy_Ai_deriv(zeta_2, GSL_PREC_SINGLE), 2);
     tmp_3_1 = pow(m1_stern, 3/2) * (1+(m0/m1))* pow(zeta_1, 1/2) * (dA1_1
- tmp_2_1);
     tmp_3_1 = pow(m2_stern, 3/2) * (1+(m0/m2))* pow(zeta_2, 1/2) * (dA1_2
- tmp_2_2);     
     alpha = 1.3 * pow(10, 60) * (lambda_m/(2*C)) * (tmp_3_1 + tmp_3_2);
     *deltaalpha = alpha - (n_im * 4 * PI/(lambda_m /100));
     *delta_n = delta_n_calc(lambda_m, *deltaalpha);
     printf("%f\t %f\t", delta_n, deltaalpha);
}
     
double delta_n_calc(double lambda_m, double deltaalpha ) { 
  double delta_n;
  double f (double W, void * params) {
  double x = *(double *) params;
  double f = x / (pow(W, 2)-pow(H*C/lambda_m , 2));
  return f; }   
       double error, result;
       double x = deltaalpha;
       double h_quer = H/(2*PI); 
       gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
        gsl_function F;
        F.function = &f;
        F.params = &x;
  
 gsl_integration_qagi(&F, 0, 1e-7, 1000, w, &result, &error);
delta_n = 2*h_quer*C*result;  
return delta_n;  
}
 
Do you know what I have done wrong?

Many thanks.
Silke