program test_ze1
  use fvn_fnlib
  implicit none

  complex(kind=dp_kind) :: z,resz
  complex(kind=dp_kind), dimension(6) :: abramowitz_x,abramowitz_y
  integer(kind=ip_kind) :: i

  abramowitz_x(1)=cmplx(5.d0,5.d0,dp_kind)
  abramowitz_y(1)=cmplx(0.906058d0,0.070209d0,dp_kind)
  abramowitz_x(2)=cmplx(-5.d0,5.d0,dp_kind)
  abramowitz_y(2)=cmplx(1.079407d0,0.143879d0,dp_kind)
  abramowitz_x(3)=cmplx(11.d0,5.d0,dp_kind)
  abramowitz_y(3)=cmplx(0.932672d0,0.026361d0,dp_kind)
  abramowitz_x(4)=cmplx(-11.d0,5.d0,dp_kind)
  abramowitz_y(4)=cmplx(1.084526d0,0.049336d0,dp_kind)
  abramowitz_x(5)=cmplx(15.d0,18.d0,dp_kind)
  abramowitz_y(5)=cmplx(0.972359d0,0.029448d0,dp_kind)
  abramowitz_x(6)=cmplx(-15.d0,18.d0,dp_kind)
  abramowitz_y(6)=cmplx(1.026183d0,0.036552d0,dp_kind)

  
  ! 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