module fvn
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! fvn : a f95 module replacement for some imsl routines
  ! it uses lapack for linear algebra
  ! it uses modified quadpack for integration
  !

  ! William Daniau 2007->today

  ! william.daniau@femto-st.fr
  !
  ! Routines naming scheme :
  !
  !           fvn_x_name
  !           where x can be  s : real 
  !                           d : real double precision
  !                           c : complex
  !                           z : double complex
  !
  !
  ! This piece of code is totally free! Do whatever you want with it. However

  ! if you find it usefull it would be kind to give credits ;-) 

  !

  ! svn version

  ! February 2008 : added fnlib to repository so there's no use to have
  !                 special functions and trigonometry here. Some functions
  !                 are then removed

  ! January 2008 : added quadratic interpolation, gamma/factorial function,
  !                 a function which return identity matrix,
  !                 evaluation of nterm chebyshev series

  ! September 2007 : added sparse system solving by interfacing umfpack

  ! June 2007 : added some complex trigonometric functions

  !
  ! TO DO LIST :
  ! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm 
  ! eigenvalues are given with no particular order.
  ! + Generic interface for fvn_x_name family  -> fvn_name
  ! + Make some parameters optional, status for example
  ! + use f95 kinds "double complex" -> complex(kind=8)
  ! + unify quadpack routines
  ! + ...
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  implicit none

  ! We define pi and i for the module
  real(kind=8),parameter :: fvn_pi = 3.141592653589793_8
  complex(kind=8),parameter :: fvn_i = (0._8,1._8)

  ! All quadpack routines are private to the module
  private ::  d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &
              dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, &
              dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, &
              dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt

  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Generic interface Definition
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  ! Identity Matrix
  interface fvn_ident
        module procedure fvn_s_ident,fvn_d_ident,fvn_c_ident,fvn_z_ident
  end interface fvn_ident
  
  ! Matrix inversion
  interface fvn_matinv
        module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv
  end interface fvn_matinv
  
  ! Determinant
  interface fvn_det
      module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det
  end interface fvn_det
  
  ! Condition
  interface fvn_matcon
      module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon
  end interface fvn_matcon
  
  ! Eigen
  interface fvn_matev
      module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev
  end interface fvn_matev
  
  ! Utility procedure find interval
  interface fvn_find_interval
      module procedure fvn_s_find_interval,fvn_d_find_interval
  end interface fvn_find_interval
  
  ! Quadratic 1D interpolation
  interface fvn_quad_interpol
      module procedure fvn_s_quad_interpol,fvn_d_quad_interpol
  end interface fvn_quad_interpol
  
  ! Quadratic 2D interpolation
  interface fvn_quad_2d_interpol
      module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol
  end interface fvn_quad_2d_interpol
  
  ! Quadratic 3D interpolation
  interface fvn_quad_3d_interpol
      module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol
  end interface fvn_quad_3d_interpol
  
  ! Akima interpolation
  interface fvn_akima
      module procedure fvn_s_akima,fvn_d_akima
  end interface fvn_akima
  
  ! Akima evaluation
  interface fvn_spline_eval
      module procedure fvn_s_spline_eval,fvn_d_spline_eval
  end interface fvn_spline_eval
  
  ! Least square polynomial
  interface fvn_lspoly
      module procedure fvn_s_lspoly,fvn_d_lspoly
  end interface fvn_lspoly
  
  ! Muller
  interface fvn_muller
      module procedure fvn_z_muller
  end interface fvn_muller
  
  ! Gauss legendre
  interface fvn_gauss_legendre
      module procedure fvn_d_gauss_legendre
  end interface fvn_gauss_legendre
  
  ! Simple Gauss Legendre integration
  interface fvn_gl_integ
      module procedure fvn_d_gl_integ
  end interface fvn_gl_integ
  
  ! Adaptative Gauss Kronrod integration f(x)
  interface fvn_integ_1_gk
      module procedure fvn_d_integ_1_gk
  end interface fvn_integ_1_gk
  
  ! Adaptative Gauss Kronrod integration f(x,y)
  interface fvn_integ_2_gk
      module procedure fvn_d_integ_2_gk
  end interface fvn_integ_2_gk
  
  ! Sparse solving
  interface fvn_sparse_solve
      module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve
  end interface fvn_sparse_solve

  contains 

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Identity Matrix
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function fvn_d_ident(n)
        implicit none
        integer(kind=4) :: n
        real(kind=8), dimension(n,n) :: fvn_d_ident
  
        real(kind=8),dimension(n*n) :: vect
        integer(kind=4) :: i
  
        vect=0._8
        vect(1:n*n:n+1) = 1._8
        fvn_d_ident=reshape(vect, shape = (/ n,n /))
  end function
  
  function fvn_s_ident(n)
        implicit none
        integer(kind=4) :: n
        real(kind=4), dimension(n,n) :: fvn_s_ident
  
        real(kind=4),dimension(n*n) :: vect
        integer(kind=4) :: i
  
       vect=0._4
        vect(1:n*n:n+1) = 1._4
        fvn_s_ident=reshape(vect, shape = (/ n,n /))
  end function
  
  function fvn_c_ident(n)
        implicit none
        integer(kind=4) :: n
        complex(kind=4), dimension(n,n) :: fvn_c_ident
  
        complex(kind=4),dimension(n*n) :: vect
        integer(kind=4) :: i
  
        vect=(0._4,0._4)
        vect(1:n*n:n+1) = (1._4,0._4)
        fvn_c_ident=reshape(vect, shape = (/ n,n /))
  end function
  
  function fvn_z_ident(n)
        implicit none
        integer(kind=4) :: n
        complex(kind=8), dimension(n,n) :: fvn_z_ident
  
        complex(kind=8),dimension(n*n) :: vect
        integer(kind=4) :: i
  
        vect=(0._8,0._8)
        vect(1:n*n:n+1) = (1._8,0._8)
        fvn_z_ident=reshape(vect, shape = (/ n,n /))
  end function

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Matrix inversion subroutines
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine fvn_s_matinv(d,a,inva,status)
      !
      ! Matrix inversion of a real matrix using BLAS and LAPACK
      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !

      implicit none

      integer, intent(in) :: d
      real, intent(in) :: a(d,d)
      real, intent(out) :: inva(d,d)

      integer, intent(out),optional :: status

  
      integer, allocatable :: ipiv(:)
      real, allocatable :: work(:)
      real twork(1)
      integer :: info
      integer :: lwork

      if (present(status)) status=1

  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call scopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      ! LU factorization using LAPACK
      call sgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call sgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call sgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  subroutine fvn_d_matinv(d,a,inva,status)
      !
      ! Matrix inversion of a double precision matrix using BLAS and LAPACK
      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !

      implicit none

      integer, intent(in), optional :: d

      double precision, intent(in) :: a(d,d)
      double precision, intent(out) :: inva(d,d)

      integer, intent(out),optional :: status

  
      integer, allocatable :: ipiv(:)
      double precision, allocatable :: work(:)
      double precision :: twork(1)
      integer :: info
      integer :: lwork

      if (present(status)) status=1

  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call dcopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      ! LU factorization using LAPACK
      call dgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call dgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call dgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  subroutine fvn_c_matinv(d,a,inva,status)
      !
      ! Matrix inversion of a complex matrix using BLAS and LAPACK
      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !

      implicit none

      integer, intent(in) :: d
      complex, intent(in) :: a(d,d)
      complex, intent(out) :: inva(d,d)

      integer, intent(out),optional :: status

  
      integer, allocatable :: ipiv(:)
      complex, allocatable :: work(:)
      complex :: twork(1)
      integer :: info
      integer :: lwork

       if (present(status)) status=1

  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call ccopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      
      ! LU factorization using LAPACK
      call cgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call cgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call cgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  subroutine fvn_z_matinv(d,a,inva,status)
      !
      ! Matrix inversion of a double complex matrix using BLAS and LAPACK
      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !

      implicit none

      integer, intent(in) :: d
      double complex, intent(in) :: a(d,d)
      double complex, intent(out) :: inva(d,d)

      integer, intent(out),optional :: status

  
      integer, allocatable :: ipiv(:)
      double complex, allocatable :: work(:)
      double complex :: twork(1)
      integer :: info
      integer :: lwork

       if (present(status)) status=1

  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call zcopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      
      ! LU factorization using LAPACK
      call zgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call zgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call zgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then

          if (present(status)) status=0

      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Determinants
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function fvn_s_det(d,a,status)
      !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! status (out) : =0 if LU factorization failed
      !

      implicit none

      integer, intent(in) :: d
      real, intent(in) :: a(d,d)

      integer, intent(out), optional :: status

      real :: fvn_s_det
      
      real, allocatable :: wc_a(:,:)
      integer, allocatable :: ipiv(:)
      integer :: info,i
      

       if (present(status)) status=1

      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call sgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then

           if (present(status)) status=0

          fvn_s_det=0.e0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_s_det=1.e0
      do i=1,d
          if (ipiv(i)==i) then
              fvn_s_det=fvn_s_det*wc_a(i,i)
          else
              fvn_s_det=-fvn_s_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  function fvn_d_det(d,a,status)
      !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! status (out) : =0 if LU factorization failed
      !

      implicit none

      integer, intent(in) :: d
      double precision, intent(in) :: a(d,d)

      integer, intent(out), optional :: status

      double precision :: fvn_d_det
      
      double precision, allocatable :: wc_a(:,:)
      integer, allocatable :: ipiv(:)
      integer :: info,i
      

       if (present(status)) status=1

      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call dgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then

           if (present(status)) status=0

          fvn_d_det=0.d0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_d_det=1.d0
      do i=1,d
          if (ipiv(i)==i) then
              fvn_d_det=fvn_d_det*wc_a(i,i)
          else
              fvn_d_det=-fvn_d_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  function fvn_c_det(d,a,status)    !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! status (out) : =0 if LU factorization failed
      !

      implicit none

      integer, intent(in) :: d
      complex, intent(in) :: a(d,d)

      integer, intent(out), optional :: status

      complex :: fvn_c_det
      
      complex, allocatable :: wc_a(:,:)
      integer, allocatable :: ipiv(:)
      integer :: info,i
      

       if (present(status)) status=1

      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call cgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then

           if (present(status)) status=0

          fvn_c_det=(0.e0,0.e0)
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_c_det=(1.e0,0.e0)
      do i=1,d
          if (ipiv(i)==i) then
              fvn_c_det=fvn_c_det*wc_a(i,i)
          else
              fvn_c_det=-fvn_c_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  function fvn_z_det(d,a,status)
      !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! det (out) : determinant
      ! status (out) : =0 if LU factorization failed
      !

      implicit none

      integer, intent(in) :: d
      double complex, intent(in) :: a(d,d)

      integer, intent(out), optional :: status

      double complex :: fvn_z_det
      
      double complex, allocatable :: wc_a(:,:)
      integer, allocatable :: ipiv(:)
      integer :: info,i
      

       if (present(status)) status=1

      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call zgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then

           if (present(status)) status=0

          fvn_z_det=(0.d0,0.d0)
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_z_det=(1.d0,0.d0)
      do i=1,d
          if (ipiv(i)==i) then
              fvn_z_det=fvn_z_det*wc_a(i,i)
          else
              fvn_z_det=-fvn_z_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Condition test
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! 1-norm 
  ! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
  ! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
  !
  subroutine fvn_s_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !

      implicit none

      integer, intent(in) :: d
      real, intent(in) :: a(d,d)
      real, intent(out) :: rcond

      integer, intent(out), optional :: status

      
      real, allocatable :: work(:)
      integer, allocatable :: iwork(:)
      real :: anorm
      real, allocatable :: wc_a(:,:) ! working copy of a
      integer :: info
      integer, allocatable :: ipiv(:)
      
      real, external :: slange
      
      

       if (present(status)) status=1

      
      anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call scopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call sgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then

           if (present(status)) status=0

          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      allocate(work(4*d))
      allocate(iwork(d))
      call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
      if (info /= 0) then

           if (present(status)) status=0

      end if
      deallocate(iwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  
  end subroutine
  
  subroutine fvn_d_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !

      implicit none

      integer, intent(in) :: d
      double precision, intent(in) :: a(d,d)
      double precision, intent(out) :: rcond

      integer, intent(out), optional :: status

      
      double precision, allocatable :: work(:)
      integer, allocatable :: iwork(:)
      double precision :: anorm
      double precision, allocatable :: wc_a(:,:) ! working copy of a
      integer :: info
      integer, allocatable :: ipiv(:)
      
      double precision, external :: dlange
      
      

       if (present(status)) status=1

      
      anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call dcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call dgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then

           if (present(status)) status=0

          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
  
      allocate(work(4*d))
      allocate(iwork(d))
      call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
      if (info /= 0) then

           if (present(status)) status=0

      end if
      deallocate(iwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  
  end subroutine
  
  subroutine fvn_c_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !

      implicit none

      integer, intent(in) :: d
      complex, intent(in) :: a(d,d)
      real, intent(out) :: rcond

      integer, intent(out), optional :: status

      
      real, allocatable :: rwork(:)
      complex, allocatable :: work(:)
      integer, allocatable :: iwork(:)
      real :: anorm
      complex, allocatable :: wc_a(:,:) ! working copy of a
      integer :: info
      integer, allocatable :: ipiv(:)
      
      real, external :: clange
      
      

       if (present(status)) status=1

      
      anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call ccopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call cgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then

           if (present(status)) status=0

          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      allocate(work(2*d))
      allocate(rwork(2*d))
      call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
      if (info /= 0) then

           if (present(status)) status=0

      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  end subroutine
  
  subroutine fvn_z_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !

      implicit none

      integer, intent(in) :: d
      double complex, intent(in) :: a(d,d)
      double precision, intent(out) :: rcond

      integer, intent(out), optional :: status

      
      double complex, allocatable :: work(:)
      double precision, allocatable :: rwork(:)
      double precision :: anorm
      double complex, allocatable :: wc_a(:,:) ! working copy of a
      integer :: info
      integer, allocatable :: ipiv(:)
      
      double precision, external :: zlange
      
      

       if (present(status)) status=1

      
      anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call zcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call zgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then

           if (present(status)) status=0

          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
  
      allocate(work(2*d))
      allocate(rwork(2*d))
      call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
      if (info /= 0) then

           if (present(status)) status=0

      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  end subroutine
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Valeurs propres/ Vecteurs propre
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  subroutine fvn_s_matev(d,a,evala,eveca,status)
      ! 
      ! integer d (in) : matrice rank
      ! real a(d,d) (in) : The Matrix
      ! complex evala(d) (out) : eigenvalues
      ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer (out) : status =0 if something went wrong
      !
      ! interfacing Lapack routine SGEEV

      implicit none

      integer, intent(in) :: d
      real, intent(in) :: a(d,d)
      complex, intent(out) :: evala(d)
      complex, intent(out) :: eveca(d,d)

      integer, intent(out), optional :: status

      
      real, allocatable :: wc_a(:,:)  ! a working copy of a
      integer :: info
      integer :: lwork
      real, allocatable :: wr(:),wi(:)
      real :: vl      ! unused but necessary for the call
      real, allocatable :: vr(:,:)
      real, allocatable :: work(:)
      real :: twork(1)
      integer i
      integer j
      

       if (present(status)) status=1

      ! making a working copy of a
      allocate(wc_a(d,d))
      !call scopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(wr(d))
      allocate(wi(d))
      allocate(vr(d,d))
      ! query optimal work size
      call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
  
      if (info /= 0) then

           if (present(status)) status=0

          deallocate(work)
          deallocate(vr)
          deallocate(wi)
          deallocate(wr)
          deallocate(wc_a)
          return
      end if
  
      ! now fill in the results
      i=1
      do while(i<=d)
          evala(i)=cmplx(wr(i),wi(i))
          if (wi(i) == 0.) then ! eigenvalue is real
              eveca(:,i)=cmplx(vr(:,i),0.)
          else ! eigenvalue is complex
              evala(i+1)=cmplx(wr(i+1),wi(i+1))
              eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
              eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
              i=i+1
          end if
          i=i+1
      enddo
      deallocate(work)
      deallocate(vr)
      deallocate(wi)
      deallocate(wr)
      deallocate(wc_a)
  
  end subroutine
  
  subroutine fvn_d_matev(d,a,evala,eveca,status)
      ! 
      ! integer d (in) : matrice rank
      ! double precision a(d,d) (in) : The Matrix
      ! double complex evala(d) (out) : eigenvalues
      ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer (out) : status =0 if something went wrong
      !
      ! interfacing Lapack routine DGEEV

      implicit none

      integer, intent(in) :: d
      double precision, intent(in) :: a(d,d)
      double complex, intent(out) :: evala(d)
      double complex, intent(out) :: eveca(d,d)

      integer, intent(out), optional :: status

      
      double precision, allocatable :: wc_a(:,:)  ! a working copy of a
      integer :: info
      integer :: lwork
      double precision, allocatable :: wr(:),wi(:)
      double precision :: vl      ! unused but necessary for the call
      double precision, allocatable :: vr(:,:)
      double precision, allocatable :: work(:)
      double precision :: twork(1)
      integer i
      integer j
      

      if (present(status)) status=1

      ! making a working copy of a
      allocate(wc_a(d,d))
      !call dcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(wr(d))
      allocate(wi(d))
      allocate(vr(d,d))
      ! query optimal work size
      call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
  
      if (info /= 0) then

           if (present(status)) status=0

          deallocate(work)
          deallocate(vr)
          deallocate(wi)
          deallocate(wr)
          deallocate(wc_a)
          return
      end if
  
      ! now fill in the results
      i=1
      do while(i<=d)
          evala(i)=dcmplx(wr(i),wi(i))
          if (wi(i) == 0.) then ! eigenvalue is real
              eveca(:,i)=dcmplx(vr(:,i),0.)
          else ! eigenvalue is complex
              evala(i+1)=dcmplx(wr(i+1),wi(i+1))
              eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
              eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
              i=i+1
          end if
          i=i+1
      enddo
  
      deallocate(work)
      deallocate(vr)
      deallocate(wi)
      deallocate(wr)
      deallocate(wc_a)
  
  end subroutine
  
  subroutine fvn_c_matev(d,a,evala,eveca,status)
      ! 
      ! integer d (in) : matrice rank
      ! complex a(d,d) (in) : The Matrix
      ! complex evala(d) (out) : eigenvalues
      ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer (out) : status =0 if something went wrong
      !
      ! interfacing Lapack routine CGEEV

      implicit none

      integer, intent(in) :: d
      complex, intent(in) :: a(d,d)
      complex, intent(out) :: evala(d)
      complex, intent(out) :: eveca(d,d)

      integer, intent(out), optional :: status

  
      complex, allocatable :: wc_a(:,:) ! a working copy of a
      integer :: info
      integer :: lwork
      complex, allocatable :: work(:)
      complex :: twork(1)
      real, allocatable :: rwork(:)
      complex :: vl   ! unused but necessary for the call

       if (present(status)) status=1

  
      ! making a working copy of a
      allocate(wc_a(d,d))
      !call ccopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      
      ! query optimal work size
      call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      allocate(rwork(2*d))
      call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
      
      if (info /= 0) then

           if (present(status)) status=0

      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(wc_a)
  
  end subroutine
  
  subroutine fvn_z_matev(d,a,evala,eveca,status)
      ! 
      ! integer d (in) : matrice rank
      ! double complex a(d,d) (in) : The Matrix
      ! double complex evala(d) (out) : eigenvalues
      ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer (out) : status =0 if something went wrong
      !
      ! interfacing Lapack routine ZGEEV

      implicit none

      integer, intent(in) :: d
      double complex, intent(in) :: a(d,d)
      double complex, intent(out) :: evala(d)
      double complex, intent(out) :: eveca(d,d)

      integer, intent(out), optional :: status

  
      double complex, allocatable :: wc_a(:,:) ! a working copy of a
      integer :: info
      integer :: lwork
      double complex, allocatable :: work(:)
      double complex :: twork(1)
      double precision, allocatable :: rwork(:)
      double complex :: vl   ! unused but necessary for the call
      

       if (present(status)) status=1

  
      ! making a working copy of a
      allocate(wc_a(d,d))
      !call zcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      ! query optimal work size
      call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      allocate(rwork(2*d))
      call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
      
      if (info /= 0) then

           if (present(status)) status=0

      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(wc_a)
  
  end subroutine

  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Quadratic interpolation of tabulated function of 1,2 or 3 variables
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  subroutine fvn_s_find_interval(x,i,xdata,n)
        implicit none
        ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
        ! xdata(n) must contains a set of increasingly ordered values 
        ! if x < xdata(1) i=0 is returned
        ! if x > xdata(n) i=n is returned
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        ! a simple dichotomy method is used
  
        real(kind=4), intent(in) :: x
        real(kind=4), intent(in), dimension(n) :: xdata
        integer(kind=4), intent(in) :: n
        integer(kind=4), intent(out) :: i
  
        integer(kind=4) :: imin,imax,imoyen
  
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        if (x == xdata(n)) then
              i=n-1
              return
        end if
  
        ! if x < xdata(1) i=0 is returned
        if (x < xdata(1)) then
              i=0
              return
        end if
  
        ! if x > xdata(n) i=n is returned
        if (x > xdata(n)) then
              i=n
              return
        end if
  
        ! here xdata(1) <= x <= xdata(n)
        imin=0
        imax=n+1
  
        do while((imax-imin) > 1)
              imoyen=(imax+imin)/2
              if (x >= xdata(imoyen)) then
                    imin=imoyen
              else
                    imax=imoyen
              end if
        end do
  
        i=imin
  
  end subroutine
  
  
  subroutine fvn_d_find_interval(x,i,xdata,n)
        implicit none
        ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
        ! xdata(n) must contains a set of increasingly ordered values 
        ! if x < xdata(1) i=0 is returned
        ! if x > xdata(n) i=n is returned
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        ! a simple dichotomy method is used
  
        real(kind=8), intent(in) :: x
        real(kind=8), intent(in), dimension(n) :: xdata
        integer(kind=4), intent(in) :: n
        integer(kind=4), intent(out) :: i
  
        integer(kind=4) :: imin,imax,imoyen
  
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        if (x == xdata(n)) then
              i=n-1
              return
        end if
  
        ! if x < xdata(1) i=0 is returned
        if (x < xdata(1)) then
              i=0
              return
        end if
  
        ! if x > xdata(n) i=n is returned
        if (x > xdata(n)) then
              i=n
              return
        end if
  
        ! here xdata(1) <= x <= xdata(n)
        imin=0
        imax=n+1
  
        do while((imax-imin) > 1)
              imoyen=(imax+imin)/2
              if (x >= xdata(imoyen)) then
                    imin=imoyen
              else
                    imax=imoyen
              end if
        end do
  
        i=imin
  
  end subroutine
  
  
  function fvn_s_quad_interpol(x,n,xdata,ydata)
        implicit none
        ! This function evaluate the value of a function defined by a set of points
        ! and values, using a quadratic interpolation
        ! xdata must be increasingly ordered
        ! x must be within xdata(1) and xdata(n) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: n
        real(kind=4), intent(in), dimension(n) :: xdata,ydata
        real(kind=4), intent(in) :: x
        real(kind=4) ::  fvn_s_quad_interpol
  
        integer(kind=4) :: iinf,base,i,j
        real(kind=4) :: p
  
        call fvn_s_find_interval(x,iinf,xdata,n)
  
        ! Settings for extrapolation
        if (iinf==0) then
              ! TODO -> Lower bound extrapolation warning
              iinf=1
        end if
  
        if (iinf==n) then
              ! TODO -> Higher bound extrapolation warning
              iinf=n-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (iinf==1) then
              base=0
        else
              base=iinf-2
        end if
  
        ! The three points we will use are :
        ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
  
        ! Straight forward Lagrange polynomial
        fvn_s_quad_interpol=0.
        do i=1,3
              !  polynome i
              p=ydata(base+i)
              do j=1,3
                    if (j /= i) then
                          p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
                    end if
              end do
              fvn_s_quad_interpol=fvn_s_quad_interpol+p
        end do
  
  end function
  
  
  function fvn_d_quad_interpol(x,n,xdata,ydata)
        implicit none
        ! This function evaluate the value of a function defined by a set of points
        ! and values, using a quadratic interpolation
        ! xdata must be increasingly ordered
        ! x must be within xdata(1) and xdata(n) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: n
        real(kind=8), intent(in), dimension(n) :: xdata,ydata
        real(kind=8), intent(in) :: x
        real(kind=8) ::  fvn_d_quad_interpol
  
        integer(kind=4) :: iinf,base,i,j
        real(kind=8) :: p
  
        call fvn_d_find_interval(x,iinf,xdata,n)
  
        ! Settings for extrapolation
        if (iinf==0) then
              ! TODO -> Lower bound extrapolation warning
              iinf=1
        end if
  
        if (iinf==n) then
              ! TODO Higher bound extrapolation warning
              iinf=n-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (iinf==1) then
              base=0
        else
              base=iinf-2
        end if
  
        ! The three points we will use are :
        ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
  
        ! Straight forward Lagrange polynomial
        fvn_d_quad_interpol=0.
        do i=1,3
              !  polynome i
              p=ydata(base+i)
              do j=1,3
                    if (j /= i) then
                          p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
                    end if
              end do
              fvn_d_quad_interpol=fvn_d_quad_interpol+p
        end do
  
  end function
  
  
  function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
        implicit none
        ! This function evaluate the value of a two variable function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata and ydata must be increasingly ordered
        ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
        ! y within ydata(1) and ydata(ny) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny
        real(kind=4), intent(in) :: x,y
        real(kind=4), intent(in), dimension(nx) :: xdata
        real(kind=4), intent(in), dimension(ny) :: ydata
        real(kind=4), intent(in), dimension(nx,ny) :: zdata
        real(kind=4) :: fvn_s_quad_2d_interpol
  
        integer(kind=4) :: ixinf,iyinf,basex,basey,i
        real(kind=4),dimension(3) :: ztmp
        !real(kind=4), external :: fvn_s_quad_interpol
  
        call fvn_s_find_interval(x,ixinf,xdata,nx)
        call fvn_s_find_interval(y,iyinf,ydata,ny)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
        ! stored in ztmp(1:3)
        do i=1,3
              ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i))
        end do
  
        ! Then we make an interpolation for y using previous interpolations
        fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
  end function
  
  
  function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
        implicit none
        ! This function evaluate the value of a two variable function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata and ydata must be increasingly ordered
        ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
        ! y within ydata(1) and ydata(ny) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny
        real(kind=8), intent(in) :: x,y
        real(kind=8), intent(in), dimension(nx) :: xdata
        real(kind=8), intent(in), dimension(ny) :: ydata
        real(kind=8), intent(in), dimension(nx,ny) :: zdata
        real(kind=8) :: fvn_d_quad_2d_interpol
  
        integer(kind=4) :: ixinf,iyinf,basex,basey,i
        real(kind=8),dimension(3) :: ztmp
        !real(kind=8), external :: fvn_d_quad_interpol
  
        call fvn_d_find_interval(x,ixinf,xdata,nx)
        call fvn_d_find_interval(y,iyinf,ydata,ny)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
        ! stored in ztmp(1:3)
        do i=1,3
              ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i))
        end do
  
        ! Then we make an interpolation for y using previous interpolations
        fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
  end function
  
  
  function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
        implicit none
        ! This function evaluate the value of a 3 variables function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata, ydata and zdata must be increasingly ordered
        ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 
        ! perform an interpolation, otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny,nz
        real(kind=4), intent(in) :: x,y,z
        real(kind=4), intent(in), dimension(nx) :: xdata
        real(kind=4), intent(in), dimension(ny) :: ydata
        real(kind=4), intent(in), dimension(nz) :: zdata
        real(kind=4), intent(in), dimension(nx,ny,nz) :: tdata
        real(kind=4) :: fvn_s_quad_3d_interpol
  
        integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
        !real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol
        real(kind=4),dimension(3,3) :: ttmp
  
        call fvn_s_find_interval(x,ixinf,xdata,nx)
        call fvn_s_find_interval(y,iyinf,ydata,ny)
        call fvn_s_find_interval(z,izinf,zdata,nz)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        if (izinf==0) then
              ! TODO -> Lower z bound extrapolation warning
              izinf=1
        end if
  
        if (izinf==nz) then
              ! TODO -> Higher z bound extrapolation warning
              izinf=nz-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        if (izinf==1) then
              basez=0
        else
              basez=izinf-2
        end if
  
        ! We first make 9 one dimensional interpolation on variable x.
        ! results are stored in ttmp
        do i=1,3
              do j=1,3
                    ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
              end do
        end do
  
        ! We then make a 2 dimensionnal interpolation on variables y and z
        fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &
              3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
  end function
  
  
  function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
        implicit none
        ! This function evaluate the value of a 3 variables function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata, ydata and zdata must be increasingly ordered
        ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 
        ! perform an interpolation, otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny,nz
        real(kind=8), intent(in) :: x,y,z
        real(kind=8), intent(in), dimension(nx) :: xdata
        real(kind=8), intent(in), dimension(ny) :: ydata
        real(kind=8), intent(in), dimension(nz) :: zdata
        real(kind=8), intent(in), dimension(nx,ny,nz) :: tdata
        real(kind=8) :: fvn_d_quad_3d_interpol
  
        integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
        !real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol
        real(kind=8),dimension(3,3) :: ttmp
  
        call fvn_d_find_interval(x,ixinf,xdata,nx)
        call fvn_d_find_interval(y,iyinf,ydata,ny)
        call fvn_d_find_interval(z,izinf,zdata,nz)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        if (izinf==0) then
              ! TODO -> Lower z bound extrapolation warning
              izinf=1
        end if
  
        if (izinf==nz) then
              ! TODO -> Higher z bound extrapolation warning
              izinf=nz-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        if (izinf==1) then
              basez=0
        else
              basez=izinf-2
        end if
  
        ! We first make 9 one dimensional interpolation on variable x.
        ! results are stored in ttmp
        do i=1,3
              do j=1,3
                    ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
              end do
        end do
  
        ! We then make a 2 dimensionnal interpolation on variables y and z
        fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &
              3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
  end function

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Akima spline interpolation and spline evaluation
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  ! Single precision
  subroutine fvn_s_akima(n,x,y,br,co)
      implicit none
      integer, intent(in)  :: n
      real, intent(in) :: x(n)
      real, intent(in) :: y(n)
      real, intent(out) :: br(n)
      real, intent(out) :: co(4,n)
      
      real, allocatable :: var(:),z(:)
      real :: wi_1,wi
      integer :: i
      real :: dx,a,b
  
      ! br is just a copy of x
      br(:)=x(:)
      
      allocate(var(n))
      allocate(z(n))
      ! evaluate the variations
      do i=1, n-1
          var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
      end do
      var(n+2)=2.e0*var(n+1)-var(n)
      var(n+3)=2.e0*var(n+2)-var(n+1)
      var(2)=2.e0*var(3)-var(4)
      var(1)=2.e0*var(2)-var(3)
    
      do i = 1, n
      wi_1=abs(var(i+3)-var(i+2))
      wi=abs(var(i+1)-var(i))
      if ((wi_1+wi).eq.0.e0) then
          z(i)=(var(i+2)+var(i+1))/2.e0
      else
          z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
      end if
      end do
      
      do i=1, n-1
          dx=x(i+1)-x(i)
          a=(z(i+1)-z(i))*dx                      ! coeff intermediaires pour calcul wd
          b=y(i+1)-y(i)-z(i)*dx                   ! coeff intermediaires pour calcul wd
          co(1,i)=y(i)
          co(2,i)=z(i)
          !co(3,i)=-(a-3.*b)/dx**2                ! méthode wd
          !co(4,i)=(a-2.*b)/dx**3                 ! méthode wd
          co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx   ! méthode JP Moreau
          co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2  !
          ! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6
          ! etrangement la fonction csval corrige et donne la bonne valeur ...
      end do
      co(1,n)=y(n)
      co(2,n)=z(n)
      co(3,n)=0.e0
      co(4,n)=0.e0
  
      deallocate(z)
      deallocate(var)
  
  end subroutine
  
  ! Double precision
  subroutine fvn_d_akima(n,x,y,br,co)
  
      implicit none
      integer, intent(in)  :: n
      double precision, intent(in) :: x(n)
      double precision, intent(in) :: y(n)
      double precision, intent(out) :: br(n)
      double precision, intent(out) :: co(4,n)
      
      double precision, allocatable :: var(:),z(:)
      double precision :: wi_1,wi
      integer :: i
      double precision :: dx,a,b
      
      ! br is just a copy of x
      br(:)=x(:)
  
      allocate(var(n))
      allocate(z(n))
      ! evaluate the variations
      do i=1, n-1
          var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
      end do
      var(n+2)=2.d0*var(n+1)-var(n)
      var(n+3)=2.d0*var(n+2)-var(n+1)
      var(2)=2.d0*var(3)-var(4)
      var(1)=2.d0*var(2)-var(3)
    
      do i = 1, n
      wi_1=dabs(var(i+3)-var(i+2))
      wi=dabs(var(i+1)-var(i))
      if ((wi_1+wi).eq.0.d0) then
          z(i)=(var(i+2)+var(i+1))/2.d0
      else
          z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
      end if
      end do
      
      do i=1, n-1
          dx=x(i+1)-x(i)
          a=(z(i+1)-z(i))*dx                      ! coeff intermediaires pour calcul wd
          b=y(i+1)-y(i)-z(i)*dx                   ! coeff intermediaires pour calcul wd
          co(1,i)=y(i)
          co(2,i)=z(i)
          !co(3,i)=-(a-3.*b)/dx**2                ! méthode wd
          !co(4,i)=(a-2.*b)/dx**3                 ! méthode wd
          co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx   ! méthode JP Moreau
          co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2  !
          ! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6
          ! etrangement la fonction csval corrige et donne la bonne valeur ...
      end do
      co(1,n)=y(n)
      co(2,n)=z(n)
      co(3,n)=0.d0
      co(4,n)=0.d0
  
      deallocate(z)
      deallocate(var)
  
  end subroutine
  
  !
  ! Single precision spline evaluation
  !
  function fvn_s_spline_eval(x,n,br,co)
      implicit none
      real, intent(in) :: x           ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
      integer, intent(in) :: n        ! number of intervals
      real, intent(in) :: br(n+1)     ! breakpoints
      real, intent(in) :: co(4,n+1)   ! spline coeeficients
      real :: fvn_s_spline_eval
      
      integer :: i
      real :: dx
      
      if (x<=br(1)) then
          i=1
      else if (x>=br(n+1)) then
          i=n
      else
      i=1
      do while(x>=br(i))
          i=i+1
      end do
      i=i-1
      end if
      dx=x-br(i)
      fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
  
  end function
  
  ! Double precision spline evaluation
  function fvn_d_spline_eval(x,n,br,co)
      implicit none
      double precision, intent(in) :: x           ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
      integer, intent(in) :: n        ! number of intervals
      double precision, intent(in) :: br(n+1)     ! breakpoints
      double precision, intent(in) :: co(4,n+1)   ! spline coeeficients
      double precision :: fvn_d_spline_eval
      
      integer :: i
      double precision :: dx
      
      
      if (x<=br(1)) then
          i=1
      else if (x>=br(n+1)) then
          i=n
      else
      i=1
      do while(x>=br(i))
          i=i+1
      end do
      i=i-1
      end if
      
      dx=x-br(i)
      fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
  
  end function

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! 
  ! Least square problem
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  !

  subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
  !
  !   Least square polynomial fitting
  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none
  
  integer, intent(in) :: np,deg
  real(kind=8), intent(in), dimension(np) :: x,y
  real(kind=8), intent(out), dimension(deg+1) :: coeff

  integer, intent(out), optional :: status

  real(kind=8), allocatable, dimension(:,:) :: mat,bmat
  real(kind=8),dimension(:),allocatable :: work
  real(kind=8),dimension(1) :: twork
  integer :: lwork,info
  
  integer :: i,j

   if (present(status)) status=1

  allocate(mat(np,deg+1),bmat(np,1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
  lwork=twork(1)
  allocate(work(int(lwork)))
  ! real call
  call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
  
  if (info /= 0) then

       if (present(status)) status=0

  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat)
  end subroutine
  
  subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)
  !
  !   Least square polynomial fitting
  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none
  
  integer, intent(in) :: np,deg
  real(kind=4), intent(in), dimension(np) :: x,y
  real(kind=4), intent(out), dimension(deg+1) :: coeff

  integer, intent(out), optional :: status

  
  real(kind=4), allocatable, dimension(:,:) :: mat,bmat
  real(kind=4),dimension(:),allocatable :: work
  real(kind=4),dimension(1) :: twork
  integer :: lwork,info
  
  integer :: i,j

   if (present(status)) status=1

  allocate(mat(np,deg+1),bmat(np,1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
  lwork=twork(1)
  allocate(work(int(lwork)))
  ! real call
  call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
  
  if (info /= 0) then

       if (present(status)) status=0

  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat)
  end subroutine
  
  
  
  
  
  
  
  
  subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status)
  !

  !   Least square polynomial fitting using singular value decomposition

  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none
  
  integer, intent(in) :: np,deg
  real(kind=8), intent(in), dimension(np) :: x,y
  real(kind=8), intent(out), dimension(deg+1) :: coeff

  integer, intent(out), optional :: status

  real(kind=8), allocatable, dimension(:,:) :: mat,bmat
  real(kind=8),dimension(:),allocatable :: work,singval
  real(kind=8),dimension(1) :: twork
  integer :: lwork,info,rank
  
  integer :: i,j

   if (present(status)) status=1

  allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
  lwork=twork(1)
  allocate(work(int(lwork)))
  ! real call
  call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
  
  if (info /= 0) then

       if (present(status)) status=0

  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat,singval)
  end subroutine

  subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status)

  !

  !   Least square polynomial fitting using singular value decomposition

  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none
  
  integer, intent(in) :: np,deg
  real(kind=4), intent(in), dimension(np) :: x,y
  real(kind=4), intent(out), dimension(deg+1) :: coeff

  integer, intent(out), optional :: status

  
  real(kind=4), allocatable, dimension(:,:) :: mat,bmat
  real(kind=4),dimension(:),allocatable :: work,singval
  real(kind=4),dimension(1) :: twork
  integer :: lwork,info,rank
  
  integer :: i,j

   if (present(status)) status=1

  allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
  lwork=twork(1)
  allocate(work(int(lwork)))
  ! real call
  call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
  
  if (info /= 0) then

       if (present(status)) status=0

  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat,singval)
  end subroutine

  !
  ! Muller
  !
  !
  !
  ! William Daniau 2007
  !
  ! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
  ! http://plato.asu.edu/ftp/other_software/muller.f
  !
  ! it can be used as a replacement for imsl routine dzanly with minor changes
  !
  !-----------------------------------------------------------------------
  !
  !   purpose             - zeros of an analytic complex function
  !                           using the muller method with deflation
  !
  !   usage               - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
  !                           infer,ier)
  !
  !   arguments    f      - a complex function subprogram, f(z), written
  !                           by the user specifying the equation whose
  !                           roots are to be found.  f must appear in
  !                           an external statement in the calling pro-
  !                           gram.
  !                eps    - 1st stopping criterion.  let fp(z)=f(z)/p
  !                           where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
  !                           and z(1),...,z(k-1) are previously found
  !                           roots.  if ((cdabs(f(z)).le.eps) .and.
  !                           (cdabs(fp(z)).le.eps)), then z is accepted
  !                           as a root. (input)
  !                eps1   - 2nd stopping criterion.  a root is accepted
  !                           if two successive approximations to a given
  !                           root agree within eps1. (input)
  !                             note. if either or both of the stopping
  !                             criteria are fulfilled, the root is
  !                             accepted.
  !                kn     - the number of known roots which must be stored
  !                           in x(1),...,x(kn), prior to entry to muller
  !                nguess - the number of initial guesses provided. these
  !                           guesses must be stored in x(kn+1),...,
  !                           x(kn+nguess).  nguess must be set equal
  !                           to zero if no guesses are provided. (input)
  !                n      - the number of new roots to be found by
  !                           muller (input)
  !                x      - a complex vector of length kn+n.  x(1),...,
  !                           x(kn) on input must contain any known
  !                           roots.  x(kn+1),..., x(kn+n) on input may,
  !                           on user option, contain initial guesses for
  !                           the n new roots which are to be computed.
  !                           if the user does not provide an initial
  !                           guess, zero is used.
  !                           on output, x(kn+1),...,x(kn+n) contain the
  !                           approximate roots found by muller.
  !                itmax  - the maximum allowable number of iterations
  !                           per root (input)
  !                infer  - an integer vector of length kn+n.  on
  !                           output infer(j) contains the number of
  !                           iterations used in finding the j-th root
  !                           when convergence was achieved.  if
  !                           convergence was not obtained in itmax
  !                           iterations, infer(j) will be greater than
  !                           itmax (output).
  !                ier    - error parameter (output)
  !                         warning error
  !                           ier = 33 indicates failure to converge with-
  !                             in itmax iterations for at least one of
  !                             the (n) new roots.
  !
  !
  !   remarks      muller always returns the last approximation for root j
  !                in x(j). if the convergence criterion is satisfied,
  !                then infer(j) is less than or equal to itmax. if the
  !                convergence criterion is not satisified, then infer(j)
  !                is set to either itmax+1 or itmax+k, with k greater
  !                than 1. infer(j) = itmax+1 indicates that muller did
  !                not obtain convergence in the allowed number of iter-
  !                ations. in this case, the user may wish to set itmax
  !                to a larger value. infer(j) = itmax+k means that con-
  !                vergence was obtained (on iteration k) for the defla-
  !                ted function
  !                              fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
  !
  !                but failed for f(z). in this case, better initial
  !                guesses might help or, it might be necessary to relax
  !                the convergence criterion.
  !
  !-----------------------------------------------------------------------
  !
  subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
       implicit none
        double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
        double complex ::   d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
                            tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
                            zero,p1,one,four,p5
        
        double complex, external :: f
        integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
                      knpng,jk,ick,nn,lm1,errcode
        double complex :: x(kn+n)
        integer :: infer(kn+n)
        
        
        data                zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
                            one/(1.0,0.0)/,four/(4.0,0.0)/, &
                            p5/(0.5,0.0)/, &
                            rzero/0.0/,rten/10.0/,rhun/100.0/, &
                            ax/0.1/,ickmax/3/,rp01/0.01/
  
              ier = 0
              if (n .lt. 1) then ! What the hell are doing here then ...
                  return
              end if
              !eps1 = rten **(-nsig)
              eps1 = min(eps1,rp01)
              
              knp1 = kn+1
              knpn = kn+n
              knpng = kn+nguess
              do i=1,knpn
                  infer(i) = 0
                  if (i .gt. knpng) x(i) = zero
              end do
              l= knp1
  
              ic=0
  rloop:      do while (l<=knpn)   ! Main loop over new roots
                  jk = 0
                  ick = 0
                  xl = x(l)
  icloop:         do
                      ic = 0
                      h = ax
                      h = p1*h
                      if (cdabs(xl) .gt. ax) h = p1*xl
  !                                  first three points are
  !                                    xl+h,  xl-h,  xl
                      rt = xl+h
                      call deflated_work(errcode)
                      if (errcode == 1) then
                          exit icloop
                      end if
  
                      z0 = fprt
                      y0 = frt
                      x0 = rt
                      rt = xl-h
                      call deflated_work(errcode)
                      if (errcode == 1) then
                          exit icloop
                      end if
  
                      z1 = fprt
                      y1 = frt
                      h = xl-rt
                      d = h/(rt-x0)
                      rt = xl
  
                      call deflated_work(errcode)
                      if (errcode == 1) then
                          exit icloop
                      end if
  
     
                      z2 = fprt
                      y2 = frt
  !                                  begin main algorithm
   iloop:             do
                          dd = one + d
                          t1 = z0*d*d
                          t2 = z1*dd*dd
                          xx = z2*dd
                          t3 = z2*d
                          bi = t1-t2+xx+t3
                          den = bi*bi-four*(xx*t1-t3*(t2-xx))
  !                                  use denominator of maximum amplitude 
                          t1 = cdsqrt(den)
                          qz = rhun*max(cdabs(bi),cdabs(t1))
                          t2 = bi + t1
                          tpq = cdabs(t2)+qz
                          if (tpq .eq. qz) t2 = zero
                          t3 = bi - t1
                          tpq = cdabs(t3) + qz
                          if (tpq .eq. qz) t3 = zero
                          den = t2
                          qz = cdabs(t3)-cdabs(t2)
                          if (qz .gt. rzero) den = t3
  !                                  test for zero denominator            
                          if (cdabs(den) .eq. rzero) then
                              call trans_rt()
                              call deflated_work(errcode)
                              if (errcode == 1) then
                                  exit icloop
                              end if
                              z2 = fprt
                              y2 = frt
                              cycle iloop
                          end if
  
  
                          d = -xx/den
                          d = d+d
                          h = d*h
                          rt = rt + h
  !                                  check convergence of the first kind  
                          if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
                              if (ic .ne. 0) then
                                  exit icloop
                              end if
                              ic = 1
                              z0 = y1
                              z1 = y2
                              z2 = f(rt)
                              xl = rt
                              ick = ick+1
                              if (ick .le. ickmax) then
                                  cycle iloop 
                              end if
  !                                  warning error, itmax = maximum
                              jk = itmax + jk
                              ier = 33
                          end if
                          if (ic .ne. 0) then
                              cycle icloop
                          end if
                          call deflated_work(errcode)
                          if (errcode == 1) then
                              exit icloop
                          end if
  
                          do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
                              !   take remedial action to induce
                              !   convergence
                              d = d*p5
                              h = h*p5
                              rt = rt-h
                              call deflated_work(errcode)
                              if (errcode == 1) then
                                  exit icloop
                              end if
                          end do
                          z0 = z1
                          z1 = z2
                          z2 = fprt
                          y0 = y1
                          y1 = y2
                          y2 = frt
                      end do iloop
                  end do icloop
          x(l) = rt
          infer(l) = jk
          l = l+1
        end do rloop
        
        contains
          subroutine trans_rt()
             tem = rten*eps1
             if (cdabs(rt) .gt. ax) tem = tem*rt
             rt = rt+tem
             d = (h+tem)*d/h
             h = h+tem
          end subroutine trans_rt
          
          subroutine deflated_work(errcode)
              ! errcode=0 => no errors
              ! errcode=1 => jk>itmax or convergence of second kind achieved
              integer :: errcode,flag
              
              flag=1
      loop1:  do while(flag==1)
                  errcode=0
                  jk = jk+1
                  if (jk .gt. itmax) then
                      ier=33
                      errcode=1
                      return
                  end if
                  frt = f(rt)
                  fprt = frt
                  if (l /= 1) then
                      lm1 = l-1
                      do i=1,lm1
                          tem = rt - x(i)
                          if (cdabs(tem) .eq. rzero) then
                          !if (ic .ne. 0) go to 15 !! ?? possible?
                              call trans_rt()
                              cycle loop1
                          end if
                          fprt = fprt/tem
                      end do
                  end if
                  flag=0
              end do loop1
   
              if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
                  errcode=1
                  return
              end if
              
          end subroutine deflated_work
        
        end subroutine
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  !   Integration
  !
  !   Only double precision coded atm
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  
  subroutine fvn_d_gauss_legendre(n,qx,qw)
  !
  ! This routine compute the n Gauss Legendre abscissas and weights
  ! Adapted from Numerical Recipes routine gauleg
  !
  ! n (in) : number of points
  ! qx(out) : abscissas
  ! qw(out) : weights
  !
  implicit none
  double precision,parameter :: pi=3.141592653589793d0
  integer, intent(in) :: n
  double precision, intent(out) :: qx(n),qw(n)
  
  integer :: m,i,j
  double precision :: z,z1,p1,p2,p3,pp
  m=(n+1)/2
  do i=1,m
      z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
  iloop:  do 
              p1=1.d0
              p2=0.d0
              do j=1,n
                  p3=p2
                  p2=p1
                  p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
              end do
              pp=dble(n)*(z*p1-p2)/(z*z-1.d0)
              z1=z
              z=z1-p1/pp
              if (dabs(z-z1)<=epsilon(z)) then
                  exit iloop
              end if
          end do iloop
      qx(i)=-z
      qx(n+1-i)=z
      qw(i)=2.d0/((1.d0-z*z)*pp*pp)
      qw(n+1-i)=qw(i)
  end do
  end subroutine
  
  
  
  subroutine fvn_d_gl_integ(f,a,b,n,res)
  !
  ! This is a simple non adaptative integration routine 
  ! using n gauss legendre abscissas and weights
  !
  !   f(in)   : the function to integrate
  !   a(in)   : lower bound
  !   b(in)   : higher bound
  !   n(in)   : number of gauss legendre pairs
  !   res(out): the evaluation of the integral
  !
  double precision,external :: f
  double precision, intent(in) :: a,b
  integer, intent(in):: n
  double precision, intent(out) :: res
  
  double precision, allocatable :: qx(:),qw(:)
  double precision :: xm,xr
  integer :: i
  
  ! First compute n gauss legendre abs and weight
  allocate(qx(n))
  allocate(qw(n))
  call fvn_d_gauss_legendre(n,qx,qw)
  
  xm=0.5d0*(b+a)
  xr=0.5d0*(b-a)
  
  res=0.d0
  
  do i=1,n
      res=res+qw(i)*f(xm+xr*qx(i))
  end do
  
  res=xr*res
  
  deallocate(qw)
  deallocate(qx)
  
  end subroutine
  
  !!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Simple and double adaptative Gauss Kronrod integration based on
  ! a modified version of quadpack ( http://www.netlib.org/quadpack
  !
  ! Common parameters :
  !
  !       key (in)
  !       epsabs
  !       epsrel
  !
  !
  !!!!!!!!!!!!!!!!!!!!!!!!
  
  subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
  !
  ! Evaluate the integral of function f(x) between a and b
  !
  ! f(in) : the function
  ! a(in) : lower bound
  ! b(in) : higher bound
  ! epsabs(in) : desired absolute error
  ! epsrel(in) : desired relative error
  ! key(in) : gauss kronrod rule
  !                     1:   7 - 15 points
  !                     2:  10 - 21 points
  !                     3:  15 - 31 points
  !                     4:  20 - 41 points
  !                     5:  25 - 51 points
  !                     6:  30 - 61 points
  !
  ! limit(in) : maximum number of subintervals in the partition of the 
  !               given integration interval (a,b). A value of 500 will give the same
  !               behaviour as the imsl routine dqdag
  !
  ! res(out) : estimated integral value
  ! abserr(out) : estimated absolute error
  ! ier(out) : error flag from quadpack routines
  !               0 : no error
  !               1 : maximum number of subdivisions allowed
  !                   has been achieved. one can allow more
  !                   subdivisions by increasing the value of
  !                   limit (and taking the according dimension
  !                   adjustments into account). however, if
  !                   this yield no improvement it is advised
  !                   to analyze the integrand in order to
  !                   determine the integration difficulaties.
  !                   if the position of a local difficulty can
  !                   be determined (i.e.singularity,
  !                   discontinuity within the interval) one
  !                   will probably gain from splitting up the
  !                   interval at this point and calling the
  !                   integrator on the subranges. if possible,
  !                   an appropriate special-purpose integrator
  !                   should be used which is designed for
  !                   handling the type of difficulty involved.
  !               2 : the occurrence of roundoff error is
  !                   detected, which prevents the requested
  !                   tolerance from being achieved.
  !               3 : extremely bad integrand behaviour occurs
  !                   at some points of the integration
  !                   interval.
  !               6 : the input is invalid, because
  !                   (epsabs.le.0 and
  !                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
  !                   or limit.lt.1 or lenw.lt.limit*4.
  !                   result, abserr, neval, last are set
  !                   to zero.
  !                   except when lenw is invalid, iwork(1),
  !                   work(limit*2+1) and work(limit*3+1) are
  !                   set to zero, work(1) is set to a and
  !                   work(limit+1) to b.
  
  implicit none
  double precision, external :: f
  double precision, intent(in) :: a,b,epsabs,epsrel
  integer, intent(in) :: key
  integer, intent(in) :: limit
  double precision, intent(out) :: res,abserr
  integer, intent(out) :: ier
  
  double precision, allocatable :: work(:)
  integer, allocatable :: iwork(:)
  integer :: lenw,neval,last
  
  ! imsl value for limit is 500
  lenw=limit*4
  
  allocate(iwork(limit))
  allocate(work(lenw))
  
  call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
  
  deallocate(work)
  deallocate(iwork)
  
  end subroutine
  
  
  
  subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
  !
  ! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)
  !
  ! f(in) : the function
  ! a(in) : lower bound
  ! b(in) : higher bound
  ! g(in) : external function describing lower bound for y
  ! h(in) : external function describing higher bound for y
  ! epsabs(in) : desired absolute error
  ! epsrel(in) : desired relative error
  ! key(in) : gauss kronrod rule
  !                     1:   7 - 15 points
  !                     2:  10 - 21 points
  !                     3:  15 - 31 points
  !                     4:  20 - 41 points
  !                     5:  25 - 51 points
  !                     6:  30 - 61 points
  !
  ! limit(in) : maximum number of subintervals in the partition of the 
  !               given integration interval (a,b). A value of 500 will give the same
  !               behaviour as the imsl routine dqdag
  !
  ! res(out) : estimated integral value
  ! abserr(out) : estimated absolute error
  ! ier(out) : error flag from quadpack routines
  !               0 : no error
  !               1 : maximum number of subdivisions allowed
  !                   has been achieved. one can allow more
  !                   subdivisions by increasing the value of
  !                   limit (and taking the according dimension
  !                   adjustments into account). however, if
  !                   this yield no improvement it is advised
  !                   to analyze the integrand in order to
  !                   determine the integration difficulaties.
  !                   if the position of a local difficulty can
  !                   be determined (i.e.singularity,
  !                   discontinuity within the interval) one
  !                   will probably gain from splitting up the
  !                   interval at this point and calling the
  !                   integrator on the subranges. if possible,
  !                   an appropriate special-purpose integrator
  !                   should be used which is designed for
  !                   handling the type of difficulty involved.
  !               2 : the occurrence of roundoff error is
  !                   detected, which prevents the requested
  !                   tolerance from being achieved.
  !               3 : extremely bad integrand behaviour occurs
  !                   at some points of the integration
  !                   interval.
  !               6 : the input is invalid, because
  !                   (epsabs.le.0 and
  !                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
  !                   or limit.lt.1 or lenw.lt.limit*4.
  !                   result, abserr, neval, last are set
  !                   to zero.
  !                   except when lenw is invalid, iwork(1),
  !                   work(limit*2+1) and work(limit*3+1) are
  !                   set to zero, work(1) is set to a and
  !                   work(limit+1) to b.
  
  implicit none
  double precision, external:: f,g,h
  double precision, intent(in) :: a,b,epsabs,epsrel
  integer, intent(in) :: key,limit
  integer, intent(out) :: ier
  double precision, intent(out) :: res,abserr
  
  
  double precision, allocatable :: work(:)
  integer, allocatable :: iwork(:)
  integer :: lenw,neval,last
  
  ! imsl value for limit is 500
  lenw=limit*4
  allocate(work(lenw))
  allocate(iwork(limit))
  
  call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
  
  deallocate(iwork)
  deallocate(work)
  end subroutine
  
  
  
  subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
  !
  ! Evaluate the single integral of function f(x,y) for y between a and b with a
  ! given x value
  !
  ! This function is used for the evaluation of the double integral fvn_d_integ_2_gk
  !
  ! f(in) : the function
  ! x(in) : x
  ! a(in) : lower bound
  ! b(in) : higher bound
  ! epsabs(in) : desired absolute error
  ! epsrel(in) : desired relative error
  ! key(in) : gauss kronrod rule
  !                     1:   7 - 15 points
  !                     2:  10 - 21 points
  !                     3:  15 - 31 points
  !                     4:  20 - 41 points
  !                     5:  25 - 51 points
  !                     6:  30 - 61 points
  !
  ! limit(in) : maximum number of subintervals in the partition of the 
  !               given integration interval (a,b). A value of 500 will give the same
  !               behaviour as the imsl routine dqdag
  !
  ! res(out) : estimated integral value
  ! abserr(out) : estimated absolute error
  ! ier(out) : error flag from quadpack routines
  !               0 : no error
  !               1 : maximum number of subdivisions allowed
  !                   has been achieved. one can allow more
  !                   subdivisions by increasing the value of
  !                   limit (and taking the according dimension
  !                   adjustments into account). however, if
  !                   this yield no improvement it is advised
  !                   to analyze the integrand in order to
  !                   determine the integration difficulaties.
  !                   if the position of a local difficulty can
  !                   be determined (i.e.singularity,
  !                   discontinuity within the interval) one
  !                   will probably gain from splitting up the
  !                   interval at this point and calling the
  !                   integrator on the subranges. if possible,
  !                   an appropriate special-purpose integrator
  !                   should be used which is designed for
  !                   handling the type of difficulty involved.
  !               2 : the occurrence of roundoff error is
  !                   detected, which prevents the requested
  !                   tolerance from being achieved.
  !               3 : extremely bad integrand behaviour occurs
  !                   at some points of the integration
  !                   interval.
  !               6 : the input is invalid, because
  !                   (epsabs.le.0 and
  !                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
  !                   or limit.lt.1 or lenw.lt.limit*4.
  !                   result, abserr, neval, last are set
  !                   to zero.
  !                   except when lenw is invalid, iwork(1),
  !                   work(limit*2+1) and work(limit*3+1) are
  !                   set to zero, work(1) is set to a and
  !                   work(limit+1) to b.
  
  implicit none
  double precision, external:: f
  double precision, intent(in) :: x,a,b,epsabs,epsrel
  integer, intent(in) :: key,limit
  integer, intent(out) :: ier
  double precision, intent(out) :: res,abserr
  
  
  double precision, allocatable :: work(:)
  integer, allocatable :: iwork(:)
  integer :: lenw,neval,last
  
  ! imsl value for limit is 500
  lenw=limit*4
  allocate(work(lenw))
  allocate(iwork(limit))
  
  call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
  
  deallocate(iwork)
  deallocate(work)
  end subroutine
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Include the modified quadpack files
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  include "fvn_quadpack/dqag_2d_inner.f"
  include "fvn_quadpack/dqk15_2d_inner.f"
  include "fvn_quadpack/dqk31_2d_outer.f"
  include "fvn_quadpack/d1mach.f"
  include "fvn_quadpack/dqk31_2d_inner.f"
  include "fvn_quadpack/dqage.f"
  include "fvn_quadpack/dqk15.f"
  include "fvn_quadpack/dqk21.f"
  include "fvn_quadpack/dqk31.f"
  include "fvn_quadpack/dqk41.f"
  include "fvn_quadpack/dqk51.f"
  include "fvn_quadpack/dqk61.f"
  include "fvn_quadpack/dqk41_2d_outer.f"
  include "fvn_quadpack/dqk41_2d_inner.f"
  include "fvn_quadpack/dqag_2d_outer.f"
  include "fvn_quadpack/dqpsrt.f"
  include "fvn_quadpack/dqag.f"
  include "fvn_quadpack/dqage_2d_outer.f"
  include "fvn_quadpack/dqage_2d_inner.f"
  include "fvn_quadpack/dqk51_2d_outer.f"
  include "fvn_quadpack/dqk51_2d_inner.f"
  include "fvn_quadpack/dqk61_2d_outer.f"
  include "fvn_quadpack/dqk21_2d_outer.f"
  include "fvn_quadpack/dqk61_2d_inner.f"
  include "fvn_quadpack/dqk21_2d_inner.f"
  include "fvn_quadpack/dqk15_2d_outer.f"

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Trigonometric functions
  !
  ! fvn_z_acos, fvn_z_asin : complex arc cosine and sine
  ! fvn_d_acosh : arc cosinus hyperbolic
  ! fvn_d_asinh : arc sinus hyperbolic
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! February 2008 
  ! All Trigonometric functions removed due to implementation of fnlib

  function fvn_z_acos(z)
      ! double complex arccos function adapted from 
      ! the c gsl library
      ! http://www.gnu.org/software/gsl/
      implicit none
      complex(kind=8) :: fvn_z_acos
      complex(kind=8) :: z
      real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
      real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
      complex(kind=8),parameter :: i=(0._8,1._8)
      real(kind=8) :: r_res,i_res

      rz=dreal(z)
      iz=dimag(z)

      if ( iz == 0._8 ) then
          fvn_z_acos=fvn_z_acos_real(rz)
          return
      end if
  
      x=dabs(rz)
      y=dabs(iz)
      r=fvn_d_hypot(x+1.,y)
      s=fvn_d_hypot(x-1.,y)
      a=0.5*(r + s)
      b=x/a
      y2=y*y
  
      if (b <= b_crossover) then
          r_res=dacos(b)
      else
          if (x <= 1.) then
              d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
              r_res=datan(dsqrt(d)/x)
          else
              apx=a+x
              d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
              r_res=datan((y*dsqrt(d))/x);
          end if
      end if
  
      if (a <= a_crossover) then
          if (x < 1.) then
              am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
          else
              am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
          end if
          i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
      else
          i_res = dlog (a + dsqrt (a*a - 1.));
      end if
      if (rz <0.) then
          r_res=fvn_pi-r_res
      end if
      i_res=-sign(1._8,iz)*i_res
      fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)
  
  end function fvn_z_acos
  
  function fvn_z_acos_real(r)
      ! return the double complex arc cosinus for a 
      ! double precision argument
      implicit none
      real(kind=8) :: r
      complex(kind=8) :: fvn_z_acos_real
  
      if (dabs(r)<=1._8) then
          fvn_z_acos_real=dcmplx(dacos(r))
          return
      end if
      if (r < 0._8) then
          fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))
      else
          fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))
      end if
  end function
  
  
  function fvn_z_asin(z)
      ! double complex arcsin function derived from 
      ! the c gsl library
      ! http://www.gnu.org/software/gsl/
      implicit none
      complex(kind=8) :: fvn_z_asin
      complex(kind=8) :: z
      real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
      real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
      real(kind=8) :: r_res,i_res

      rz=dreal(z)
      iz=dimag(z)

      if ( iz == 0._8 ) then
          ! z is real
          fvn_z_asin=fvn_z_asin_real(rz)
          return
      end if
      
      x=dabs(rz)
      y=dabs(iz)
      r=fvn_d_hypot(x+1.,y)
      s=fvn_d_hypot(x-1.,y)
      a=0.5*(r + s)
      b=x/a
      y2=y*y
  
      if (b <= b_crossover) then
          r_res=dasin(b)
      else
          if (x <= 1.) then
              d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
              r_res=datan(x/dsqrt(d))
          else
              apx=a+x
              d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
              r_res=datan(x/(y*dsqrt(d)));
          end if
      end if
  
      if (a <= a_crossover) then
          if (x < 1.) then
              am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
          else
              am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
          end if
          i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
      else
          i_res = dlog (a + dsqrt (a*a - 1.));
      end if
      r_res=sign(1._8,rz)*r_res
      i_res=sign(1._8,iz)*i_res
      fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)
  
  end function fvn_z_asin
  
  function fvn_z_asin_real(r)
      ! return the double complex arc sinus for a 
      ! double precision argument
      implicit none
      real(kind=8) :: r
      complex(kind=8) :: fvn_z_asin_real
  
      if (dabs(r)<=1._8) then
          fvn_z_asin_real=dcmplx(dasin(r))
          return
      end if
      if (r < 0._8) then
          fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))
      else
          fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))
      end if
  end function fvn_z_asin_real
  
  function fvn_d_acosh(r)
      ! return the arc hyperbolic cosine
      implicit none
      real(kind=8) :: r
      real(kind=8) :: fvn_d_acosh
      if (r >=1) then
          fvn_d_acosh=dlog(r+dsqrt(r*r-1))
      else
          !! TODO : Better error handling!!!!!!
          stop "Argument to fvn_d_acosh lesser than 1"
      end if
  end function fvn_d_acosh
  
  function fvn_d_asinh(r)
      ! return the arc hyperbolic sine
      implicit none
      real(kind=8) :: r
      real(kind=8) :: fvn_d_asinh
      fvn_d_asinh=dlog(r+dsqrt(r*r+1))
  end function fvn_d_asinh
  
  function fvn_d_hypot(a,b)
      implicit none
      ! return the euclidian norm of vector(a,b)
      real(kind=8) :: a,b
      real(kind=8) :: fvn_d_hypot
      fvn_d_hypot=dsqrt(a*a+b*b)
  end function

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! SPARSE RESOLUTION
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! 
  ! Sparse resolution is done by interfaçing Tim Davi's UMFPACK
  ! http://www.cise.ufl.edu/research/sparse/SuiteSparse/
  ! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig
  !
  ! Solve Ax=B using UMFPACK
  !
  ! Where A is a sparse matrix given in its triplet form 
  ! T -> non zero elements
  ! Ti,Tj -> row and column index (1-based) of the given elt
  ! n : rank of matrix A
  ! nz : number of non zero elts
  !
  ! fvn_*_sparse_solve
  ! * = zl : double complex + integer(8)
  ! * = zi : double complex + integer(4)
  !
  subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
  implicit none
  integer(8), intent(in) :: n,nz
  complex(8),dimension(nz),intent(in) :: T
  integer(8),dimension(nz),intent(in)  :: Ti,Tj
  complex(8),dimension(n),intent(in) :: B
  complex(8),dimension(n),intent(out) :: x
  integer(8), intent(out) :: status
  
  integer(8),dimension(:),allocatable :: wTi,wTj
  real(8),dimension(:),allocatable :: Tx,Tz
  real(8),dimension(:),allocatable :: Ax,Az
  integer(8),dimension(:),allocatable :: Ap,Ai
  integer(8) :: symbolic,numeric
  real(8),dimension(:),allocatable :: xx,xz,bx,bz
  real(8),dimension(90) :: info
  real(8),dimension(20) :: control
  integer(8) :: sys
  
  
  status=0
  
  ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
  ! Tx and Tz are the real and imaginary parts of T
  allocate(wTi(nz),wTj(nz))
  allocate(Tx(nz),Tz(nz))
  Tx=dble(T)
  Tz=aimag(T)
  wTi=Ti-1
  wTj=Tj-1
  allocate(Ax(nz),Az(nz))
  allocate(Ap(n+1),Ai(nz))
  
  ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
  call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
  ! if status is not zero a problem has occured
  if (status /= 0) then
      write(*,*) "Problem during umfpack_zl_triplet_to_col"
  endif
  
  ! Define defaults control values
  call umfpack_zl_defaults(control)
  
  ! Symbolic analysis
  call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during symbolic analysis"
      status=info(1)
  endif
  
  ! Numerical factorization
  call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during numerical factorization"
      status=info(1)
  endif
  
  ! free the C symbolic pointer
  call umfpack_zl_free_symbolic (symbolic)
  
  allocate(bx(n),bz(n),xx(n),xz(n))
  bx=dble(B)
  bz=aimag(B)
  sys=0
  ! sys may be used to define type of solving -> see umfpack.h
  
  ! Solving
  call umfpack_zl_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during solving"
      status=info(1)
  endif
  
  
  ! free the C numeric pointer
  call umfpack_zl_free_numeric (numeric)
  
  x=dcmplx(xx,xz)
  
  deallocate(bx,bz,xx,xz)
  deallocate(Ax,Az)
  deallocate(Tx,Tz)
  deallocate(wTi,wTj)
  end subroutine

  subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
  implicit none
  integer(4), intent(in) :: n,nz
  complex(8),dimension(nz),intent(in) :: T
  integer(4),dimension(nz),intent(in)  :: Ti,Tj
  complex(8),dimension(n),intent(in) :: B
  complex(8),dimension(n),intent(out) :: x
  integer(4), intent(out) :: status
  
  integer(4),dimension(:),allocatable :: wTi,wTj
  real(8),dimension(:),allocatable :: Tx,Tz
  real(8),dimension(:),allocatable :: Ax,Az
  integer(4),dimension(:),allocatable :: Ap,Ai
  !integer(8) :: symbolic,numeric
  integer(4),dimension(2) :: symbolic,numeric
  ! As symbolic and numeric are used to store a C pointer, it is necessary to
  ! still use an integer(8) for 64bits machines
  ! An other possibility : integer(4),dimension(2) :: symbolic,numeric
  real(8),dimension(:),allocatable :: xx,xz,bx,bz
  real(8),dimension(90) :: info
  real(8),dimension(20) :: control
  integer(4) :: sys
  
  status=0
  ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
  ! Tx and Tz are the real and imaginary parts of T
  allocate(wTi(nz),wTj(nz))
  allocate(Tx(nz),Tz(nz))
  Tx=dble(T)
  Tz=aimag(T)
  wTi=Ti-1
  wTj=Tj-1
  allocate(Ax(nz),Az(nz))
  allocate(Ap(n+1),Ai(nz))
  
  ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
  call umfpack_zi_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
  ! if status is not zero a problem has occured
  if (status /= 0) then
      write(*,*) "Problem during umfpack_zl_triplet_to_col"
  endif
  
  ! Define defaults control values
  call umfpack_zi_defaults(control)
  
  ! Symbolic analysis
  call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during symbolic analysis"
      status=info(1)
  endif
  
  ! Numerical factorization
  call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during numerical factorization"
      status=info(1)
  endif
  
  ! free the C symbolic pointer
  call umfpack_zi_free_symbolic (symbolic)
  
  allocate(bx(n),bz(n),xx(n),xz(n))
  bx=dble(B)
  bz=aimag(B)
  sys=0
  ! sys may be used to define type of solving -> see umfpack.h
  
  ! Solving
  call umfpack_zi_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during solving"
      status=info(1)
  endif
  
  ! free the C numeric pointer
  call umfpack_zi_free_numeric (numeric)
  
  x=dcmplx(xx,xz)
  
  deallocate(bx,bz,xx,xz)
  deallocate(Ax,Az)
  deallocate(Tx,Tz)
  deallocate(wTi,wTj)
  end subroutine

  
  
  
  
  subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
  implicit none
  integer(8), intent(in) :: n,nz
  real(8),dimension(nz),intent(in) :: T
  integer(8),dimension(nz),intent(in)  :: Ti,Tj
  real(8),dimension(n),intent(in) :: B
  real(8),dimension(n),intent(out) :: x
  integer(8), intent(out) :: status
  
  integer(8),dimension(:),allocatable :: wTi,wTj
  real(8),dimension(:),allocatable :: A
  integer(8),dimension(:),allocatable :: Ap,Ai
  !integer(8) :: symbolic,numeric
  integer(8) :: symbolic,numeric
  real(8),dimension(90) :: info
  real(8),dimension(20) :: control
  integer(8) :: sys
  
  status=0
  ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
  allocate(wTi(nz),wTj(nz))
  wTi=Ti-1
  wTj=Tj-1
  allocate(A(nz))
  allocate(Ap(n+1),Ai(nz))
  
  ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
  call umfpack_dl_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
  ! if status is not zero a problem has occured
  if (status /= 0) then
      write(*,*) "Problem during umfpack_dl_triplet_to_col"
  endif
  
  ! Define defaults control values
  call umfpack_dl_defaults(control)
  
  ! Symbolic analysis
  call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during symbolic analysis"
      status=info(1)
  endif
  
  ! Numerical factorization
  call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during numerical factorization"
      status=info(1)
  endif
  
  ! free the C symbolic pointer
  call umfpack_dl_free_symbolic (symbolic)
  
  sys=0
  ! sys may be used to define type of solving -> see umfpack.h
  
  ! Solving
  call umfpack_dl_solve (sys, Ap, Ai, A, x, B, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during solving"
      status=info(1)
  endif
  
  ! free the C numeric pointer
  call umfpack_dl_free_numeric (numeric)
  
  deallocate(A)
  deallocate(wTi,wTj)
  end subroutine
  
  
  
  
  
  
  subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
  implicit none
  integer(4), intent(in) :: n,nz
  real(8),dimension(nz),intent(in) :: T
  integer(4),dimension(nz),intent(in)  :: Ti,Tj
  real(8),dimension(n),intent(in) :: B
  real(8),dimension(n),intent(out) :: x
  integer(4), intent(out) :: status
  
  integer(4),dimension(:),allocatable :: wTi,wTj
  real(8),dimension(:),allocatable :: A
  integer(4),dimension(:),allocatable :: Ap,Ai
  !integer(8) :: symbolic,numeric
  integer(4),dimension(2) :: symbolic,numeric
  ! As symbolic and numeric are used to store a C pointer, it is necessary to
  ! still use an integer(8) for 64bits machines
  ! An other possibility : integer(4),dimension(2) :: symbolic,numeric
  real(8),dimension(90) :: info
  real(8),dimension(20) :: control
  integer(4) :: sys
  
  status=0
  ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
  allocate(wTi(nz),wTj(nz))
  wTi=Ti-1
  wTj=Tj-1
  allocate(A(nz))
  allocate(Ap(n+1),Ai(nz))
  
  ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
  call umfpack_di_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status)
  ! if status is not zero a problem has occured
  if (status /= 0) then
      write(*,*) "Problem during umfpack_di_triplet_to_col"
  endif
  
  ! Define defaults control values
  call umfpack_di_defaults(control)
  
  ! Symbolic analysis
  call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during symbolic analysis"
      status=info(1)
  endif
  
  ! Numerical factorization
  call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during numerical factorization"
      status=info(1)
  endif
  
  ! free the C symbolic pointer
  call umfpack_di_free_symbolic (symbolic)
  
  sys=0
  ! sys may be used to define type of solving -> see umfpack.h
  
  ! Solving
  call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info)
  ! info(1) should be zero
  if (info(1) /= 0) then
      write(*,*) "Problem during solving"
      status=info(1)
  endif
  
  ! free the C numeric pointer
  call umfpack_di_free_numeric (numeric)
  
  deallocate(A)
  deallocate(wTi,wTj)
  end subroutine

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Special Functions
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! February 2008 
  ! All Special functions removed due to implementation of fnlib
  !

  function fvn_d_lngamma(x)
        ! This function returns ln(gamma(x))
        ! adapted from Numerical Recipes
        implicit none
        real(kind=8) :: x
        real(kind=8) :: fvn_d_lngamma
  
        real(kind=8) :: ser,stp,tmp,y,cof(6)
        integer(kind=4) :: i
  
        cof = (/ 76.18009172947146d0,-86.50532032941677d0, &
                    24.01409824083091d0,-1.231739572450155d0, &
                    .1208650973866179d-2,-.5395239384953d-5 /)
        stp = 2.5066282746310005d0
  
        tmp=x+5.5d0
        tmp=(x+0.5d0)*log(tmp)-tmp

        ser=1.000000000190015d0

        y=x
  
        do i=1,6
              y=y+1.d0
              ser=ser+cof(i)/y
        end do
        fvn_d_lngamma=tmp+log(stp*ser/x)
  end function
  
  function fvn_d_factorial(n)
        ! This function returns factorial(n) as a real(8)
        ! adapted from Numerical Recipes
        ! real value is calculated for integers lower than 32
        implicit none
        integer(kind=4) :: n
        real(kind=8) :: fvn_d_factorial
  
        integer(kind=4) :: j
  
        fvn_d_factorial=1.
  
        if (n < 0) then
              write(*,*) "Factorial of a negative integer"
              stop
        end if
  
        if (n == 0) then
              return
        end if
  
        if (n <= 32) then 
              do j=1,n
                    fvn_d_factorial=fvn_d_factorial*j
              end do
              return
        else
              fvn_d_factorial=exp(fvn_d_lngamma(dble(n)+1.))
              return
        end if
  end function
  
  function fvn_d_csevl(x,a,n)
        implicit none
        ! This function evaluate the n-term chebyshev series a at x
        ! directly adapted from http://www.netlib.org/fn
        real(kind=8) :: x
        real(kind=8), dimension(n) :: a
        integer(kind=4) :: n
        real(kind=8) :: fvn_d_csevl
  
        real(kind=8) :: twox, b0, b1, b2
        integer(kind=4) :: i,ni
  
        twox = 2.0d0*x
        b1 = 0.d0
        b0 = 0.d0
        do i=1,n
          b2 = b1
          b1 = b0
          ni = n - i + 1
          b0 = twox*b1 - b2 + a(ni)
        end do
  
        fvn_d_csevl = 0.5d0 * (b0-b2)
  
  end function
  
  function fvn_s_csevl(x,a,n)
        implicit none
        ! This function evaluate the n-term chebyshev series a at x
        ! directly adapted from http://www.netlib.org/fn
        real(kind=4) :: x
        real(kind=4), dimension(n) :: a
        integer(kind=4) :: n
        real(kind=4) :: fvn_s_csevl
  
        real(kind=4) :: twox, b0, b1, b2
        integer(kind=4) :: i,ni
  
        twox = 2.0d0*x
        b1 = 0.d0
        b0 = 0.d0
        do i=1,n
          b2 = b1
          b1 = b0
          ni = n - i + 1
          b0 = twox*b1 - b2 + a(ni)
        end do
  
        fvn_s_csevl = 0.5d0 * (b0-b2)
  
  end function

  end module fvn