diff --git a/stable/fvnlib.f90 b/stable/fvnlib.f90
new file mode 100644
index 0000000..c4eceac
--- /dev/null
+++ b/stable/fvnlib.f90
@@ -0,0 +1,1779 @@
+
+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
+! 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 ;-) Nevertheless, you
+! may give credits to quadpack authors. 
+!
+! Version 1.1
+!
+! 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
+! 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
+            
+            
+contains 
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! 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
+    !
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    real, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    real, allocatable :: work(:)
+    real twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double precision, intent(in) :: a(d,d)
+    double precision, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    double precision, allocatable :: work(:)
+    double precision :: twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    complex, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    complex, allocatable :: work(:)
+    complex :: twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double complex, intent(in) :: a(d,d)
+    double complex, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    double complex, allocatable :: work(:)
+    double complex :: twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    real :: fvn_s_det
+    
+    real, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    double precision, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    double precision :: fvn_d_det
+    
+    double precision, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    complex :: fvn_c_det
+    
+    complex, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    double complex, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    double complex :: fvn_z_det
+    
+    double complex, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    real, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double precision, intent(in) :: a(d,d)
+    double precision, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    real, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double complex, intent(in) :: a(d,d)
+    double precision, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    complex, intent(out) :: evala(d)
+    complex, intent(out) :: eveca(d,d)
+    integer, intent(out) :: 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
+    
+    ! 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
+        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
+    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) :: 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
+    
+    ! 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
+        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
+
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    complex, intent(out) :: evala(d)
+    complex, intent(out) :: eveca(d,d)
+    integer, intent(out) :: 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
+
+    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
+        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
+
+    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) :: 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
+    
+    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
+        status=0
+    end if
+    deallocate(rwork)
+    deallocate(work)
+    deallocate(wc_a)
+
+end subroutine
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! 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
+
+
+!
+! 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"
+
+
+
+
+
+end module fvn
\ No newline at end of file
diff --git a/testing/fvnlib.f90 b/testing/fvnlib.f90
new file mode 100644
index 0000000..c4eceac
--- /dev/null
+++ b/testing/fvnlib.f90
@@ -0,0 +1,1779 @@
+
+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
+! 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 ;-) Nevertheless, you
+! may give credits to quadpack authors. 
+!
+! Version 1.1
+!
+! 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
+! 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
+            
+            
+contains 
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! 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
+    !
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    real, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    real, allocatable :: work(:)
+    real twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double precision, intent(in) :: a(d,d)
+    double precision, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    double precision, allocatable :: work(:)
+    double precision :: twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    complex, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    complex, allocatable :: work(:)
+    complex :: twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double complex, intent(in) :: a(d,d)
+    double complex, intent(out) :: inva(d,d)
+    integer, intent(out) :: status
+
+    integer, allocatable :: ipiv(:)
+    double complex, allocatable :: work(:)
+    double complex :: twork(1)
+    integer :: info
+    integer :: lwork
+
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    real :: fvn_s_det
+    
+    real, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    double precision, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    double precision :: fvn_d_det
+    
+    double precision, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    complex :: fvn_c_det
+    
+    complex, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    double complex, intent(in) :: a(d,d)
+    integer, intent(out) :: status
+    double complex :: fvn_z_det
+    
+    double complex, allocatable :: wc_a(:,:)
+    integer, allocatable :: ipiv(:)
+    integer :: info,i
+    
+    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
+        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
+    !
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    real, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double precision, intent(in) :: a(d,d)
+    double precision, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    real, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+    !
+    integer, intent(in) :: d
+    double complex, intent(in) :: a(d,d)
+    double precision, intent(out) :: rcond
+    integer, intent(out) :: 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
+    
+    
+    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
+        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
+        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
+
+    integer, intent(in) :: d
+    real, intent(in) :: a(d,d)
+    complex, intent(out) :: evala(d)
+    complex, intent(out) :: eveca(d,d)
+    integer, intent(out) :: 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
+    
+    ! 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
+        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
+    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) :: 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
+    
+    ! 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
+        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
+
+    integer, intent(in) :: d
+    complex, intent(in) :: a(d,d)
+    complex, intent(out) :: evala(d)
+    complex, intent(out) :: eveca(d,d)
+    integer, intent(out) :: 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
+
+    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
+        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
+
+    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) :: 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
+    
+    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
+        status=0
+    end if
+    deallocate(rwork)
+    deallocate(work)
+    deallocate(wc_a)
+
+end subroutine
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+! 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
+
+
+!
+! 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"
+
+
+
+
+
+end module fvn
\ No newline at end of file