program test_ze1
  use fvn_fnlib
  implicit none
  
  complex(8) :: z,resz
  complex(8), dimension(6) :: abramowitz_x,abramowitz_y
  integer(4) :: i
  
  abramowitz_x(1)=dcmplx(5.,5.)
  abramowitz_y(1)=dcmplx(0.906058d0,0.070209d0)
  abramowitz_x(2)=dcmplx(-5.,5.)
  abramowitz_y(2)=dcmplx(1.079407d0,0.143879d0)
  abramowitz_x(3)=dcmplx(11.,5.)
  abramowitz_y(3)=dcmplx(0.932672d0,0.026361d0)
  abramowitz_x(4)=dcmplx(-11.,5.)
  abramowitz_y(4)=dcmplx(1.084526d0,0.049336d0)
  abramowitz_x(5)=dcmplx(15.,18.)
  abramowitz_y(5)=dcmplx(0.972359d0,0.029448d0)
  abramowitz_x(6)=dcmplx(-15.,18.)
  abramowitz_y(6)=dcmplx(1.026183d0,0.036552d0)
  
  ! Value of z*exp(z)*e1(z) from Abramowitz & Stegun
  ! to cover the different cases of the algo
  !
  ! Case (5,5) modulus <10 with positive real part
  ! Case (-5,5) modulus <10 with negative real part
  ! Case (11,5) modulus between 10 and 20 with positive real part
  ! Case (-11,5) modulus between 10 and 20 with negative real part
  ! Case (15,18) modulus > 20 with positive real part
  ! Case (-15,18) modulus > 20 with negative real part
  
  do i=1,6
    z=abramowitz_x(i)
    resz=z*exp(z)*e1(z)
    write(*,*) "Tabulated : ",abramowitz_y(i)
    write(*,*) "Calculated :",resz
    write(*,*)
  end do
  end program