voigt function algorithm for GSL attached: First contact

[michael crescimanno <mcrescim at cc dot ysu dot edu>]

Hi GSL friends,

   I have enjoyed using GSL and recently for a project needed a routine to
compute the voigt function and its first moment. I've checked the
attached code with IDL's voigt (which has only the function, no moments; IDL
agrees with the attached up to 7 digits over large range) and would like to
know if it might be worthwhile to submit it to GSL for inclusion in an
updated version since I believe the present package doesn't have a voigt
function. Code here is based on a relatively new 1994 released series by
Gubner, though it may have older antecedents in the BC (before computer)
literature.

   Let me know what I need to do to bring this up to the high level 
(and to put it in
the correct format; I'm just a physicist and not much of a C-programmer) 
I've found in
other GSL routines and I can resubmit. Peace,

michael



--
====================================================================
Michael Crescimanno                              mcrescim@cc.ysu.edu
Professor of Physics                             Tel: (330)-941-7109
Department of Physics and Astronomy
Youngstown State University                      Fax: (330)-941-3121
Youngstown, OH 44555-2001	     http://www.as.ysu.edu/~mcrescim
====================================================================

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
// Voigt.c : computes voigt function and first moment of the voigt function. 
//
//   voight is from convolution of gaussian and lorentzian
/*   This algorithm based on John A Gubner's 
    "A new series for approximating the Voigt Function" J. Phys. A: 
    Math. Gen. {\bf 27} L745-L749 (1994). Very useful and converges fast. 
*/ 
/* voight is defined as h(z,\alpha) = {{1}\over{\pi \sqrt{2\pi}} 
\int_-\infty^\infty dt e^{-(z-\alpha t)^2/2} /(1+t^2) 
and gubner has shown that 
h_{W,N} (z,\alpha) = \sqrt{2/pi} 1/W \sum_{-N}^{N} e^{-2(pi/W)^2 m^2} I_n(z,\lpha, W) 
in which I_n(z,\alpha, W) is defined by
I_n(z, \alpha, W) = Re [ (1-(-1)^n e^{-W(\alpha+iz)/2})/(\alpha + i(z-n 2\pi/W))]. 
*/ 
/* First Moment; let h(z,\alpha) \Delta_avg =  {{1}\over{\pi \sqrt{2\pi}} 
\int_-\infty^\infty dt t e^{-(z-\alpha t)^2/2} /(1+t^2) implies in terms of h defined above and normalized 

\Delta_avg = 1/\alpha [z + {{\partial_z h(z,\alpha)}\over{h}}]
and h and \Delta_avg is what this routine computes. 
*/ 

int

main (int argc, char *argv[])
{
  int N=1000, j; 
  double z, pi=4.*atan(1.0), w , sign, q, h, Deltaavg, alpha, converge2;
  gsl_complex A, B, C, D, X, Y, i, rho01,rho02, rho12, O1, O2, previous1, previous2, temp1, sum1, sum2;
  i=gsl_complex_rect(0.,1.);
  //inpute section; from command line
  z = -(double)(atoi(argv[1])/1000.);  // in the atomic context, this is v_0/sigma_MB, where v_0 is the laser detuning and sigma_MB is the maxwell-boltzmann width
  alpha = (double)(atoi(argv[2])/1000.); //this is the lorenztian half-width (gamma) divided by the sigma_MB. 
  sum1=gsl_complex_rect(0.,0.); //for h
  sum2=sum1; // for \Delta_avg
  sign = 1.;
  w = 150.;
  q = 2.*pi*pi/w/w;
  Y = gsl_complex_rect(alpha,z);
  C = gsl_complex_exp(gsl_complex_mul_real(Y, -w/2.));
  for (j=-N; j<N+1; j++)
    {
      sign = pow(-1.,j); 
      D = gsl_complex_rect(alpha, z-2.*pi*j/w);
      A = gsl_complex_add_real(gsl_complex_mul_real(C, -sign), 1.);
      // compute h
      A = gsl_complex_mul_real(gsl_complex_div(A, D), exp(-q*j*j)); 
      sum1 = gsl_complex_add(sum1, A);
      //compute \Delta_avg
      A = gsl_complex_add_real(gsl_complex_mul_real(D,w/2.),1.);
      B = gsl_complex_mul(gsl_complex_mul_real(C, sign), A);
      D = gsl_complex_div(gsl_complex_mul(D,D),i);
      B = gsl_complex_div(gsl_complex_add_real(B,-1.), D);
      B = gsl_complex_mul_real(B, exp(-q*j*j));
      sum2 = gsl_complex_add(sum2, B);
    }
  // POST PROCESSING AREA
  h = GSL_REAL(sum1);
  Deltaavg = (z+GSL_REAL(sum2)/h)/alpha; 
  Deltaavg = Deltaavg+z/alpha; //remove the overall detuning to find the excess from the average atom's point-of-view. 
  h = h*sqrt(2./pi)/w;
  h = sqrt(2.*pi)*h; // this line to make it agree with the convention used in IDL, for comparison only...remove? 
  converge2 = GSL_REAL(B)/GSL_REAL(sum2);
  printf("%3.5e   %3.5e    %3.5e   %3.5e   %3.5e\n",z, alpha,  h, Deltaavg, converge2);
  return 0;
}