Numerical Methods

Outline

Introduction

Linear
Equations

Interpolation

Numerical
Solutions

Computer
Measurement

      

Fortran library

Fortran has the reputation for delivering the fastest execution speeds and often scientists working on computationally intensive problems use Fortran. Many Fortran libraries containing efficient routines are available. Below is a library of routines that were discussed in this course.

module numeric_lib
  ! Fortran library containing the routines discussed in the Numerical Methods course at the TU Graz
  ! Created by Markus Krammer

  implicit none

  double precision, parameter :: m_pi=4.d0*atan(1.d0) !=3.141592653589793238462643383279502884197169399375105820974944d0

contains

  subroutine jacobi(a,n,tmax,eigvec,error)
    double precision, intent(inout) :: a(:,:),eigvec(:,:)
    integer, intent(in) :: n,tmax
    integer, intent(out) :: error
    ! solves a symmetric real eigenvalue problem a*eigvec=E*eigvec where E are the eigenvalues
    ! a(:,:):      input matrix, eigenvaluses after call in the main diagonal of a(:,:), only upper triangle needed
    ! n:           number of rows / columns of the matrix
    ! tmax:        maximum steps of iteration
    ! eigvec(:,:): eigenvectors in the columns: eigenvector i = eigvec(:,i)
    ! error:       0->convergence reached, 1->iteration failed

    integer :: i1,i2,i3,t
    double precision :: sum_a,limit_a,theta,g,h,tfac,sine,cosine,store_a

    error=1
    if(n<2) then
       eigvec(1,1)=1.0d0
       return
    end if
    eigvec=0.0d0
    do i1=1,n
       eigvec(i1,i1)=1.0d0
    end do

    t=0
    do while(t.le.tmax)
       sum_a=0.0d0
       do i1=1,n-1
          do i2=i1+1,n
             sum_a=sum_a+2.0d0*a(i1,i2)**2
          end do
       end do
       if(sum_a.eq.0.0d0) then
          error=0
          exit
       end if
       limit_a=sqrt(sum_a)/n**2
       do i1=1,n-1
          do i2=i1+1,n
             g=1.0d2*abs(a(i1,i2))
             if(abs(a(i1,i2)).gt.limit_a) then
                h=a(i1,i1)-a(i2,i2)
                if((abs(h)+g).eq.abs(h)) then
                   tfac=a(i1,i2)/h
                else
                   theta=h/(2.0d0*a(i1,i2))
                   tfac=1.0d0/(abs(theta)+sqrt(1.0d0+theta**2))
                   if(theta.lt.0.0d0) then
                      tfac=-tfac
                   end if
                end if
                cosine=1.0d0/sqrt(tfac**2+1.0d0)
                sine=tfac*cosine
                do i3=i1+1,i2-1
                   store_a=a(i1,i3)
                   a(i1,i3)=cosine*a(i1,i3)+sine*a(i3,i2)
                   a(i3,i2)=cosine*a(i3,i2)-sine*store_a
                end do
                do i3=i2+1,n
                   store_a=a(i1,i3)
                   a(i1,i3)=cosine*a(i1,i3)+sine*a(i2,i3)
                   a(i2,i3)=cosine*a(i2,i3)-sine*store_a
                end do
                do i3=1,i1-1
                   store_a=a(i3,i1)
                   a(i3,i1)=cosine*a(i3,i1)+sine*a(i3,i2)
                   a(i3,i2)=cosine*a(i3,i2)-sine*store_a
                end do
                store_a=a(i1,i1)
                a(i1,i1)=cosine**2*a(i1,i1)+2.0d0*cosine*sine*a(i1,i2)+sine**2*a(i2,i2)
                a(i2,i2)=cosine**2*a(i2,i2)-2.0d0*cosine*sine*a(i1,i2)+sine**2*store_a
                a(i1,i2)=0.0d0
                do i3=1,n
                   store_a=eigvec(i3,i2)
                   eigvec(i3,i2)=cosine*eigvec(i3,i2)-sine*eigvec(i3,i1)
                   eigvec(i3,i1)=cosine*eigvec(i3,i1)+sine*store_a
                end do
             end if
          end do
       end do
       t=t+1
    end do
    if(error.ne.0) then
       print '(a,i0,a)','ERROR in jacobi: convergence not reached after ',t,' iterations'
    end if

  end subroutine jacobi

  subroutine tql2(n,d,e,z,error)
    integer, intent(in) :: n
    double precision, intent(inout) :: d(:),e(:),z(:,:)
    integer, intent(out) :: error
    ! calculate the eigenvalues and eigenvectors of a symmetric tridiagonal matrix by the ql method,
    ! a full symmetric matrix can also be diagonalized if tred2(...) is used to reduce the problem
    ! to a tridiagonal problem
    ! n:      number of rows / columns in the matrix
    ! d(:):   conatins the diagonal matrix elements as input and the ascending sorted eigenvalues as output
    ! e(:):   contains the subdiagonal elements, e(1) is arbitraty and e(2:n) holds the (n-1) values,
    !         e is destroyed during the calculations
    ! z(:,:): holds the transformation matrix produced by tred2(:) as input and the orthonormal
    !         eigenvectors as output, sorting is the same as for the eigenvalues,
    !         if the problem is tridiagonal, z has to be the unity matrix as input
    ! error:  0->convergence reached, i>0->iteration failed during the calculation of the i-th eigenvalue,
    !                                      eigenvalues and -vectors 1 to i-1 are correct but not sorted

    integer i1,i2,i3,i4,i5,i41,i42
    integer, allocatable :: indx(:)
    double precision c1,c2,c3,di41,ei41,f,g,h,p,r,s1,s2,tst1,tst2
    integer, parameter :: maxtry=30
    logical :: conv

    if(n.eq.1) then
       error=0
       return
    end if
    e(1:(n-1))=e(2:n)
    f=0.0d0
    tst1=0.0d0
    e(n)=0.0d0
    do i4=1,n
       error=i4
       h=abs(d(i4))+abs(e(i4))
       if(tst1.lt.h) tst1=h
       do i5=i4,n
          tst2=tst1+abs(e(i5))
          if(tst2.eq.tst1) exit
       end do
       if(i5.ne.i4) then
          conv=.false.
          do i2=1,maxtry
             i41=i4+1
             i42=i41+1
             g=d(i4)
             p=(d(i41)-g)/(2.0d0*e(i4))
             r=sign(sqrt(p**2+1.0d0),p)+p
             d(i4)=e(i4)/r
             d(i41)=e(i4)*r
             di41=d(i41)
             h=g-d(i4)
             if(i42.le.n) then
                do i1=i42,n
                   d(i1)=d(i1)-h
                end do
             end if
             f=f+h
             p=d(i5)
             c1=1.0d0
             c2=c1
             ei41=e(i41)
             s1=0.0d0
             do i1=i5-1,i4,-1
                c3=c2
                c2=c1
                s2=s1
                g=c1*e(i1)
                h=c1*p
                r=sqrt(p**2+e(i1)**2)
                e(i1+1)=s1*r
                s1=e(i1)/r
                c1=p/r
                p=c1*d(i1)-s1*g
                d(i1+1)=h+s1*(c1*g+s1*d(i1))
                do i3=1,n
                   h=z(i3,i1+1)
                   z(i3,i1+1)=s1*z(i3,i1)+c1*h
                   z(i3,i1)=c1*z(i3,i1)-s1*h
                end do
             end do
             p=-s1*s2*c3*ei41*e(i4)/di41
             e(i4)=s1*p
             d(i4)=c1*p
             tst2=tst1+abs(e(i4))
             if(tst2.le.tst1) then
                conv=.true.
                exit
             end if
          end do
          if(.not.conv) then
             print '(a,i0,a)','ERROR in tql2: only ',error-1,' eigenvalues were determined'
             return
          end if
       end if
       d(i4)=d(i4)+f
    end do
    allocate(indx(n))
    call hpsort_indx(n,d,indx)
    z=z(:,indx)
    deallocate(indx)

  end subroutine tql2

  subroutine tred2(n,a,d,e)
    integer, intent(in) :: n
    double precision, intent(inout) :: a(:,:)
    double precision, intent(out) :: d(:),e(:)
    ! reduce a symmetric real eigenvalue problem a*eigvec=E*eigvec (where E are the eigenvalues)
    ! to a symmetric tridiagonal eigenvalue problem which can be solved by tql2(...) using 
    ! orthogonal similarity transformation
    ! n:      number of rows / columns of the matrix a
    ! a(:,:): symmetric matrix of the eigenvalueproblem as input and transformation matrix z(:,:)
    !         for tql2(...) as output, only the lower triangle has to be provided
    ! d(:):   diagonal elements of the tridiagonal matrix
    ! e(:):   subdiagonal elements of the tridiagonal matrix on sites 2 to n

    integer :: i1,i2,i3,i4,i2p1
    double precision :: f,g,h,hh,scale

    d=a(n,:)
    if(n.eq.1) then
       return
    end if
    do i1=n,2,-1
       i4=i1-1
       h=0.0d0
       scale=0.0d0
       if(i4.ge.2) then
          do i3=1,i4
             scale=scale+abs(d(i3))
          end do
       end if
       if(scale.eq.0.0d0) then
          e(i1)=d(i4)
          do i2=1,i4
             d(i2)=a(i4,i2)
             a(i1,i2)=0.0d0
             a(i2,i1)=0.0d0
          end do
       else
          do i3=1,i4
             d(i3)=d(i3)/scale
             h=h+d(i3)**2
          end do
          f=d(i4)
          g=-sign(sqrt(h),f)
          e(i1)=scale*g
          h=h-f*g
          d(i4)=f-g
          do i2=1,i4
             e(i2)=0.0d0
          end do
          do i2=1,i4
             f=d(i2)
             a(i2,i1)=f
             g=e(i2)+a(i2,i2)*f
             i2p1=i2+1
             if(i4.ge.i2p1) then
                do i3=i2p1,i4
                   g=g+a(i3,i2)*d(i3)
                   e(i3)=e(i3)+a(i3,i2)*f
                end do
             end if
             e(i2)=g
          end do
          f=0.0d0
          do i2=1,i4
             e(i2)=e(i2)/h
             f=f+e(i2)*d(i2)
          end do
          hh=f/(2.0d0*h)
          do i2=1,i4
             e(i2)=e(i2)-hh*d(i2)
          end do
          do i2=1,i4
             f=d(i2)
             g=e(i2)
             do i3=i2,i4
                a(i3,i2)=a(i3,i2)-f*e(i3)-g*d(i3)
             end do
             d(i2)=a(i4,i2)
             a(i1,i2)=0.0d0
          end do
       end if
       d(i1)=h
    end do
    do i1=2,n
       i4=i1-1
       a(n,i4)=a(i4,i4)
       a(i4,i4)=1.0d0
       h=d(i1)
       if(h.ne.0.0d0) then
          do i3=1,i4
             d(i3)=a(i3,i1)/h
          end do
          do i2=1,i4
             g=0.0d0
             do i3=1,i4
                g=g+a(i3,i1)*a(i3,i2)
             end do
             do i3=1,i4
                a(i3,i2)=a(i3,i2)-g*d(i3)
             end do
          end do
       end if
       do i3=1,i4
          a(i3,i1)=0.0d0
       end do
    end do
    do i1=1,n
       d(i1)=a(n,i1)
       a(n,i1)=0.0d0
    end do
    a(n,n)=1.0d0
    e(1)=0.0d0

  end subroutine tred2

  subroutine tdiag(n,a,eigval)
    integer, intent(in) :: n
    double precision, intent(inout) :: a(:,:)
    double precision, intent(out) :: eigval(:)
    ! solve real symmetric eigenvalue problem a*x=eigval*x with tred2(...) and tql2(...)
    ! n:         number of rows / columns of the matrix a
    ! a(:,:):    real symmetric matrix of the eigenvalue problem as input and eigenvectors as output
    !            eigenvectors are found in the columns: eigenvector i = a(:,i)
    ! eigval(:): eigenvalues of the initial matrix a

    double precision, allocatable :: e(:)
    integer :: error

    allocate(e(n))
    call tred2(n,a,eigval,e)
    call tql2(n,eigval,e,a,error)
    deallocate(e)

  end subroutine tdiag

  subroutine reduc(n,a,b,error)
    integer, intent(in) :: n
    double precision, intent(inout) :: a(:,:),b(:,:)
    integer, intent(out) :: error
    ! reduce the generalized symmetric eigenvalue problem a*x=E*b*x with eigenvectors x and eigenvalues E
    ! to a standard eigenvalue problem using cholesky factorization of matrix b
    ! n:      number of rows / columns of matrix a and b
    ! a(:,:): input:  real symmetric matrix a(:,:) (only full upper triangle needs to be provided)
    !         output: symmetric matrix of the standard eigenvalue problem in the full lower triangle, the 
    !                 strict upper triangle is not changed
    ! b(:,:): input:  real symmetric and positiv definite matrix b(:,:) (only full upper triangle needs to
    !                 be provided)
    !         output: cholesky factor l (symmetric matrix) in the full lower diagonal, the strict upper
    !                 triangle is not changed
    ! error:  0->reduction was successful, 1->b(:,:) is not positive definite

    integer :: i1,i2,i3,i1m1,i2m1
    double precision :: x,y
    double precision, allocatable :: dl(:)

    allocate(dl(n))
    error=0
    do i1=1,n
       i1m1=i1-1
       do i2=i1,n
          x=b(i1,i2)
          if(i1.ne.1) then
             do i3=1,i1m1
                x=x-b(i1,i3)*b(i2,i3)
             end do
          end if
          if(i2.eq.i1) then
             if(x.le.0.0d0) then
                print '(a)','ERROR in reduc: b is not positive definite'
                error=1
                return
             end if
             y=sqrt(x)
             dl(i1)=y
          else
             b(i2,i1)=x/y
          end if
       end do
    end do
    do i1=1,n
       i1m1=i1-1
       y=dl(i1)
       do i2=i1,n
          x=a(i1,i2)
          if(i1.gt.1) then
             do i3=1,i1m1
                x=x-b(i1,i3)*a(i2,i3)
             end do
          end if
          a(i2,i1)=x/y
       end do
    end do
    do i2=1,n
       i2m1=i2-1
       do i1=i2,n
          x=a(i1,i2)
          if(i1.ne.i2) then
             do i3=i2,i1-1
                x=x-a(i3,i2)*b(i1,i3)
             end do
          end if
          if(i2.gt.1) then
             do i3=1,i2m1
                x=x-a(i2,i3)*b(i1,i3)
             end do
          end if
          a(i1,i2)=x/dl(i1)
       end do
    end do
    do i1=1,n
       b(i1,i1)=dl(i1)
    end do
    deallocate(dl)

  end subroutine reduc

  subroutine eigsys(n,a,b,eigval)
    integer, intent(in) :: n
    double precision, intent(inout) :: a(:,:),b(:,:)
    double precision, intent(out) :: eigval(:)
    ! solves a generalized eigenvalue problem a*x=eigval*b*x with eigenvectors x
    ! n:         number of rows / columns of the matricies a(:,:) and b(:,:)
    ! a(:,:):    input:  real symmetric matrix a(:,:)
    !            output: eigenvectors in the columns: eigenvector i = a(:,i)
    ! b(:,:):    real symmetric positiv defined matrix b(:,:) as input
    !            values modified during diagonalization
    ! eigval(:): eigenvalues sorted corresponding to the eigenvectors

    integer :: i1,i2,i3,error
    double precision :: sum

    call reduc(n,a,b,error)
    if(error.ne.0) return
    call tdiag(n,a,eigval)
    do i3=1,n
       a(n,i3)=a(n,i3)/b(n,n)
       do i1=n-1,1,-1
          sum=a(i1,i3)
          do i2=i1+1,n
             sum=sum-b(i2,i1)*a(i2,i3)
          end do
          a(i1,i3)=sum/b(i1,i1)
       end do
    end do

  end subroutine eigsys

  subroutine gausei(n,ndiag,s,diag,f,tmax,irel,acc,sol,error,w,t0,t_out)
    integer, intent(in) :: n,ndiag,s(:),tmax,irel
    double precision, intent(in) :: diag(:,:),f(:),acc
    double precision, intent(out) :: sol(:)
    integer, intent(out) :: error
    double precision, optional, intent(inout) :: w
    integer, optional, intent(in) :: t0
    integer, optional, intent(out) :: t_out
    ! solves the equation system diag*sol=f for sparse matricies
    ! n:         number of rows / columns of the matrix
    ! ndiag:     number of filled diagonals
    ! s(:):      position of diagonals, 0 is the main diagonal, +/-i is the i-th upper/lower secondary diagonal
    ! diag(:,:): matrix stored in sparse form
    ! f(:):      inhomogeniety vector
    ! tmax:      maximum number of iterations
    ! irel:      1->relative accuracy, else->absolute accuracy
    ! acc:       desired accuracy
    ! sol(:):    solution vector
    ! error:     0->convergence reached, 1->iteration failed
    ! w:         optimizationparameter omega
    ! t0:        step where an omega optimization should take place
    ! t_out:     iterations needed to finish

    integer ::  main_diag,t,i1,i2,i3,t0_intern
    double precision :: s1,dx,sum_dx,old_sum_dx,err_dx,w_intern

    if(present(w)) then
       w_intern=w
    else
       w_intern=1.0d0
    end if
    if(present(t0)) then
       t0_intern=t0
       w_intern=1.0d0
    else
       t0_intern=tmax+2
    end if

    error=1
    main_diag=0
    do i1=1,ndiag
       if(s(i1).eq.0) then
          main_diag=i1
       end if
    end do
    if(main_diag.eq.0) then
       print '(a)','FATAL ERROR in gausei: no main diagonal found'
       return
    end if
    do i1=1,n
       if(diag(i1,main_diag).eq.0.0d0) then
          print '(a)','FATAL ERROR in gausei: main diagonal contains zeros'
          return
       end if
    end do
    sol=0.0d0

    old_sum_dx=0.0d0
    sum_dx=0.0d0
    t=0
    do while((t.lt.tmax).and.(error.eq.1))
       error=0
       do i1=1,n
          s1=0.0d0
          do i2=1,ndiag
             i3=s(i2)+i1
             if((i2.ne.main_diag).and.(i3.le.n)) then
                s1=s1+diag(i1,i2)*sol(i3)
             end if
          end do
          dx=w_intern*((s1-f(i1))/diag(i1,main_diag)+sol(i1))
          if(t.eq.t0_intern-1) then
             old_sum_dx=old_sum_dx+dx**2  
          else if(t.eq.t0_intern) then
             sum_dx=sum_dx+dx**2
          end if
          sol(i1)=sol(i1)-dx
          if(irel.eq.1) then
             err_dx=abs(dx/sol(i1))
          else
             err_dx=abs(dx)
          end if
          if(err_dx.gt.acc) then
             error=1
          end if
       end do
       if(t.eq.t0_intern) then
          w_intern=2.0d0/(1.0d0+sqrt(1.0d0-sqrt(sum_dx/old_sum_dx)))
       end if
       t=t+1
    end do
    if(error.eq.1) then
       print '(a,i0,a)','ERROR in gausei: convergence not reached after ',tmax,' iterations'
    end if
    if(present(t_out)) then
       t_out=t
    end if
    if(present(w)) then
       w=w_intern
    end if

  end subroutine gausei

  subroutine ludcmp(a,n,indx,d,khad)
    double precision, intent(inout) :: a(:,:)
    integer, intent(in) :: n
    integer, intent(out) :: indx(:)
    double precision, intent(out) :: d,khad
    ! perform a LU decomposition of a real quadratic matrix
    ! a(:,:):  input matrix and LU decomposited output matrix
    ! n:       number of rows / columns of the matrix
    ! indx(:): index of the swaps lines
    ! d:       calculated determinant
    ! khad:    calculated Hadramard condition number

    integer :: i1,i2,i3,imax
    double precision :: sum,aamax,dum
    double precision, parameter :: ctiny=1.0d-20

    khad=1.d0
    do i1=1,n
       sum=0.0d0
       do i2=1,n
          sum=sum+a(i1,i2)**2
       end do
       khad=khad*sqrt(sum)
    end do

    d=1.d0
    do i2=1,n
       do i1=1,i2-1
          sum=a(i1,i2)
          do i3=1,i1-1
             sum=sum-a(i1,i3)*a(i3,i2)
          end do
          a(i1,i2)=sum
       end do
       aamax=0.0d0
       do i1=i2,n
          sum=a(i1,i2)
          do i3=1,i2-1
             sum=sum-a(i1,i3)*a(i3,i2)
          end do
          a(i1,i2)=sum
          dum=abs(sum)
          if(dum.ge.aamax) then
             imax=i1
             aamax=dum
          end if
       end do
       if(i2.ne.imax) then
          do i3=1,n
             dum=a(imax,i3)
             a(imax,i3)=a(i2,i3)
             a(i2,i3)=dum
          end do
          d=-d
       end if
       indx(i2)=imax
       if(a(i2,i2).eq.0.0d0) then
          a(i2,i2)=ctiny
       end if
       dum=1.0d0/a(i2,i2)
       do i1=i2+1,n
          a(i1,i2)=a(i1,i2)*dum
       end do
       d=d*a(i2,i2)
    end do

    khad=abs(d)/khad

  end subroutine ludcmp

  subroutine lubksb(a,n,indx,b)
    double precision, intent(in) :: a(:,:)
    integer, intent(in) :: n,indx(:)
    double precision, intent(inout) :: b(:)
    ! calculate the solution vector a*x=b for a LU decomposed matrix a
    ! a(:,:):  LU decomposited matrix a
    ! n:       number of rows / columns of the matrix
    ! indx(:): index of the swaped lines
    ! b(:):    inhomogeniety vector b as input and solutionvector x as output

    integer :: i1,i2,ii,ll
    double precision :: sum

    ii=0
    do i1=1,n
       ll=indx(i1)
       sum=b(ll)
       b(ll)=b(i1)
       if(ii.ne.0) then
          do i2=ii,i1-1
             sum=sum-a(i1,i2)*b(i2)
          end do
       else if(sum.ne.0.0d0) then
          ii=i1
       end if
       b(i1)=sum
    end do
    do i1=n,1,-1
       sum=b(i1)
       do i2=i1+1,n
          sum=sum-a(i1,i2)*b(i2)
       end do
       b(i1)=sum/a(i1,i1)
    end do

  end subroutine lubksb

  subroutine luinv(a,n,det,khad,woff)
    double precision, intent(inout) :: a(:,:)
    integer, intent(in) :: n
    double precision, intent(out) :: det,khad
    logical, optional, intent(in) :: woff
    ! invert a real quadratic matrix a with LU decomposition
    ! a(:,:): input matrix and inverted output matrix
    ! n:      number of rows / columns of the matrix
    ! det:    derterminant of the matris
    ! khad:   Hadramard condition number of the matrix
    ! woff:   .true.->warnings are turned off, .false. or not present->warnings are displayed

    integer :: i1
    integer, allocatable :: indx(:)
    double precision, allocatable :: mat_LU(:,:),b(:)

    allocate(indx(n))
    allocate(mat_LU(n,n))
    allocate(b(n))
    mat_LU=a
    call ludcmp(mat_LU,n,indx,det,khad)
    if((.not.present(woff)).or.(.not.woff)) then
       if(khad.lt.1.0d-2) then
          print '(a,f12.8)','WARNING in luinv: Hadramard condition number is too low: ',khad
       else if(khad.lt.1.0d-1) then
          print '(a,f12.8)','WARNING in luinv: Hadramard condition number is low: ',khad
       end if
    end if
    do i1=1,n
       b=0.0d0
       b(i1)=1.0d0
       call lubksb(mat_LU,n,indx,b)
       a(:,i1)=b
    end do
    deallocate(indx)
    deallocate(mat_LU)
    deallocate(b)

  end subroutine luinv

  subroutine lusol(a,n,b,x,woff)
    double precision, intent(in) :: a(:,:),b(:)
    integer, intent(in) :: n
    double precision, intent(out) :: x(:)
    logical, optional, intent(in) :: woff
    ! solve the equation system a*x=b of a real quadratic matrix a with LU decomposition
    ! a(:,:): input matrix
    ! n:      number of rows / columns of the matrix
    ! b(:):   inhomogeniety vector
    ! x(:):   solution vector
    ! woff:   .true.->warnings are turned off, .false. or not present->warnings are displayed

    integer, allocatable :: indx(:)
    double precision :: det,khad
    double precision, allocatable :: mat_LU(:,:)

    allocate(indx(n))
    allocate(mat_LU(n,n))
    mat_LU=a
    call ludcmp(mat_LU,n,indx,det,khad)
    if((.not.present(woff)).or.(.not.woff)) then
       if(khad.lt.1.0d-2) then
          print '(a,f12.8)','WARNING in lusol: Hadramard condition number is too low: ',khad
       else if(khad.lt.1.0d-1) then
          print '(a,f12.8)','WARNING in lusol: Hadramard condition number is low: ',khad
       end if
    end if
    x=b
    call lubksb(mat_LU,n,indx,x)
    deallocate(indx)
    deallocate(mat_LU)

  end subroutine lusol

  subroutine trid(b,a,c,r,n)
    double precision, intent(inout) :: b(:),r(:)
    double precision, intent(in) :: a(:),c(:)
    integer, intent(in) :: n
    ! solving an equation system of a tridiagonal matrix 
    ! b(:): main diaglnal of the tridaigonal matrix, changed during process
    ! a(:): first lower minor diagonal stored in a(2:n)
    ! c(:): first upper minor diagonal stored in c(1:(n-1))
    ! r(:): inhomogeniety vector as input and solution vector as output
    ! n:    number of rows / columns of the matrix

    double precision :: m
    integer :: i1

    if(b(1).eq.0.0d0) then
       print '(a)','ERROR in trid: matrix is singular'
       return
    end if
    if(n.lt.2) then
       r(1)=r(1)/b(1)
       return
    end if
    do i1=2,n
       m=a(i1)/b(i1-1)
       b(i1)=b(i1)-m*c(i1-1)
       if(b(i1).eq.0.0d0) then
          print '(a)','ERROR in trid: matrix is singular or close to singularity'
          return
       end if
       r(i1)=r(i1)-m*r(i1-1)
    end do
    r(n)=r(n)/b(n)
    do i1=(n-1),1,-1
       r(i1)=(r(i1)-c(i1)*r(i1+1))/b(i1)
    end do

  end subroutine trid

  subroutine hpsort(n,ra)
    integer, intent(in) :: n
    double precision, intent(inout) :: ra(:)
    ! sort data in vector ra with heapsort algorithm and return the sorted data in vector ra
    ! n:     number of elements in vector ra
    ! ra(:): vector to be sorted

    integer :: i1,i2,imax,isort
    double precision  :: rra

    if(n.lt.2) then
       return
    end if
    isort=n/2+1
    imax=n
    do
       if(isort.gt.1) then
          isort=isort-1
          rra=ra(isort)
       else
          rra=ra(imax)
          ra(imax)=ra(1)
          imax=imax-1
          if(imax.eq.1) then
             ra(1)=rra
             exit
          end if
       end if
       i1=isort
       i2=isort*2
       do while(i2.le.imax)
          if((i2.lt.imax).and.(ra(i2).lt.ra(i2+1))) then
             i2=i2+1
          end if
          if(rra.lt.ra(i2)) then
             ra(i1)=ra(i2)
             i1=i2
             i2=i2*2
          else
             i2=imax+1
          end if
       end do
       ra(i1)=rra
    end do

  end subroutine hpsort

  subroutine hpsort_indx(n,ra,indx)
    integer, intent(in) :: n
    double precision, intent(inout) :: ra(:)
    integer, intent(out) :: indx(:)
    ! sort data in vector ra with heapsort algorithm and return the sorted data in vector ra
    ! n:       number of elements in vector ra
    ! ra(:):   vector to be sorted
    ! indx(:): index where the sorted elements were found in the old array

    integer :: i1,i2,imax,isort,iindx
    double precision  :: rra
    
    do i1=1,n
       indx(i1)=i1
    end do
    if(n.lt.2) then
       return
    end if
    isort=n/2+1
    imax=n
    do
       if(isort.gt.1) then
          isort=isort-1
          rra=ra(isort)
          iindx=indx(isort)
       else
          rra=ra(imax)
          ra(imax)=ra(1)
          iindx=indx(imax)
          indx(imax)=indx(1)
          imax=imax-1
          if(imax.eq.1) then
             ra(1)=rra
             indx(1)=iindx
             exit
          end if
       end if
       i1=isort
       i2=isort*2
       do while(i2.le.imax)
          if((i2.lt.imax).and.(ra(i2).lt.ra(i2+1))) then
             i2=i2+1
          end if
          if(rra.lt.ra(i2)) then
             ra(i1)=ra(i2)
             indx(i1)=indx(i2)
             i1=i2
             i2=i2*2
          else
             i2=imax+1
          end if
       end do
       ra(i1)=rra
       indx(i1)=iindx
    end do

  end subroutine hpsort_indx

  subroutine hpsort_indx_only(n,ra_in,indx)
    integer, intent(in) :: n
    double precision, intent(in) :: ra_in(:)
    integer, intent(out) :: indx(:)
    ! sort data in vector ra with heapsort algorithm and return the indices of sorting in indx
    ! n:        number of elements in vector ra
    ! ra_in(:): vector to be sorted
    ! indx(:):  index where the sorted elements were found in the old array

    integer :: i1,i2,imax,isort,iindx
    double precision  :: rra,ra(n)
    
    do i1=1,n
       ra(i1)=ra_in(i1)
       indx(i1)=i1
    end do
    if(n.lt.2) then
       return
    end if
    isort=n/2+1
    imax=n
    do
       if(isort.gt.1) then
          isort=isort-1
          rra=ra(isort)
          iindx=indx(isort)
       else
          rra=ra(imax)
          ra(imax)=ra(1)
          iindx=indx(imax)
          indx(imax)=indx(1)
          imax=imax-1
          if(imax.eq.1) then
             ra(1)=rra
             indx(1)=iindx
             exit
          end if
       end if
       i1=isort
       i2=isort*2
       do while(i2.le.imax)
          if((i2.lt.imax).and.(ra(i2).lt.ra(i2+1))) then
             i2=i2+1
          end if
          if(rra.lt.ra(i2)) then
             ra(i1)=ra(i2)
             indx(i1)=indx(i2)
             i1=i2
             i2=i2*2
          else
             i2=imax+1
          end if
       end do
       ra(i1)=rra
       indx(i1)=iindx
    end do

  end subroutine hpsort_indx_only

  subroutine nestint_sing(fct,x0,x1,acc_x,acc_fct,xzero)
    interface
       double precision function fct(x)
         double precision, intent(in) :: x
       end function fct
    end interface
    double precision, intent(in) :: x0,x1
    double precision, intent(in) :: acc_x,acc_fct
    double precision, intent(out) :: xzero
    ! find one root of a function fct in the interval [x0,x1] with nested intervals
    ! fct:     pure function declared by the input
    ! x0:      lower boundary of the search interval
    ! x1:      higher boundary of the search interval
    ! acc_x:   desired absolute accuracy for xzero
    ! acc_fct: required minimum absolute value of the function
    ! xzero:   found root

    double precision :: fct0,fct1,x00,x11

    if(x0.ge.x1) then
       print '(a)','ERROR in nestint_sing: x0 is not less than x1'
       xzero=x0
       return
    end if
    fct0=fct(x0)
    fct1=fct(x1)
    if((fct0*fct1).gt.0.0d0) then
       print '(a)','ERROR in nestint_sing: no root in the given interval'
       xzero=x0
       return
    end if
    if((fct0.eq.0.0d0).and.(fct1.eq.0.0d0)) then
       print '(a)','WARNING in nestint_sing: the function value on both sides of the search interval are zero'
       xzero=(x0+x1)/2.0d0
       return
    end if
    x00=x0
    x11=x1
    do
       xzero=(x00+x11)/2.0d0
       fct1=fct(xzero)
       if((abs(fct1).lt.acc_fct).and.((x11-x00).lt.acc_x)) then
          xzero=(x00+x11)/2.0d0
          exit
       end if
       if((fct0*fct1).le.0.0d0) then
          x11=xzero
       else
          fct0=fct1
          x00=xzero
       end if
    end do

  end subroutine nestint_sing

  subroutine nestint_mult(fct,xstart,xend,step,acc_x,acc_fct,nzerosmax,xzeros,nzeros)
    interface
       double precision function fct(x)
         double precision, intent(in) :: x
       end function fct
    end interface
    double precision, intent(in) :: xstart,xend,step,acc_x,acc_fct
    integer, intent(in) :: nzerosmax
    double precision, intent(out) :: xzeros(:)
    integer, intent(out) :: nzeros
    ! find roots of a function fct in the interval [x0,x1] with nested intervals
    ! fct:       pure function declared by the input
    ! xstart:    lower boundary of the search interval
    ! xend:      higher boundary of the search interval
    ! step:      stepsize for scanning the search interval
    ! acc_x:     desired absolute accuracy for xzero
    ! acc_fct:   required minimum absolute value of the function
    ! nzerosmax: maximum number of roots that can be found
    ! xzeros(:): found roots
    ! nzeros:    found number of roots

    double precision :: x0,x1,fct0,fct1

    x0=xstart
    fct0=fct(x0)
    nzeros=0
    do while(x0.lt.xend)
       x1=x0+step
       fct1=fct(x1)
       if((fct0*fct1).le.0.0d0) then
          nzeros=nzeros+1
          call nestint_sing(fct,x0,x1,acc_x,acc_fct,xzeros(nzeros))
          if(nzeros.ge.nzerosmax) then
             print '(a)','WARNING in nestint_mult: maximum number of roots found'
             exit
          end if
       end if
       x0=x1
       fct0=fct1
    end do
    xzeros(nzeros+1:nzerosmax)=0.0d0

  end subroutine nestint_mult

  subroutine gauleg(x1,x2,x,w,n)
    double precision, intent(in) :: x1,x2
    double precision, intent(out) :: x(:),w(:)
    integer, intent(in) :: n
    ! calculating the roots of the gauss legendre polynome and the weighting factors for a gauss quadrature
    ! x1:   start of the interval
    ! x2:   end of the interval
    ! x(:): gauss legendre points
    ! w(:): weighting factors
    ! n:    number of points to calculate

    integer :: i1,i2,m
    double precision :: z1,z,xm,xl,pp,p1,p2,p3
    double precision, parameter :: m_eps=3.0d-14

    m=(n+1)/2
    xm=0.5d0*(x2+x1)
    xl=0.5d0*(x2-x1)
    do i1=1,m
       z=cos(m_pi*(i1-0.25d0)/(n+0.5d0))
       do
          p1=1.0d0
          p2=0.0d0
          do i2=1,n
             p3=p2
             p2=p1
             p1=((2.0d0*i2-1.0d0)*z*p2-(i2-1.0d0)*p3)/i2
          end do
          pp=n*(z*p1-p2)/(z**2-1.0d0)
          z1=z
          z=z1-p1/pp
          if(abs(z-z1).lt.m_eps) then
             exit
          end if
       end do
       x(i1)=xm-xl*z
       x(n+1-i1)=xm+xl*z
       w(i1)=2.0d0*xl/((1.0d0-z**2)*pp**2)
       w(n+1-i1)=w(i1)
    end do

  end subroutine gauleg

  subroutine gauint(fct,x0,x1,irel,acc,ostart,oinc,omax,intfct,error,woff)
    interface
       double precision function fct(x)
         double precision, intent(in) :: x
       end function fct
    end interface
    double precision, intent(in) :: x0,x1,acc
    integer, intent(in) :: irel,ostart,oinc,omax
    double precision, intent(out) :: intfct
    integer, intent(out) :: error
    logical, optional, intent(in) :: woff
    ! integrate the function fct from x0 to x1 with gauss legendre quadrature
    ! fct:    function to integrate
    ! x0:     start of integration interval
    ! x1:     end of integration interval
    ! irel:   1->relative accuracy, else->absolute accuracy
    ! acc:    desired accuracy
    ! ostart: starting order of the quadrature
    ! oinc:   increment of the order of the quadrature
    ! omax:   maximum order of the quadrature
    ! intfct: value of the integral
    ! error:  0->convergence reached, 1->quadrature failed

    integer :: order,i1
    double precision :: oldintfct,err_acc
    double precision, allocatable :: w(:),z(:)

    error=1
    if(ostart.lt.1) then
       print '(a)','ERROR in gauint: starting order < 1'
       return
    end if
    if(oinc.lt.1) then
       print '(a)','ERROR in gauint: order increment < 1'
       return
    end if
    if(omax.lt.ostart) then
       print '(a)','ERROR in gauint: starting order < maximum order'
       return
    end if
    order=ostart
    do while(order.le.omax)
       allocate(z(order))
       allocate(w(order))
       call gauleg(x0,x1,z,w,order)
       intfct=0.0d0
       do i1=1,order
          intfct=intfct+w(i1)*fct(z(i1))
       end do
       deallocate(z)
       deallocate(w)
       if(order.ne.ostart) then
          if((oldintfct.eq.0.0d0).and.(intfct.eq.0.0d0)) then
             if((.not.present(woff)).or.(.not.woff)) print '(a)','WARNING in gauint: integral is zero for every evaluation'
             error=0
             exit
          end if
          err_acc=abs(oldintfct-intfct)
          if(irel.eq.1) then
             err_acc=err_acc/abs(oldintfct)
          end if
          if(err_acc.lt.acc) then
             error=0
             exit
          end if
       end if
       oldintfct=intfct
       order=order+oinc
    end do
    if(error.ne.0) then
       print '(a)','ERROR in gauint: convergence not reached'
    end if

  end subroutine gauint

  subroutine trapzd(fct,x0,x1,intfct,j)
    interface
       double precision function fct(x)
         double precision, intent(in) :: x
       end function fct
    end interface
    double precision, intent(in) :: x0,x1
    double precision, intent(inout) :: intfct
    integer, intent(in) :: j
    ! integrate the function fct from x0 to x1 with trapez methode, j has to start from 1
    ! fct:    function to integrate
    ! x0:     start of integration interval
    ! x1:     end of integration interval
    ! intfct: value of the integral
    ! j:      2**(j-2) new values are added to the integral, j=1 is the initialization for the margins
    !         of the integration interval, j has to increase by 1 every time it is called

    integer :: i1,nint
    double precision :: del,x,sum

    if(j.le.1) then
       intfct=0.5d0*(x1-x0)*(fct(x0)+fct(x1))
    else
       nint=2**(j-2)
       del=(x1-x0)/dble(nint)
       x=x0+0.5d0*del
       sum=0.0d0
       do i1=1,nint
          sum=sum+fct(x)
          x=x+del
       end do
       intfct=0.5d0*(intfct+sum*(x1-x0)/dble(nint))
    end if

  end subroutine trapzd

  subroutine romb(fct,x0,x1,irel,acc,omin,omax,intfct,error)
    interface
       double precision function fct(x)
         double precision, intent(in) :: x
       end function fct
    end interface
    double precision, intent(in) :: x0,x1,acc
    integer, intent(in) :: irel,omin,omax
    double precision, intent(out) :: intfct
    integer, intent(out) :: error
    ! integrate the function fct from x0 to x1 with gauss legendre quadrature
    ! fct:    function to integrate
    ! x0:     start of integration interval
    ! x1:     end of integration interval
    ! irel:   1->relative accuracy, else->absolute accuracy
    ! acc:    desired accuracy
    ! omin:   order where the error calculation starts
    ! omax:   maximum order of the quadrature
    ! intfct: value of the integral
    ! error:  0->convergence reached, 1->quadrature failed

    integer :: i1,ocurr,omax_intern
    double precision :: fac,err_acc
    double precision, allocatable :: quad(:),quadold(:)

    if(omax.lt.3) then
       omax_intern=3
    else
       omax_intern=omax
    end if
    allocate(quad(2:omax_intern))
    allocate(quadold(2:omax_intern))
    error=1
    call trapzd(fct,x0,x1,quad(2),1)
    do ocurr=3,omax_intern
       quadold(2:(ocurr-1))=quad(2:(ocurr-1))
       call trapzd(fct,x0,x1,quad(2),ocurr-1)
       fac=1.0d0
       do i1=2,ocurr-1
          fac=fac*ocurr
          quad(i1+1)=(dble(fac)*quad(i1)-quadold(i1))/(dble(fac-1))
       end do
       if(ocurr.ge.omin) then
          err_acc=abs(quadold(ocurr-1)-quad(ocurr))
          if(irel.eq.1) err_acc=err_acc/abs(quadold(ocurr-1))
          if(err_acc.lt.acc) then
             intfct=quad(ocurr)
             error=0
             exit
          end if
       end if
    end do
    if(error.ne.0) then
       print '(a)','ERROR in romb: convergence not reached'
    end if
    deallocate(quad)
    deallocate(quadold)

  end subroutine romb

  subroutine mrqcof(fct,x,y,sig,ndata,na,a,alpha,beta,chisq)
    interface
       subroutine fct(x,a,y,dyda)
         double precision, intent(in) :: x,a(:)
         double precision, intent(out) :: y,dyda(:)
       end subroutine fct
    end interface
    double precision, intent(in) :: x(:),y(:),sig(:),a(:)
    integer, intent(in) :: ndata,na
    double precision, intent(out) :: alpha(:,:),beta(:),chisq
    ! calculate the Marquardt coefficients 
    ! fct:        model function provided by the user
    !             x:       x value where the function should be evaluated
    !             a(:):    model parameters
    !             y:       function value
    !             dyda(:): derivatives of the function with respect to the model parameters
    ! x(:):       measured x data
    ! y(:):       measured y data
    ! sig(:):     uncertainty of the measured y data
    ! ndata:      number of data points
    ! na:         number of fit parameters
    ! a(:):       guessed fitparameters as input and improved fit parameters as output
    ! alpha(:,:): matrix of the linearized equation system
    ! beta(:):    inhomogeniety vector of the linearized equation system
    ! chisq:      calculated chi-squared

    integer :: i1,i2,i3
    double precision :: ymod,g,dy
    double precision, allocatable :: dyda(:)

    allocate(dyda(na))
    alpha=0.0d0
    beta=0.0d0
    chisq=0.0d0
    do i1=1,ndata
       call fct(x(i1),a,ymod,dyda)
       g=1.0d0/(sig(i1)**2)
       dy=y(i1)-ymod
       do i2=1,na
          do i3=1,i2
             alpha(i2,i3)=alpha(i2,i3)+g*dyda(i2)*dyda(i3)
          end do
          beta(i2)=beta(i2)+g*dy*dyda(i2)
       end do
       chisq=chisq+g*dy**2
    end do
    do i1=2,na
       do i2=1,i1-1
          alpha(i2,i1)=alpha(i1,i2)
       end do
    end do
    deallocate(dyda)

  end subroutine mrqcof

  subroutine mrqmin(fct,x,y,sig,ndata,na,a,alpha,covar,beta,chisq,ochisq,alambda,succmar)
    interface
       subroutine fct(x,a,y,dyda)
         double precision, intent(in) :: x,a(:)
         double precision, intent(out) :: y,dyda(:)
       end subroutine fct
    end interface
    double precision, intent(in) :: x(:),y(:),sig(:)
    integer, intent(in) :: ndata,na
    double precision, intent(inout) :: a(:),alpha(:,:),beta(:),ochisq,alambda
    double precision, intent(out) :: chisq,covar(:,:)
    logical :: succmar
    ! performs one step of a chi-squared optimization of a non linear fit function 
    ! with a Gauss-Newton-Marquardt algorithm - the function fct must provide the 
    ! functionvalues y and the derivatives with respect to the fit parameters
    ! fct:        model function provided by the user
    !             x:       x value where the function should be evaluated
    !             a(:):    model parameters
    !             y:       function value
    !             dyda(:): derivatives of the function with respect to the model parameters
    ! x(:):       measured x data
    ! y(:):       measured y data
    ! sig(:):     uncertainty of the measured y data
    ! ndata:      number of data points
    ! na:         number of fit parameters
    ! a(:):       guessed fitparameters as input and improved fit parameters as output
    ! alpha(:,:): internal matrix, the matrix from the last iteration step is needed
    ! covar(:,:): covariance matrix calculated if alambda is exactly 0.0d0
    ! beta(:):    internal vector, the vector from the last iteration step is needed
    ! chisq:      current chi-squared
    ! ochisq:     chi-squared from the last iteration step
    ! alambda:    Marquardt parameter lambda
    !             alambda < 0.0d0: the preparatory step is done
    !             alambda > 0.0d0: normal gnm iteration is performed
    !             alambda = 0.0d0: the covariance matrix is calculated
    ! succmar:    true: iteration was succesful - false: iteration failed

    integer :: i1
    double precision :: det,khad
    double precision, allocatable :: atry_da(:),vector(:),normal(:,:)
    character(len=50) :: formatstring

    if(alambda.lt.0.0d0) then
       alambda=0.3d-3
       call mrqcof(fct,x,y,sig,ndata,na,a,alpha,beta,chisq)
       ochisq=chisq
       print '(a)','iteration       chi**2        parameters a'
       write(formatstring,'(a,i0,a)') '(i9,e13.4,',na,'e13.4)'
       print formatstring,0,chisq,a
    end if
    if(alambda==0.0d0) then
       covar=alpha
       call luinv(covar,na,det,khad,woff=.true.)
    else
       allocate(normal(na,na))
       normal=alpha
       do i1=1,na
          normal(i1,i1)=alpha(i1,i1)*(1.0d0+alambda)
       end do
       allocate(atry_da(na))
       call lusol(normal,na,beta,atry_da,woff=.true.)
       atry_da=a+atry_da
       allocate(vector(na))
       call mrqcof(fct,x,y,sig,ndata,na,atry_da,normal,vector,chisq)
       if(chisq.le.ochisq) then
          alambda=alambda*0.2d0
          ochisq=chisq
          succmar=.true.
          alpha=normal
          beta=vector
          a=atry_da
       else
          alambda=alambda*5.0d0
          succmar=.false.
       end if
       deallocate(atry_da)
       deallocate(vector)
       deallocate(normal)
    end if

  end subroutine mrqmin

  subroutine mrqiter(fct,x,y,sig,acc,tmax,ndata,na,a,da,covar,modvar,modvar_confid,error)
    interface
       subroutine fct(x,a,y,dyda)
         double precision, intent(in) :: x,a(:)
         double precision, intent(out) :: y,dyda(:)
       end subroutine fct
    end interface
    double precision, intent(in) :: x(:),y(:),sig(:),acc
    integer, intent(in) :: ndata,na,tmax
    double precision, intent(inout) :: a(:)
    double precision, intent(out) :: da(:),covar(:,:),modvar,modvar_confid(2)
    integer, intent(out) :: error
    ! fit a function fct with model parameters a to data (x,y+/-sig)
    ! fct:           model function provided by the user
    !                x:       x value where the function should be evaluated
    !                a(:):    model parameters
    !                y:       function value
    !                dyda(:): derivatives of the function with respect to the model parameters
    ! x(:):          x data
    ! y(:):          y data
    ! sig(:):        error of the y data
    ! acc:           desired accuracy of the model parameters
    ! tmax:          maximum number of iterations
    ! ndata:         number of data points
    ! na:            number of model parameters
    ! a(:):          model parameters
    ! da(:):         error of the model parameters
    ! covar(:,:):    normed covariance matrix
    ! modvar:        variance of the model function
    ! modvar_confid: confindence interval of the variance of the model function
    ! error:       0->convergence reached, 1->iteration failed

    integer :: t,i1,i2
    double precision :: chisq,ochisq,alambda,delmax
    double precision, allocatable :: alpha(:,:),beta(:),aold(:)
    logical :: succmar
    character(len=50) :: formatstring

    error=1
    alambda=-1.0d0
    allocate(alpha(na,na))
    allocate(beta(na))
    allocate(aold(na))
    aold=a
    do t=1,tmax
       call mrqmin(fct,x,y,sig,ndata,na,a,alpha,covar,beta,chisq,ochisq,alambda,succmar)
       if(succmar) then
          write(formatstring,'(a,i0,a)') '(i9,e13.4,',na,'e13.4)'
          print formatstring,t,chisq,a
          delmax=maxval(abs((a-aold)/a))
          if(delmax.lt.acc) then
             error=0
             alambda=0.0d0
             call mrqmin(fct,x,y,sig,ndata,na,a,alpha,covar,beta,chisq,ochisq,alambda,succmar)
             modvar=chisq/dble(ndata-na)
             modvar_confid(2)=sqrt(2.0d0/dble(ndata-na))
             modvar_confid(1)=-modvar_confid(2)
             modvar_confid=modvar_confid+1.0d0
             do i1=1,na
                da(i1)=sqrt(covar(i1,i1))
                covar(i1,i1)=1.0d0
                do i2=i1+1,na
                   covar(i1,i2)=covar(i1,i2)/da(i1)
                   covar(i2,i1)=covar(i1,i2)
                end do
             end do
             print '(/,a,/)','mrqiter converged'
             print '(a,f6.3,a,2(f5.3,a))','variance of the model = ',modvar,' (should be between '&
                  &,modvar_confid(1),' and ',modvar_confid(2),')'
             if((modvar.gt.modvar_confid(1)).and.(modvar.lt.modvar_confid(2))) then
                print '(a)','good model function'
             else
                print '(a)','adapt model function!'
             end if
             exit
          else
             aold=a
          end if
       else
          print '(i9,a)',t,'   iteration was not successful'
       end if
    end do
    if(error.ne.0) then
       print '(a)','ERROR in mrqiter: convergence not reached'
    end if

    deallocate(alpha)
    deallocate(beta)
    deallocate(aold)

  end subroutine mrqiter

  subroutine spline(n,x,y,a,b,c,d)
    integer, intent(in) :: n
    double precision, intent(in) :: x(:),y(:)
    double precision, intent(out) :: a(:),b(:),c(:),d(:)
    ! performs a spline interpolation  y=a+b(x-x0)+c(x-x0)**2+d(x-x0)**3 with natural splines
    ! n:    number of datapoints that are interpolated
    ! x(:): data points that are interpolated
    ! y(:): function values at the data points x(:)
    ! a(:): coefficient a
    ! b(:): coefficient b
    ! c(:): coefficient c
    ! d(:): coefficient d

    double precision, allocatable :: h(:),r(:),md(:),ld(:),ud(:),inhom(:)
    integer :: i1,nm1,nm2

    nm1=n-1
    nm2=n-2
    allocate(h(nm1))
    allocate(r(nm1))
    allocate(md(nm2))
    allocate(ld(nm2))
    allocate(ud(nm2))
    allocate(inhom(nm2))

    do i1=1,nm1
       if(x(i1).ge.x(i1+1)) then
          print '(a)','ERROR in spline: x values are not strictly monotonically increasing'
          return
       end if
       h(i1)=x(i1+1)-x(i1)
       r(i1)=y(i1+1)-y(i1)
    end do
    do i1=1,n-2
       md(i1)=2.0d0*(h(i1)+h(i1+1))
       ld(i1)=h(i1)
       ud(i1)=h(i1+1)
       inhom(i1)=3.0d0*(r(i1+1)/h(i1+1)-r(i1)/h(i1))
    end do
    call trid(md,ld,ud,inhom,nm2)
    c(1)=0.0d0
    c(2:nm1)=inhom
    do i1=1,nm2
       a(i1)=y(i1)
       b(i1)=r(i1)/h(i1)-h(i1)*(c(i1+1)+2.0d0*c(i1))/3.0d0
       d(i1)=(c(i1+1)-c(i1))/(3.0d0*h(i1))
    end do
    a(nm1)=y(nm1)
    b(nm1)=r(nm1)/h(nm1)-h(nm1)*2.0d0*c(nm1)/3.0d0
    d(nm1)=-c(nm1)/(3.0d0*h(nm1))

    deallocate(h)
    deallocate(r)
    deallocate(md)
    deallocate(ld)
    deallocate(ud)
    deallocate(inhom)

  end subroutine spline

  subroutine splval(n,x,a,b,c,d,xx,yy,dyy,ddyy)
    integer, intent(in) :: n
    double precision, intent(in) :: x(:),a(:),b(:),c(:),d(:),xx
    double precision, intent(out) :: yy
    double precision, optional, intent(out) :: dyy,ddyy
    ! calculates the function value of the spline interpolation y=a+b(x-x0)+c(x-x0)**2+d(x-x0)**3
    ! n:    number of datapoints that were interpolated
    ! x(:): data points that were interpolated
    ! a(:): coefficient a
    ! b(:): coefficient b
    ! c(:): coefficient c
    ! d(:): coefficient d
    ! xx:   interpolation point
    ! yy:   function value at the interpolation point xx
    ! dyy:  value of the first derivative of the function at the interpolation point xx
    ! ddyy: value of the second derivative of the function at the interpolation point xx

    integer, parameter :: div_expol=10
    integer :: i1,i1l,i1u
    double precision :: max_expol,z

    max_expol=(x(n)-x(1))/dble(n*div_expol)
    if((xx.lt.(x(1)-max_expol)).or.(xx.gt.(x(n)+max_expol))) then
       print '(a)','ERROR in splval: xx is not in the interpolation interval'
       return
    end if
    i1u=n+1
    i1l=0
    do while(i1l.lt.i1u-1)
       i1=(i1u+i1l)/2
       if(xx.ge.x(i1)) then
          i1l=i1
       else
          i1u=i1
       end if
    end do
    i1=i1l
    if(i1.le.0) i1=1
    if(i1.ge.n) i1=n-1
    z=xx-x(i1)
    yy=((d(i1)*z+c(i1))*z+b(i1))*z+a(i1)
    if(present(dyy)) dyy=(3.0d0*d(i1)*z+2.0d0*c(i1))*z+b(i1)
    if(present(ddyy)) ddyy=6.0d0*d(i1)*z+2.0d0*c(i1)

  end subroutine splval

  subroutine fft_pow2(a,m,inv)
    complex(kind=kind(0.0d0)), intent(inout) :: a(:)
    integer, intent(in) :: m,inv
    ! perform a fast fouriere transformation by the methode of Cooley and Tukey
    ! a(:): data to transform as input and transformed data as output
    !       a_new(1)=sum(a_old)
    ! m:    2**m is the number ot data points
    ! inv:  1->fouriere transform, else->inverse fouriere transform

    integer :: n,nd2,i1,i2,i3,i4,ie1,ie2,ipie
    complex(kind=kind(0.0d0)) :: u,w,t
    double precision :: ang

    n=2**m
    nd2=n/2
    i2=1
    do i1=1,n-1
       if(i1.lt.i2) then
          t=a(i2)
          a(i2)=a(i1)
          a(i1)=t
       end if
       i3=nd2
       do while(i3.lt.i2)
          i2=i2-i3
          i3=i3/2
       end do
       i2=i2+i3
    end do
    ie1=1
    do i4=1,m
       ie2=ie1
       ie1=ie1+ie1
       u=(1.0d0,0.0d0)
       ang=m_pi/dble(ie2)
       w=cmplx(cos(ang),-sin(ang),kind(0.0d0))
       if(inv.eq.1) w=conjg(w)
       do i2=1,ie2
          do i1=i2,n,ie1
             ipie=i1+ie2
             t=a(ipie)*u
             a(ipie)=a(i1)-t
             a(i1)=a(i1)+t
          end do
          u=u*w
       end do
    end do
    if(inv.ne.1) then
       a=a/dble(n)
    end if

  end subroutine fft_pow2

  subroutine ft_n(a,n,inv)
    complex(kind=kind(0.0d0)), intent(inout) :: a(:)
    integer, intent(in) :: n,inv
    ! perform a O(n**2) process fouriere transformation
    ! a(:): data to transform as input and transformed data as output
    !       a_new(1)=sum(a_old)
    ! n:    n is the number of data points, it does not have to be a power of 2
    ! inv:  1->fouriere transform, else->inverse fouriere transform

    integer :: i1,i2
    complex(kind=kind(0.0d0)), allocatable :: aintern(:),expfac(:)
    double precision :: ang

    allocate(aintern(n))
    aintern=0.0d0
    allocate(expfac(n))
    ang=2.0d0*m_pi/dble(n)
    expfac(1)=(1.0d0,0.0d0)
    expfac(2)=cmplx(cos(ang),-sin(ang),kind(0.0d0))
    if(inv.eq.1) expfac(2)=conjg(expfac(2))
    do i1=3,n
       expfac(i1)=expfac(2)*expfac(i1-1)
    end do
    do i1=1,n
       do i2=1,n
          aintern(i2)=aintern(i2)+a(i1)*expfac(i1)**(i2-1)
       end do
    end do
    a=aintern
    if(inv.ne.1) a=a/dble(n)

  end subroutine ft_n

  subroutine ft(a,n,inv,info)
    complex(kind=kind(0.0d0)), intent(inout) :: a(:)
    integer, intent(in) :: n,inv
    logical, optional, intent(in) :: info
    ! perform a fouriere transformation, if n=2**m fast fouriere transformation is done
    ! a(:): data to transform as input and transformed data as output
    !       a_new(1)=sum(a_old)
    ! n:    n is the number of data points, it does not have to be a power of 2
    ! inv:  1->fouriere transform, else->inverse fouriere transform
    ! info: .true.->information about the used type of fouriere transformation displayed

    integer :: m

    m=nint(log(dble(n))/log(2.0d0))
    if(2**m.eq.n) then
       if(present(info).and.info) print '(a)','n is a power of 2 -> use fft'
       call fft_pow2(a,m,inv)
    else
       if(present(info).and.info) print '(a)','n is not a power of 2 -> use slow ft'
       call ft_n(a,n,inv)
    end if

  end subroutine ft

  subroutine odeint_ro(derivs,jacobn,ystart,nvar,x1,x2,acc,h1,hmin,nsav,nmax&
       &,xfield,yfield,error,dxsav,savstart)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    interface
       subroutine jacobn(x,y,dfdx,dfdy)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dfdx(:),dfdy(:,:)
       end subroutine jacobn
    end interface
    double precision, intent(inout) :: ystart(:)
    integer, intent(in) :: nvar,nmax
    integer, intent(out) :: nsav,error
    double precision, intent(in) :: x1,x2,acc,h1,hmin
    double precision, intent(out) :: xfield(:),yfield(:,:)
    double precision, optional, intent(in) :: dxsav
    logical, optional, intent(in) :: savstart
    ! solve an initial value problem with the method of Rosenbrock
    ! especially for stiff differential equations
    ! derivs:      calculates the derivatives of the first order differential equation
    !              x:       x value
    !              y(:):    function values
    !              dydx(:): derivatives
    ! jacobn:      calculates the second derivatives of the first order differential equation
    !              x:         x value
    !              y(:):      function values
    !              dfdx(:):   derivatives of f = dydx with respect to x
    !              dfdy(:,:): derivatives of f = dydx with respect to y
    ! ystart(:):   initial values of the initial value problem as input and last values as output
    ! nvar:        number of first order differential equations
    ! x1:          starting point x1
    ! x2:          final point x2
    ! acc:         desired accuracy of one step
    ! h1:          starting step size
    ! hmin:        minimum step size
    ! nsav:        number of saved steps
    ! nmax:        maximum number of steps
    ! xfield(:):   calculated x values
    ! yfield(:,:): calculated y values, first dimension: different y variables, second dimension: different x values
    ! error:       0->convergence reached, 1->iteration failed
    ! dxsav:       mimimum distance between two x values that should be saved
    ! savstart:    .true.-> start values are saved in xfield and yfield, else they are not

    double precision, parameter :: ctiny=1.0d-30
    integer :: nstp,nstpmax
    double precision :: h,hdid,hnext,x,xsav,dxsav_intern
    double precision, allocatable :: y(:),yscal(:),dydx(:)

    error=1
    allocate(y(nvar))
    allocate(yscal(nvar))
    allocate(dydx(nvar))
    x=x1
    h=sign(h1,x2-x1)
    y=ystart
    if(present(savstart).and.savstart) then
       nsav=1
       nstpmax=nmax-1
       xfield(nsav)=x
       xsav=x
       yfield(:,nsav)=ystart
    else
       nsav=0
       nstpmax=nmax
    end if
    if(.not.present(dxsav)) then
       dxsav_intern=0.0d0
    else
       dxsav_intern=abs(dxsav)
    end if
    do nstp=1,nstpmax
       call derivs(x,y,dydx)
       yscal=abs(y)+abs(h*dydx)+ctiny
       if((x+h-x2)*(x+h-x1).gt.0.0d0) h=x2-x
       call stiff(derivs,jacobn,y,dydx,nvar,x,h,acc,yscal,hdid,hnext)
       if(abs(x-xsav).gt.dxsav_intern) then
          nsav=nsav+1
          xfield(nsav)=x
          xsav=x
          yfield(:,nsav)=y
       end if
       if((x-x2)*(x2-x1).ge.0.0d0) then
          ystart=y
          if(xsav.ne.x) then
             nsav=nsav+1
             xfield(nsav)=x
             xsav=x
             yfield(:,nsav)=ystart
          end if
          error=0
          exit
       end if
       if(abs(hnext).lt.hmin) then
          print '(a)','WARNING in odeint_ro: minimum step size underrun'
       end if
       h=hnext
    end do
    if(error.ne.0) then
       print '(a)','ERROR in odeint_ro: too many steps, final point not reached'
    end if
    deallocate(y)
    deallocate(yscal)
    deallocate(dydx)

  end subroutine odeint_ro

  subroutine stiff(derivs,jacobn,y,dydx,n,x,htry,acc,yscal,hdid,hnext)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    interface
       subroutine jacobn(x,y,dfdx,dfdy)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dfdx(:),dfdy(:,:)
       end subroutine jacobn
    end interface
    double precision, intent(inout) :: y(:),dydx(:),x
    integer, intent(in) :: n
    double precision, intent(in) :: htry,acc,yscal(:)
    double precision, intent(out) :: hdid,hnext
    ! one single step for the method of Rosenbrock
    ! derivs:  calculates the derivatives of the first order differential equation
    !          x:       x value
    !          y(:):    function values
    !          dydx(:): derivatives
    ! jacobn:  calculates the second derivatives of the first order differential equation
    !          x:         x value
    !          y(:):      function values
    !          dfdx(:):   derivatives of f = dydx with respect to x
    !          dfdy(:,:): derivatives of f = dydx with respect to y
    ! y(:):    initial values of the initial value problem as input and next values as output
    ! dydx(:): derivatives according to the differntial equation
    ! n:       number of first order differential equations
    ! x:       current point as input and next point as output
    ! htry:    suggested step width
    ! acc:     desired accuracy of this step
    ! yscal:   scaled y values needed for the error estimation
    ! hdid:    performed step size
    ! hnext:   suggested step size for the next step

    double precision, parameter :: safety=0.9d0,shrnk=0.4d0,pshrnk=-1.0d0/4.0d0&
         &,grow=4.0d0,pgrow=-1.0d0/5.0d0,errcongrow=(grow/safety)**(1.0d0/pgrow)&
         &,errconshrnk=(shrnk/safety)**(1.0d0/pshrnk),gam=1.0d0/2.0d0&
         &,a21=2.0d0&
         &,a31=48.0d0/25.0d0,a32=6.0d0/25.0d0&
         &,c21=-8.0d0&
         &,c31=372.0d0/25.0d0,c32=12.0d0/5.0d0&
         &,c41=-112.0d0/125.0d0,c42=-54.0d0/125.0d0,c43=-2.0d0/5.0d0&
         &,b1=19.0d0/9.0d0,b2=1.0d0/2.0d0,b3=25.0d0/108.0d0,b4=125.0d0/108.0d0&
         &,e1=17.0d0/54.0d0,e2=7.0d0/36.0d0,e3=0.0d0,e4=125.0d0/108.0d0&
         &,c1x=1.0d0/2.0d0,c2x=-3.0d0/2.0d0,c3x=121.0d0/50.0d0,c4x=29.0d0/250.0d0&
         &,a2x=1.0d0,a3x=3.0d0/5.0d0
    integer :: i1
    integer, allocatable :: indx(:)
    double precision :: det,khad,errmax,h,xsav
    double precision, allocatable :: a(:,:),dfdx(:),dfdy(:,:),dysav(:),err(:)&
         &,g1(:),g2(:),g3(:),g4(:),ysav(:)
    logical :: stepok

    allocate(indx(n))
    allocate(a(n,n))
    allocate(dfdx(n))
    allocate(dfdy(n,n))
    allocate(dysav(n))
    allocate(err(n))
    allocate(g1(n))
    allocate(g2(n))
    allocate(g3(n))
    allocate(g4(n))
    allocate(ysav(n))
    xsav=x
    ysav=y
    dysav=dydx
    call jacobn(xsav,ysav,dfdx,dfdy)
    h=htry
    stepok=.false.
    do while(.not.stepok)
       a=-dfdy
       do i1=1,n
          a(i1,i1)=1.0d0/(gam*h)+a(i1,i1)
       end do
       call ludcmp(a,n,indx,det,khad)
       g1=dysav+h*c1x*dfdx
       call lubksb(a,n,indx,g1)
       y=ysav+a21*g1
       x=xsav+a2x*h
       call derivs(x,y,dydx)
       g2=dydx+h*c2x*dfdx+c21*g1/h
       call lubksb(a,n,indx,g2)
       y=ysav+a31*g1+a32*g2
       x=xsav+a3x*h
       call derivs(x,y,dydx)
       g3=dydx+h*c3x*dfdx+(c31*g1+c32*g2)/h
       call lubksb(a,n,indx,g3)
       g4=dydx+h*c4x*dfdx+(c41*g1+c42*g2+c43*g3)/h
       call lubksb(a,n,indx,g4)
       y=ysav+b1*g1+b2*g2+b3*g3+b4*g4
       err=e1*g1+e2*g2+e3*g3+e4*g4
       x=xsav+h
       if(x.eq.xsav)print '(a)','ERROR in stiff: stepsize not significant'
       errmax=maxval(abs(err/yscal))
       errmax=errmax/acc
       if(errmax.le.1.0d0) then
          hdid=h
          if(errmax.gt.errcongrow) then
             hnext=safety*h*errmax**pgrow
          else
             hnext=grow*h
          endif
          stepok=.true.
       else
          if(errmax.lt.errconshrnk) then
             hnext=safety*h*errmax**pshrnk
          else
             hnext=shrnk*h
          endif
          h=hnext
       endif
    end do
    deallocate(indx)
    deallocate(a)
    deallocate(dfdx)
    deallocate(dfdy)
    deallocate(dysav)
    deallocate(err)
    deallocate(g1)
    deallocate(g2)
    deallocate(g3)
    deallocate(g4)
    deallocate(ysav)

  end subroutine stiff

  subroutine odeint_rk(derivs,ystart,nvar,x1,x2,acc,h1,hmin,nsav,nmax,xfield,yfield,error&
       &,dxsav,savstart,rkchoose)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    double precision, intent(inout) :: ystart(:)
    integer, intent(in) :: nvar,nmax
    integer, intent(out) :: nsav,error
    double precision, intent(in) :: x1,x2,acc,h1,hmin
    double precision, intent(out) :: xfield(:),yfield(:,:)
    double precision, optional, intent(in) :: dxsav
    logical, optional, intent(in) :: savstart
    integer, optional, intent(in) :: rkchoose
    ! solve an initial value problem with the Runge-Kutta method, different methods can be choosen
    ! derivs:      calculates the derivatives of the first order differential equation
    !              x:       x value
    !              y(:):    function values
    !              dydx(:): derivatives
    ! ystart(:):   initial values of the initial value problem as input and last values as output
    ! nvar:        number of first order differential equations
    ! x1:          starting point x1
    ! x2:          final point x2
    ! acc:         desired accuracy of one step
    ! h1:          starting step size
    ! hmin:        minimum step size
    ! nsav:        number of saved steps
    ! nmax:        maximum number of steps
    ! xfield(:):   calculated x values
    ! yfield(:,:): calculated y values, first dimension: function values, second dimension: different x values
    ! error:       0->convergence reached, 1->iteration failed
    ! dxsav:       mimimum distance between two x values that should be saved
    ! savstart:    .true.-> start values are saved in xfield and yfield, else they are not
    ! rkchoose:    choose the desired Runge-Kutta method
    !              1->Runge-Kutta-Fehlberg
    !              2->Cash-Karp
    !              3->Dormand-Prince
    !              4->classical Runge-Kutta 4th order
    !              else or not present->Runge-Kutta-Fehlberg

    double precision, parameter :: ctiny=1.0d-30
    integer :: nstp,nstpmax,rk_intern
    double precision :: h,hnext,x,xsav,dxsav_intern
    double precision, allocatable :: y(:),yscal(:),dydx(:)

    error=1
    if(.not.present(rkchoose)) then
       rk_intern=1
    else if(rkchoose.eq.2) then
       rk_intern=2
    else if(rkchoose.eq.3) then
       rk_intern=3
    else if(rkchoose.eq.4) then
       rk_intern=4
    else
       rk_intern=1
    end if
    allocate(y(nvar))
    allocate(yscal(nvar))
    allocate(dydx(nvar))
    x=x1
    h=sign(h1,x2-x1)
    y=ystart
    if(present(savstart).and.savstart) then
       nsav=1
       nstpmax=nmax-1
       xfield(nsav)=x
       xsav=x
       yfield(:,nsav)=ystart
    else
       nsav=0
       nstpmax=nmax
    end if
    if(.not.present(dxsav)) then
       dxsav_intern=0.0d0
    else
       dxsav_intern=abs(dxsav)
    end if
    do nstp=1,nstpmax
       call derivs(x,y,dydx)
       yscal=abs(y)+abs(h*dydx)+ctiny
       if((x+h-x2)*(x+h-x1).gt.0.0d0) h=x2-x
       call rkqc(derivs,y,dydx,nvar,x,h,acc,yscal,hnext,rk_intern)
       if(abs(x-xsav).gt.dxsav_intern) then
          nsav=nsav+1
          xfield(nsav)=x
          xsav=x
          yfield(:,nsav)=y
       end if
       if((x-x2)*(x2-x1).ge.0.0d0) then
          ystart=y
          if(xsav.ne.x) then
             nsav=nsav+1
             xfield(nsav)=x
             xsav=x
             yfield(:,nsav)=ystart
          end if
          error=0
          exit
       end if
       if(abs(hnext).lt.hmin) then
          print '(a,3e16.6)','WARNING in odeint_rk: minimum step size underrun',h,hnext,x-x2
       end if
       h=hnext
    end do
    if(error.ne.0) then
       print '(a)','ERROR in odeint_rk: too many steps, final point not reached'
    end if
    deallocate(y)
    deallocate(yscal)
    deallocate(dydx)

  end subroutine odeint_rk

  subroutine rkqc(derivs,y,dydx,nvar,x,htry,acc,yscal,hnext,rkchoose)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    double precision, intent(inout) :: y(:),dydx(:),x
    integer, intent(in) :: nvar
    double precision, intent(in) :: htry,acc,yscal(:)
    double precision, intent(out) :: hnext
    integer, intent(in) :: rkchoose
    ! one single step of a Runge-Kutta method, different methods can be choosen
    ! derivs:   calculates the derivatives of the first order differential equation
    !           x:       x value
    !           y(:):    function values
    !           dydx(:): derivatives
    ! y(:):     initial values of the initial value problem as input and next values as output
    ! dydx(:):  derivatives according to the differntial equation
    ! nvar:     number of first order differential equations
    ! x:        current point as input and next point as output
    ! htry:     suggested step width
    ! acc:      desired accuracy of this step
    ! yscal:    scaled y values needed for the error estimation
    ! hnext:    suggested step size for the next step
    ! rkchoose: which Runge-Kutta method is chosen
    !           1->Runge-Kutta-Fehlberg
    !           2->Cash-Karp
    !           3->Dormand-Prince
    !           4->classical Runge-Kutta 4th order

    double precision, parameter :: safety=0.9d0,shrnk=0.4d0,pshrnk=-1.0d0/4.0d0,grow=4.0d0,pgrow=-1.0d0/5.0d0
    double precision, parameter :: errcongrow=(grow/safety)**(1.0d0/pgrow),errconshrnk=(shrnk/safety)**(1.0d0/pshrnk)
    double precision :: h,xsav,errmax
    double precision, allocatable :: ysav(:),dysav(:),yerr(:)
    logical :: stepok

    xsav=x
    allocate(ysav(nvar))
    allocate(dysav(nvar))
    allocate(yerr(nvar))
    ysav=y
    dysav=dydx
    h=htry
    stepok=.false.
    do while(.not.stepok)
       select case(rkchoose)
       case(1)
          call rkf(derivs,y,dydx,nvar,x,h,yerr)
       case(2)
          call rkck(derivs,y,dydx,nvar,x,h,yerr)
       case(3)
          call rkdp(derivs,y,dydx,nvar,x,h,yerr)
       case(4)
          call rkc4err(derivs,y,dydx,nvar,x,h,yerr)
       end select
       errmax=maxval(abs(yerr/yscal))
       errmax=errmax/acc
       if(errmax.gt.1.0d0) then
          if(errmax.gt.errconshrnk) then
             hnext=shrnk*h
          else
             hnext=safety*h*errmax**pshrnk
          end if
          h=hnext
          x=xsav
          y=ysav
          dydx=dysav
       else
          stepok=.true.
          if(errmax.lt.errcongrow) then
             hnext=grow*h
          else
             hnext=safety*h*errmax**pgrow
          end if
       end if
    end do
    deallocate(ysav)
    deallocate(dysav)
    deallocate(yerr)

  end subroutine rkqc

  subroutine rkf(derivs,y,dydx,nvar,x,h,yerr)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    double precision, intent(inout) :: y(:),dydx(:),x
    integer, intent(in) :: nvar
    double precision, intent(in) :: h
    double precision, intent(out) :: yerr(:)
    ! one single step of a Runge-Kutta methode in the Fehlberg version of error calculation
    ! derivs:  calculates the derivatives of the first order differential equation
    !          x:       x value
    !          y(:):    function values
    !          dydx(:): derivatives
    ! y(:):    initial values of the initial value problem as input and next values as output
    ! dydx(:): derivatives according to the differntial equation
    ! nvar:    number of first order differential equations
    ! x:       current point as input and next point as output
    ! h:       step size
    ! yerr(:): calculated estimated error

    double precision :: xsav
    double precision, allocatable :: ak2(:),ak3(:),ak4(:),ak5(:),ak6(:),ytemp(:)
    double precision, parameter :: a2=0.25d0,a3=3.0d0/8.0d0,a4=12.0d0/13.0d0,a5=1.0d0,a6=0.5d0&
         &,b21=0.25d0&
         &,b31=3.0d0/32.0d0,b32=9.0d0/32.0d0&
         &,b41=1932.0d0/2197.0d0,b42=-7200.0d0/2197.0d0,b43=7296.0d0/2197.0d0&
         &,b51=439.0d0/216.0d0,b52=-8.0d0,b53=3680.0d0/513.0d0,b54=-845.0d0/4104.0d0&
         &,b61=-8.0d0/27.0d0,b62=2.0d0,b63=-3544.0d0/2565.0d0,b64=1859.0d0/4104.0d0,b65=-11.0d0/40.0d0&
         &,c1=16.0d0/135.0d0,c3=6656.0d0/12825.0d0,c4=28561.0d0/56430.0d0,c5=-9.0d0/50.0d0,c6=2.0d0/55.0d0&
         &,dc1=c1-25.0d0/216.0d0,dc3=c3-1408.0d0/2565.0d0,dc4=c4-2197.0d0/4104.0d0,dc5=c5+1.0d0/5.0d0,dc6=c6

    allocate(ak2(nvar))
    allocate(ak3(nvar))
    allocate(ak4(nvar))
    allocate(ak5(nvar))
    allocate(ak6(nvar))
    allocate(ytemp(nvar))
    xsav=x
    ytemp=y+b21*h*dydx
    x=xsav+a2*h
    call derivs(x,ytemp,ak2)
    ytemp=y+h*(b31*dydx+b32*ak2)
    x=xsav+a3*h
    call derivs(x,ytemp,ak3)
    ytemp=y+h*(b41*dydx+b42*ak2+b43*ak3)
    x=xsav+a4*h
    call derivs(x,ytemp,ak4)
    ytemp=y+h*(b51*dydx+b52*ak2+b53*ak3+b54*ak4)
    x=xsav+a5*h
    call derivs(x,ytemp,ak5)
    ytemp=y+h*(b61*dydx+b62*ak2+b63*ak3+b64*ak4+b65*ak5)
    x=xsav+a6*h
    call derivs(x,ytemp,ak6)
    y=y+h*(c1*dydx+c3*ak3+c4*ak4+c5*ak5+c6*ak6)
    yerr=h*(dc1*dydx+dc3*ak3+dc4*ak4+dc5*ak5+dc6*ak6)
    x=xsav+h
    deallocate(ak2)
    deallocate(ak3)
    deallocate(ak4)
    deallocate(ak5)
    deallocate(ak6)
    deallocate(ytemp)

  end subroutine rkf

  subroutine rkck(derivs,y,dydx,nvar,x,h,yerr)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    double precision, intent(inout) :: y(:),dydx(:),x
    integer, intent(in) :: nvar
    double precision, intent(in) :: h
    double precision, intent(out) :: yerr(:)
    ! one single step of a Runge-Kutta methode in the Cash-Karp version of error calculation
    ! derivs:  calculates the derivatives of the first order differential equation
    !          x:       x value
    !          y(:):    function values
    !          dydx(:): derivatives
    ! y(:):    initial values of the initial value problem as input and next values as output
    ! dydx(:): derivatives according to the differntial equation
    ! nvar:    number of first order differential equations
    ! x:       current point as input and next point as output
    ! h:       step size
    ! yerr(:): calculated estimated error

    double precision :: xsav
    double precision, allocatable :: ak2(:),ak3(:),ak4(:),ak5(:),ak6(:),ytemp(:)
    double precision, parameter :: a2=0.2d0,a3=0.3d0,a4=0.6d0,a5=1.0d0,a6=0.875d0&
         &,b21=0.2d0&
         &,b31=3.0d0/40.0d0,b32=9.0d0/40.0d0&
         &,b41=0.3d0,b42=-0.9d0,b43=1.2d0&
         &,b51=-11.0d0/54.0d0,b52=2.50d0,b53=-70.0d0/27.0d0,b54=35.0d0/27.0d0&
         &,b61=1631.0d0/55296.0d0,b62=175.0d0/512.0d0,b63=575.0d0/13824.0d0,b64=44275.0d0/110592.0d0,b65=253.0d0/4096.0d0&
         &,c1=37.0d0/378.0d0,c3=250.0d0/621.0d0,c4=125.0d0/594.0d0,c6=512.0d0/1771.0d0&
         &,dc1=c1-2825.0d0/27648.0d0,dc3=c3-18575.0d0/48384.0d0,dc4=c4-13525.0d0/55296.0d0,dc5=-277.0d0/14336.0d0,dc6=c6-0.25d0

    allocate(ak2(nvar))
    allocate(ak3(nvar))
    allocate(ak4(nvar))
    allocate(ak5(nvar))
    allocate(ak6(nvar))
    allocate(ytemp(nvar))
    xsav=x
    ytemp=y+b21*h*dydx
    x=xsav+a2*h
    call derivs(x,ytemp,ak2)
    ytemp=y+h*(b31*dydx+b32*ak2)
    x=xsav+a3*h
    call derivs(x,ytemp,ak3)
    ytemp=y+h*(b41*dydx+b42*ak2+b43*ak3)
    x=xsav+a4*h
    call derivs(x,ytemp,ak4)
    ytemp=y+h*(b51*dydx+b52*ak2+b53*ak3+b54*ak4)
    x=xsav+a5*h
    call derivs(x,ytemp,ak5)
    ytemp=y+h*(b61*dydx+b62*ak2+b63*ak3+b64*ak4+b65*ak5)
    x=xsav+a6*h
    call derivs(x,ytemp,ak6)
    y=y+h*(c1*dydx+c3*ak3+c4*ak4+c6*ak6)
    yerr=h*(dc1*dydx+dc3*ak3+dc4*ak4+dc5*ak5+dc6*ak6)
    x=xsav+h
    deallocate(ak2)
    deallocate(ak3)
    deallocate(ak4)
    deallocate(ak5)
    deallocate(ak6)
    deallocate(ytemp)

  end subroutine rkck

  subroutine rkdp(derivs,y,dydx,nvar,x,h,yerr)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    double precision, intent(inout) :: y(:),dydx(:),x
    integer, intent(in) :: nvar
    double precision, intent(in) :: h
    double precision, intent(out) :: yerr(:)
    ! one single step of a Runge-Kutta methode in the Dormand-Prince version of error calculation
    ! derivs:  calculates the derivatives of the first order differential equation
    !          x:       x value
    !          y(:):    function values
    !          dydx(:): derivatives
    ! y(:):    initial values of the initial value problem as input and next values as output
    ! dydx(:): derivatives according to the differntial equation
    ! nvar:    number of first order differential equations
    ! x:       current point as input and next point as output
    ! h:       step size
    ! yerr(:): calculated estimated error

    double precision :: xsav
    double precision, allocatable :: ak2(:),ak3(:),ak4(:),ak5(:),ak6(:),ak7(:),ytemp(:)
    double precision, parameter :: a2=0.2d0,a3=0.3d0,a4=0.8d0,a5=8.0d0/9.0d0,a6=1.0d0,a7=1.0d0&
         &,b21=0.2d0&
         &,b31=3.0d0/40.0d0,b32=9.0d0/40.0d0&
         &,b41=44.0d0/45.0d0,b42=-56.0d0/15.0d0,b43=32.0d0/9.0d0&
         &,b51=19372.0d0/6561.0d0,b52=-25360.0d0/2187.0d0,b53=64448.0d0/6561.0d0,b54=-212.0d0/729.0d0&
         &,b61=9017.0d0/3168.0d0,b62=-355.0d0/33.0d0,b63=46732.0d0/5247.0d0,b64=49.0d0/176.0d0,b65=-5103.0d0/18656.0d0&
         &,b71=35.0d0/384.0d0,b73=500.0d0/1113.0d0,b74=125.0d0/192.0d0,b75=-2187.0d0/6784.0d0,b76=11.0d0/84.0d0&
         &,c1=5179.0d0/57600.0d0,c3=7571.0d0/16695.0d0,c4=393.0d0/640.0d0,c5=-92097.0d0/339200.0d0,c6=187.0d0/2100.0d0&
         &,c7=1.0d0/40.0d0&
         &,dc1=c1-35.0d0/384.0d0,dc3=c3-500.0d0/1113.0d0,dc4=c4-125.0d0/192.0d0,dc5=c5+2187.0d0/6784.0d0&
         &,dc6=c6-11.0d0/84.0d0,dc7=c7

    allocate(ak2(nvar))
    allocate(ak3(nvar))
    allocate(ak4(nvar))
    allocate(ak5(nvar))
    allocate(ak6(nvar))
    allocate(ak7(nvar))
    allocate(ytemp(nvar))
    xsav=x
    ytemp=y+b21*h*dydx
    x=xsav+a2*h
    call derivs(x,ytemp,ak2)
    ytemp=y+h*(b31*dydx+b32*ak2)
    x=xsav+a3*h
    call derivs(x,ytemp,ak3)
    ytemp=y+h*(b41*dydx+b42*ak2+b43*ak3)
    x=xsav+a4*h
    call derivs(x,ytemp,ak4)
    ytemp=y+h*(b51*dydx+b52*ak2+b53*ak3+b54*ak4)
    x=xsav+a5*h
    call derivs(x,ytemp,ak5)
    ytemp=y+h*(b61*dydx+b62*ak2+b63*ak3+b64*ak4+b65*ak5)
    x=xsav+a6*h
    call derivs(x,ytemp,ak6)
    ytemp=y+h*(b71*dydx+b73*ak3+b74*ak4+b75*ak5+b76*ak6)
    x=xsav+a7*h
    call derivs(x,ytemp,ak7)
    y=y+h*(c1*dydx+c3*ak3+c4*ak4+c5*ak5+c6*ak6+c7*ak7)
    yerr=h*(dc1*dydx+dc3*ak3+dc4*ak4+dc5*ak5+dc6*ak6+dc7*ak7)
    x=xsav+h
    deallocate(ak2)
    deallocate(ak3)
    deallocate(ak4)
    deallocate(ak5)
    deallocate(ak6)
    deallocate(ak7)
    deallocate(ytemp)

  end subroutine rkdp

  subroutine rkc4err(derivs,y,dydx,nvar,x,h,yerr)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    double precision, intent(inout) :: y(:),dydx(:),x
    integer, intent(in) :: nvar
    double precision, intent(in) :: h
    double precision, intent(out) :: yerr(:)
    ! calssical Runge-Kutte method 4th order, error calculation by h/2 method
    ! derivs:  calculates the derivatives of the first order differential equation
    !          x:       x value
    !          y(:):    function values
    !          dydx(:): derivatives
    ! y(:):    initial values of the initial value problem as input and next values as output
    ! dydx(:): derivatives according to the differntial equation
    ! nvar:    number of first order differential equations
    ! x:       current point as input and next point as output
    ! h:       step size
    ! yerr(:): calculated estimated error

    double precision :: hhalf,xsav
    double precision, allocatable :: ysav(:),dysav(:)

    xsav=x
    allocate(ysav(nvar))
    allocate(dysav(nvar))
    ysav=y
    dysav=dydx
    hhalf=0.5d0*h
    call rkc4(derivs,y,dydx,nvar,x,hhalf)
    call derivs(x,y,dydx)
    call rkc4(derivs,y,dydx,nvar,x,hhalf)
    if(x.eq.xsav) then
       print '(a)','WARNING in rkc4err: half stepsize not significant'
       x=xsav+h
       if(x.eq.xsav) then
          print '(a)','ERROR in rkc4err: stepsize not significant'
       end if
    end if
    call rkc4(derivs,ysav,dysav,nvar,xsav,h)
    yerr=y-ysav
    deallocate(ysav)
    deallocate(dysav)

  end subroutine rkc4err

  subroutine rkc4(derivs,y,dydx,nvar,x,h)
    interface
       subroutine derivs(x,y,dydx)
         double precision, intent(in) :: x,y(:)
         double precision, intent(out) :: dydx(:)
       end subroutine derivs
    end interface
    double precision, intent(inout) :: y(:),dydx(:),x
    integer, intent(in) :: nvar
    double precision, intent(in) :: h
    ! classical Runge-Kutte method 4th order
    ! derivs:  calculates the derivatives of the first order differential equation
    !          x:       x value
    !          y(:):    function values
    !          dydx(:): derivatives
    ! y(:):    initial values of the initial value problem as input and next values as output
    ! dydx(:): derivatives according to the differntial equation
    ! nvar:    number of first order differential equations
    ! x:       current point as input and next point as output
    ! h:       step size

    double precision, allocatable :: ysav(:),yki(:)

    allocate(ysav(nvar))
    allocate(yki(nvar))
    ysav=y
    y=y+h/6.0d0*dydx
    yki=ysav+h/2.0d0*dydx
    x=x+h/2.0d0
    call derivs(x,yki,dydx)
    y=y+h/3.0d0*dydx
    yki=ysav+h/2.0d0*dydx
    call derivs(x,yki,dydx)
    y=y+h/3.0d0*dydx
    yki=ysav+h*dydx
    x=x+h/2.0d0
    call derivs(x,yki,dydx)
    y=y+h/6.0d0*dydx
    deallocate(ysav)
    deallocate(yki)

  end subroutine rkc4

  real function getnan()
    ! returns the indefinite value nan (not a number) for every data type
    ! getnan: returns nan

    real :: dum

    dum=0.0
    getnan=dum/dum

  end function getnan

  real function getinf()
    ! returns the indefinite value nan (not a number) for every data type
    ! getnan: returns nan

    real :: dum

    dum=0.0
    getinf=1.0/dum

  end function getinf

  integer function factorial(n)
    integer, intent(in) :: n
    ! calculate n! (n factorial)
    ! n:         integer number as input
    ! factorial: n! as output

    integer :: i1

    i1=1
    factorial=product((/(i1,i1=1,n)/))
    
  end function factorial

  integer function nCr(n,k)
    integer, intent(in) :: n,k
    ! calculate the binomial coefficient (n,k)=n!/(k!(n-k)!)
    ! n:   integer n as input
    ! k:   integer k as input
    ! nCr: binomial coefficient as output

    integer :: i1

    nCr=1
    if(k.gt.n-k) then
       do i1=1,n-k
          nCr=nCr*(i1+k)/i1
       end do
    else
       do i1=1,k
          nCr=nCr*(i1+n-k)/i1
       end do
    end if

  end function nCr

  subroutine randinit(seed)
    integer, intent(inout) :: seed
    ! initialize the random number generator rand() with an arbitrary seed derived out of the system time
    ! seed: arbitrary integer number out of which the seed is generated, odd number recommended
    !       as output the used seed is handed over

    integer :: seed1
    real :: seedout

    call system_clock(count=seed1)
    seed1=seed1+time()
    seed=seed1*2+seed
    if(seed.lt.0) seed=-seed
    seedout=rand(seed)

  end subroutine randinit

  subroutine randinit_const(seed)
    integer, intent(in) :: seed
    ! initialize the random number generator rand() with a certain seed
    ! seed: seed used for the initialization

    real :: seedout

    seedout=rand(seed)

  end subroutine randinit_const

  real function rand_0()
    ! create a random number from a uniform distribution in the interval (0,1]

    rand_0=1.0-rand()

  end function rand_0

  real function randn()
    ! creates a gaussian distributed random number with mean 0 and sigma 1 (normal distribution)
    ! code from http://www.netlib.org/random/ downloaded on 06.10.2014
    ! part of the original comment from the downloaded code:
    ! Some of the code is adapted from Dagpunar's book:
    !     Dagpunar, J. 'Principles of random variate generation'
    !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
    ! Version 1.13, 2 October 2000
    ! Changes from version 1.01
    ! 1. The random_order, random_Poisson & random_binomial routines have been
    !    replaced with more efficient routines.
    ! 2. A routine, seed_random_number, has been added to seed the uniform random
    !    number generator.   This requires input of the required number of seeds
    !    for the particular compiler from a specified I/O unit such as a keyboard.
    ! 3. Made compatible with Lahey's ELF90.
    ! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.
    ! 5. INTENT for array f corrected in random_mvnorm.
    !
    !     Author: Alan Miller
    !             CSIRO Division of Mathematical & Information Sciences
    !             Private Bag 10, Clayton South MDC
    !             Clayton 3169, Victoria, Australia
    !     Phone: (+61) 3 9545-8016      Fax: (+61) 3 9545-8080
    !     e-mail: amiller @ bigpond.net.au
    !
    ! Adapted from the following Fortran 77 code
    !      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
    !      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
    !      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
    !
    !  The function random_normal() returns a normally distributed pseudo-random
    !  number with zero mean and unit variance.
    !
    !  The algorithm uses the ratio of uniforms method of A.J. Kinderman
    !  and J.F. Monahan augmented with quadratic bounding curves.

    real, parameter :: s=0.449871,t=-0.386595,a=0.19600,b=0.25472,r1=0.27597,r2=0.27846
    real :: u,v,x,y,q
    ! Generate P = (u,v) uniform in rectangle enclosing acceptance region
    do
       u=rand()
       v=rand()
       v=1.7156*(v-0.5)
       ! Evaluate the quadratic form
       x=u-s
       y=abs(v)-t
       q=x**2+y*(a*y-b*x)
       ! Accept P if inside inner ellipse
       if(q.lt.r1) exit
       ! Reject P if outside outer ellipse
       if(q.gt.r2) cycle
       ! Reject P if outside acceptance region or except if inside acceptance region
       if(v**2.lt.-4.0*LOG(u)*u**2) exit
    end do
    ! Return ratio of P's coordinates as the normal deviate
    randn=v/u

  end function randn

  real function randfermihot()
    ! get a random number of the x > 0 part of 1/(1+exp(x)) by inversion method

    real, parameter :: log2=log(2.0)

    randfermihot=-log(exp(log2*(1.0-rand()))-1.0)

  end function randfermihot

  subroutine random_init(seed)
    integer, intent(in) :: seed
    ! initialize the random number generator random_number() with an arbitrary seed derived out of the system time
    ! seed: arbitrary integer number out of which the seed is generated, odd number recommended

    integer :: i1,n,clock
    integer, allocatable :: seed1(:)

    call random_seed(size=n)
    allocate(seed1(n))
    call system_clock(count=clock)
    seed1=seed+clock+37*(/(i1-1,i1=1,n)/)
    call random_seed(put=seed1)
    deallocate(seed1)

  end subroutine random_init

  subroutine random_init_const(seed)
    integer, intent(in) :: seed
    ! initialize the random number generator random_number() with a certain seed
    ! seed: seed used for the initialization

    integer :: i1,n
    integer, allocatable :: seed1(:)

    call random_seed(size=n)
    allocate(seed1(n))
    seed1=seed+37*(/(i1-1,i1=1,n)/)
    call random_seed(put=seed1)
    deallocate(seed1)

  end subroutine random_init_const

  real function random_0()
    ! create a random number from a uniform distribution in the interval (0,1]

    call random_number(random_0)
    random_0=1.0-random_0

  end function random_0

  real function random_1()
    ! create a random number from a uniform distribution in the interval [0,1)

    call random_number(random_1)

  end function random_1

  double precision function random_exp(rate)
    ! create a random number from an exponential distribution with specified rate parameter
    ! in the interval (0,\infty)
    ! rate: rate parameter p(random_exp)=1/rate*exp(-rate*random_exp)
    double precision :: rate

    real :: random_01

    do
       call random_number(random_01)
       if(random_01.ne.0.0) exit
    end do
    random_exp=-log(real(random_01,8))/rate

  end function random_exp

  real function random_n()
    ! creates a gaussian distributed random number with mean 0 and sigma 1 (normal distribution)
    ! code from http://www.netlib.org/random/ downloaded on 06.10.2014
    ! part of the original comment from the downloaded code:
    ! Some of the code is adapted from Dagpunar's book:
    !     Dagpunar, J. 'Principles of random variate generation'
    !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
    ! Version 1.13, 2 October 2000
    ! Changes from version 1.01
    ! 1. The random_order, random_Poisson & random_binomial routines have been
    !    replaced with more efficient routines.
    ! 2. A routine, seed_random_number, has been added to seed the uniform random
    !    number generator.   This requires input of the required number of seeds
    !    for the particular compiler from a specified I/O unit such as a keyboard.
    ! 3. Made compatible with Lahey's ELF90.
    ! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.
    ! 5. INTENT for array f corrected in random_mvnorm.
    !
    !     Author: Alan Miller
    !             CSIRO Division of Mathematical & Information Sciences
    !             Private Bag 10, Clayton South MDC
    !             Clayton 3169, Victoria, Australia
    !     Phone: (+61) 3 9545-8016      Fax: (+61) 3 9545-8080
    !     e-mail: amiller @ bigpond.net.au
    !
    ! Adapted from the following Fortran 77 code
    !      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
    !      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
    !      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
    !
    !  The function random_normal() returns a normally distributed pseudo-random
    !  number with zero mean and unit variance.
    !
    !  The algorithm uses the ratio of uniforms method of A.J. Kinderman
    !  and J.F. Monahan augmented with quadratic bounding curves.

    real, parameter :: s=0.449871,t=-0.386595,a=0.19600,b=0.25472,r1=0.27597,r2=0.27846
    real :: u,v,x,y,q
    ! Generate P = (u,v) uniform in rectangle enclosing acceptance region
    do
       call random_number(u)
       call random_number(v)
       v=1.7156*(v-0.5)
       ! Evaluate the quadratic form
       x=u-s
       y=abs(v)-t
       q=x**2+y*(a*y-b*x)
       ! Accept P if inside inner ellipse
       if(q.lt.r1) exit
       ! Reject P if outside outer ellipse
       if(q.gt.r2) cycle
       ! Reject P if outside acceptance region or except if inside acceptance region
       if(v**2.lt.-4.0*LOG(u)*u**2) exit
    end do
    ! Return ratio of P's coordinates as the normal deviate
    random_n=v/u

  end function random_n

  real function random_fermihot()
    ! get a random number of the x > 0 part of 1/(1+exp(x)) by inversion method

    real, parameter :: log2=log(2.0)

    random_fermihot=-log(exp(log2*(random_0()))-1.0)

  end function random_fermihot

  subroutine linfit(x,y,n,k,d)
    double precision, intent(in) :: x(:),y(:)
    integer, intent(in) :: n
    double precision, optional, intent(out) :: k,d
    ! perform a linear fit of the given data y=kx+d
    ! x(:): x data for the fit
    ! y(:): y data for the fit
    ! n:    number of data points
    ! k:    gradient of the linear fit
    ! d:    shift of the linear fit

    double precision :: S_x,S_y,S_xx,S_xy,S_yy,Delta
    S_x=sum(x)
    S_y=sum(y)
    S_xx=sum(x*x)
    S_xy=sum(x*y)
    S_yy=sum(y*y)
    Delta=dble(n)*S_xx-S_x*S_x

    if(present(k)) k=(dble(n)*S_xy-S_x*S_y)/Delta
    if(present(d)) d=(S_xx*S_y-S_xy*S_x)/Delta

  end subroutine linfit

  subroutine linfit_error(x,y,sig,n,a,da,covar,modvar,modvar_confid)
    double precision, intent(in) :: x(:),y(:),sig(:)
    integer, intent(in) :: n
    double precision, intent(out) :: a(2)
    double precision, optional, intent(out) :: da(2),covar(2,2),modvar,modvar_confid(2)
    ! perform a linear fit of the given data y=a(1)x+a(2)
    ! x(:):          x data
    ! y(:):          y data
    ! sig(:):        error of the y data
    ! n:             number of data points
    ! a(2):          model parameters
    ! da(2):         error of the model parameters
    ! covar(2,2):    normed covariance matrix
    ! modvar:        variance of the model function
    ! modvar_confid: confindence interval of the variance of the model function

    double precision :: ysig,ysig_2,x_m,xx_m,xy_m,y_m,C_xx,C_xy,chisq
    integer :: i1

    ysig_2=1/sum(1/sig**2)
    ysig=sqrt(ysig_2)
    x_m=sum(x/sig**2)*ysig_2
    xx_m=sum(x**2/sig**2)*ysig_2
    C_xx=xx_m-x_m**2
    y_m=sum(y/sig**2)*ysig_2
    xy_m=sum(x*y/sig**2)*ysig_2
    C_xy=xy_m-x_m*y_m

    a(1)=C_xy/C_xx
    a(2)=y_m-x_m*a(1)
    if(present(da)) then
       da(1)=ysig_2/C_xx
       da(2)=sqrt(ysig_2*xx_m/C_xx)
    end if
    if(present(covar)) then
       covar=1.0d0
       covar(1,2)=-ysig_2*x_m/C_xx/sqrt(da(1)*da(2))
       covar(2,1)=covar(1,2)
    end if
    if(present(modvar)) then
       chisq=0.0d0
       do i1=1,n
          chisq=chisq+((y(i1)-a(1)*x(i1)-a(2))/sig(i1))**2
       end do
       modvar=chisq/dble(n-2)
    end if
    if(present(modvar_confid)) then
       modvar_confid(2)=sqrt(2.0d0/dble(n-2))
       modvar_confid(1)=-modvar_confid(2)
       modvar_confid=modvar_confid+1.0d0
    end if

  end subroutine linfit_error

  subroutine norm_vector(vector,n,norm)
    double precision, intent(inout) :: vector(:)
    integer, intent(in) :: n
    double precision, optional, intent(out) :: norm
    ! norm the vector
    ! vector(:): vector to norm
    ! n:         size of the vector
    ! norm:      calculated norm

    integer :: i1
    double precision :: norm_intern

    norm_intern=0.0d0
    do i1=1,n
       norm_intern=norm_intern+vector(i1)*vector(i1)
    end do
    norm_intern=sqrt(norm_intern)
    vector=vector/norm_intern
    if(present(norm)) norm=norm_intern

  end subroutine norm_vector

  subroutine lanczos(sparse_mat,index_mat,n_mat,n_tot,n_eigval,eigval,acc,error&
       &,eigvec,skip_diag,max_steps,ghost_state_level,startvec)
    double precision, intent(in) :: sparse_mat(:),acc
    integer, intent(in) :: index_mat(:,:),n_mat,n_tot
    integer, intent(inout) :: n_eigval
    integer, intent(out) :: error
    double precision, intent(out) :: eigval(:)
    double precision, optional, intent(out) :: eigvec(:,:)
    integer, optional, intent(in) :: skip_diag,max_steps
    double precision, optional, intent(in) :: ghost_state_level,startvec(:)
    ! find maximum absolute value eigenvalues of a huge symmetric matrix by the method of Lanczos
    ! sparse_mat(:):     matrix to diagonalize in spares form, hermiticity considered
    ! index_mat(:,1):    row indices of the spares matrix
    ! index_mat(:,2):    column indices of the spares matrix
    ! n_mat:             number of rows/columns of the matrix
    ! n_tot:             total number of non vanishing elements in the sparse matrix
    ! n_eigval:          number of eigenvalues that should be converged and calculated as input
    !                    and number of calculated eigenvalues as output
    ! eigval(:):         calculated eigenvalues
    ! acc:               desired accuracy
    ! error:             -1: convergence not reached in max_steps Lanczos steps
    !                    1: possibly a ghost state appeared
    !                    2: possibly a ghost state appeared in the first step
    !                    0: convergence reached
    ! eigvec(:,:):       calculated eigenvectors, row: basis dimension, column: different energies
    ! skip_diag:         skipping Lanczos steps between two convergence checks, > 0
    ! max_steps:         maximum number of Lanczos steps
    ! ghost_state_level: minimum level for the coupling constant in the Lanczos matrix before subspace is diagonalized
    ! startvec(:):       start vector for Lanczos method

    integer :: i1,lanczos_step,skip_diag_intern,max_steps_intern,last_diag,seed
    double precision, allocatable :: x_n_1(:),x_n_2(:),mat_diag(:),mat_coupling(:)
    double precision, allocatable :: eigval_lanczos(:),eigvec_lanczos(:,:)
    double precision :: ghost_state_level_intern

    if(present(skip_diag)) then
       skip_diag_intern=skip_diag
    else
       skip_diag_intern=4
    end if
    if(present(max_steps)) then
       max_steps_intern=max_steps
    else
       max_steps_intern=20
    end if
    if(present(ghost_state_level)) then
       ghost_state_level_intern=ghost_state_level
    else
       ghost_state_level_intern=1.0d-9
    end if
    error=-1

    print '(/,a)','calculating Lanczos eigenvalues'
    allocate(x_n_1(n_mat))
    if(present(startvec)) then
       x_n_1=startvec
    else
       seed=34123451
       call randinit(seed)
       call get_lanczos_vec_x_0(x_n_1,n_mat)
    end if
    allocate(x_n_2(n_mat))
    allocate(mat_diag(max_steps_intern+1))
    allocate(mat_coupling(max_steps_intern+1))
    mat_coupling(1)=0.0d0
    call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
    if(mat_coupling(2).lt.ghost_state_level_intern) error=2
    lanczos_step=1
    last_diag=lanczos_step/skip_diag_intern
    allocate(eigval_lanczos(n_eigval))
    allocate(eigvec_lanczos(max_steps_intern+1,n_eigval))
    eigval_lanczos=0.0d0
    do while((lanczos_step.lt.max_steps_intern).and.(error.lt.0))
       lanczos_step=lanczos_step+1
       call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
            &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
       if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
          error=1
          exit
       end if
       if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
          call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
               &,eigval_lanczos,eigvec_lanczos,error,acc)
          if(error.ge.0) exit
          last_diag=lanczos_step/skip_diag_intern
       end if
       if(lanczos_step.ge.max_steps_intern) exit
       lanczos_step=lanczos_step+1
       call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
            &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
       if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
          error=1
          exit
       end if
       if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
          call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
               &,eigval_lanczos,eigvec_lanczos,error,acc)
          if(error.ge.0) exit
          last_diag=lanczos_step/skip_diag_intern
       end if
    end do
    if(error.eq.0) then
       print '(/,a)','lanczos converged'
    else if(error.eq.1) then
       print '(/,a)','WARNING in lanczos: possible ghost state appeared'
       call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step-1,n_eigval&
            &,eigval_lanczos,eigvec_lanczos,error,acc)
       error=1
       if(lanczos_step.lt.n_eigval) then
          eigval_lanczos(lanczos_step+1:n_eigval)=getnan()
          if(present(eigvec)) then
             eigvec(:,lanczos_step+1:n_eigval)=getnan()
          end if
          n_eigval=lanczos_step
       end if
       lanczos_step=lanczos_step-1
    else if(error.eq.2) then
       print '(/,a)','WARNING in lanczos: possible ghost state appeared in the first step'
       eigval_lanczos(1)=mat_diag(1)
       if(n_eigval.gt.1) then
          eigval_lanczos(2:n_eigval)=getnan()
          if(present(eigvec)) then
             eigvec(:,2:n_eigval)=getnan()
          end if
          n_eigval=1
       end if
       if(present(eigvec)) then
          eigvec(:,1)=x_n_1
       end if
    else
       print '(/,a,i0,a)','ERROR in lanczos: no convergence after ',max_steps_intern,' iterations'
    end if
    eigval=eigval_lanczos
    deallocate(eigval_lanczos)

    if(present(eigvec).and.((error.eq.0).or.(error.eq.1))) then
       print '(/,a)','calculating Lanczos eigenvectors'
       if(present(startvec)) then
          x_n_1=startvec
       else
          call randinit_const(seed)
          call get_lanczos_vec_x_0(x_n_1,n_mat)
       end if
       do i1=1,n_eigval
          eigvec(:,i1)=x_n_1*eigvec_lanczos(1,i1)
       end do
       call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
       do i1=1,n_eigval
          eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(2,i1)
       end do
       max_steps_intern=lanczos_step
       lanczos_step=1
       do while(lanczos_step.lt.max_steps_intern)
          lanczos_step=lanczos_step+1
          call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
               &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
          do i1=1,n_eigval
             eigvec(:,i1)=eigvec(:,i1)+x_n_1*eigvec_lanczos(lanczos_step+1,i1)
          end do
          if(lanczos_step.ge.max_steps_intern) exit
          lanczos_step=lanczos_step+1
          call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
               &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
          do i1=1,n_eigval
             eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(lanczos_step+1,i1)
          end do
       end do
    end if

    deallocate(eigvec_lanczos)
    deallocate(x_n_1)
    deallocate(x_n_2)
    deallocate(mat_diag)
    deallocate(mat_coupling)

  contains

    subroutine get_lanczos_vec_x_0(x_0,n)
      double precision, intent(out) :: x_0(:)
      integer, intent(in) :: n
      ! initialize starting vector for Lanczos method with random numbers
      ! x_0(:): starting vector
      ! n:      size of the vector

      integer :: i1

      do i1=1,n
         x_0(i1)=rand()
      end do
      call norm_vector(x_0,n)

    end subroutine get_lanczos_vec_x_0

    subroutine get_lanczos_vec_x_1(x_0,sparse_mat,index_mat,n_mat,n_tot,x_1,eigval_0,eigval_1,coupling_1)
      double precision, intent(in) :: x_0(:),sparse_mat(:)
      integer, intent(in) :: index_mat(:,:),n_mat,n_tot
      double precision, intent(out) :: x_1(:),eigval_0,eigval_1,coupling_1
      ! perform the first Lanczos step
      ! x_0(:):         starting vector
      ! sparse_mat(:):  matrix of the eigenvalue problem
      ! index_mat(:,1): row indices of the sparse matrix
      ! index_mat(:,2): column indices of the sparse matrix
      ! n_mat:          number of rows/columns of the matrix
      ! n_tot:          total number of non vanishing elements in the sparse matrix
      ! x_1(:):         calculated first Lanczos vector
      ! eigval_0:       first diagonal term in the Lanczos matrix
      ! eigval_0:       second diagonal term in the Lanczos matrix
      ! coupling_1:     first offdiagonal term in the Lanczos matrix

      integer :: i1

      call get_expval(sparse_mat,index_mat,n_tot,x_0,eigval_0)
      x_1=-eigval_0*x_0
      do i1=1,n_tot
         if(index_mat(i1,1).eq.index_mat(i1,2)) then
            x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,1))
         else
            x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,2))
            x_1(index_mat(i1,2))=x_1(index_mat(i1,2))+sparse_mat(i1)*x_0(index_mat(i1,1))
         end if
      end do
      call norm_vector(x_1,n_mat,coupling_1)
      call get_expval(sparse_mat,index_mat,n_tot,x_1,eigval_1)

    end subroutine get_lanczos_vec_x_1

    subroutine get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n,x_n1,eigval_n,coupling_n&
         &,eigval_n1,coupling_n1)
      double precision, intent(in) :: sparse_mat(:),eigval_n,coupling_n
      integer, intent(in) :: index_mat(:,:),n_mat,n_tot
      double precision, intent(inout) :: x_n(:),x_n1(:)
      double precision, intent(out) :: eigval_n1,coupling_n1
      ! perform a general Lanczos step
      ! sparse_mat(:):  matrix of the eigenvalue problem
      ! index_mat(:,1): row indices of the sparse matrix
      ! index_mat(:,2): column indices of the sparse matrix
      ! n_mat:          number of rows/columns of the matrix
      ! n_tot:          total number of non vanishing elements in the sparse matrix
      ! x_n(:):         vector for the n-th interation
      ! x_n1(:):        vector for the (n-1)-th iteration as input and (n+1)-th interation as output
      ! eigval_n:       n-th diagonal term in the Lanczos matrix
      ! coupling_n:     n-th offdiagonal term in the Lanczos matrix
      ! eigval_n1:      (n+1)-th diagonal term in the Lanczos matrix
      ! coupling_n1:    (n+1)-th offdiagonal term in the Lanczos matrix

      integer :: i1

      x_n1=-eigval_n*x_n-coupling_n*x_n1
      do i1=1,n_tot
         if(index_mat(i1,1).eq.index_mat(i1,2)) then
            x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,1))
         else
            x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,2))
            x_n1(index_mat(i1,2))=x_n1(index_mat(i1,2))+sparse_mat(i1)*x_n(index_mat(i1,1))
         end if
      end do
      call norm_vector(x_n1,n_mat,coupling_n1)
      call get_expval(sparse_mat,index_mat,n_tot,x_n1,eigval_n1)

    end subroutine get_lanczos_vec_x_n_plus_1

    subroutine get_expval(sparse_mat,index_mat,n_tot,vector,expval)
      double precision, intent(in) :: sparse_mat(:),vector(:)
      integer, intent(in) :: index_mat(:,:),n_tot
      double precision, intent(out) :: expval
      ! calculate the expectation value of a vector and the sparse matrix (hermiticity used)
      ! sparse_mat(:):  matrix to calculate the expectation value
      ! index_mat(:,1): row indices of the sparse matrix
      ! index_mat(:,2): column indices of the sparse matrix
      ! n_tot:          total number of non vanishing elements in the sparse matrix
      ! vector(:):      vector to calculate the expectation value
      ! expval:         expectation value
      
      integer :: i1

      expval=0.0d0
      do i1=1,n_tot
         if(index_mat(i1,1).eq.index_mat(i1,2)) then
            expval=expval+vector(index_mat(i1,1))*sparse_mat(i1)*vector(index_mat(i1,2))
         else
            expval=expval+2.0d0*vector(index_mat(i1,1))*sparse_mat(i1)*vector(index_mat(i1,2))
         end if
      end do

    end subroutine get_expval

    subroutine  diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step&
         &,n_eigval,eigval_lanczos,eigvec_lanczos,error,acc)
      double precision, intent(in) :: mat_diag(:),mat_coupling(:),acc
      integer, intent(in) :: lanczos_step,n_eigval
      double precision, intent(inout) :: eigval_lanczos(:)
      double precision, intent(out) :: eigvec_lanczos(:,:)
      integer, intent(out) :: error
      ! calculate the lowest eigenvalues in the spaned subspace and check if they converged
      ! mat_diag(:):         main diagonal of the Lanczos matrix
      ! mat_coupling(:):     first minor diagonal of the Lanczos matrix
      ! lanczos_step:        current number of executed Lanczos steps
      ! n_eigval:            number of eigenvalues that should be calculated
      ! eigval_lanczos(:):   calculated lowest eigenvalues of the Lanczos matrix
      ! eigvec_lanczos(:,:): eigenvectors of the calculated lowest eigenvalues of the Lanczos matrix
      ! error:               unchanged: no convergence for this check
      !                      0: convergence reached
      ! acc:                 desired accuracy

      double precision, allocatable :: lanczos_diag(:),lanczos_offdiag(:),eigvec(:,:)
      double precision :: error_eigval
      integer :: i1,error_diag

      allocate(eigvec(lanczos_step+1,lanczos_step+1))
      allocate(lanczos_diag(lanczos_step+1))
      allocate(lanczos_offdiag(lanczos_step+1))
      eigvec=0.0d0
      do i1=1,lanczos_step+1
         eigvec(i1,i1)=1.0d0
      end do
      lanczos_diag=mat_diag
      lanczos_offdiag=mat_coupling
      call tql2(lanczos_step+1,lanczos_diag,lanczos_offdiag,eigvec,error_diag)
      deallocate(lanczos_offdiag)
      error_eigval=maxval(abs((eigval_lanczos-lanczos_diag(1:n_eigval))/eigval_lanczos))
      eigval_lanczos=lanczos_diag(1:n_eigval)
      eigvec_lanczos=eigvec(:,1:n_eigval)
      if(error_eigval.lt.acc) error=0
      print '(a,i3,a,f22.15)','Iteration ',lanczos_step,', Lanczos lowest eigenvalue: ',eigval_lanczos(1)
      deallocate(lanczos_diag)
      deallocate(eigvec)

    end subroutine diagonalize_mat_and_check_finished

  end subroutine lanczos

  subroutine norm_vector_complex(vector,n,norm)
    complex(kind=kind(0.0d0)), intent(inout) :: vector(:)
    integer, intent(in) :: n
    double precision, optional, intent(out) :: norm
    ! norm the vector, complex version
    ! vector(:): vector to norm
    ! n:         size of the vector
    ! norm:      calculated norm

    integer :: i1
    double precision :: norm_intern

    norm_intern=0.0d0
    do i1=1,n
       norm_intern=norm_intern+abs(vector(i1))**2
    end do
    norm_intern=sqrt(norm_intern)
    vector=vector/norm_intern
    if(present(norm)) norm=norm_intern

  end subroutine norm_vector_complex

  subroutine lanczos_complex(sparse_mat,index_mat,n_mat,n_tot,n_eigval,eigval,acc,error&
       &,eigvec,skip_diag,max_steps,ghost_state_level,startvec)
    complex(kind=kind(0.0d0)), intent(in) :: sparse_mat(:)
    integer, intent(in) :: index_mat(:,:),n_mat,n_tot
    integer, intent(inout) :: n_eigval
    double precision, intent(in) :: acc
    integer, intent(out) :: error
    double precision, intent(out) :: eigval(:)
    complex(kind=kind(0.0d0)), optional, intent(out) :: eigvec(:,:)
    integer, optional, intent(in) :: skip_diag,max_steps
    double precision, optional, intent(in) :: ghost_state_level
    complex(kind=kind(0.0d0)), optional, intent(in) :: startvec(:)
    ! find maximum absolute value eigenvalues of a huge Hermitian symmetric matrix by the method of Lanczos, complex version
    ! sparse_mat(:):     matrix to diagonalize in spares form, hermiticity considered
    ! index_mat(:,1):    row indices of the spares matrix
    ! index_mat(:,2):    column indices of the spares matrix
    ! n_mat:             number of rows/columns of the matrix
    ! n_tot:             total number of non vanishing elements in the sparse matrix
    ! n_eigval:          number of eigenvalues that should be converged and calculated as input
    !                    and number of calculated eigenvalues as output
    ! eigval(:):         calculated eigenvalues
    ! acc:               desired accuracy
    ! error:             -1: convergence not reached in max_steps Lanczos steps
    !                    1: possibly a ghost state appeared
    !                    2: possibly a ghost state appeared in the first step
    !                    0: convergence reached
    ! eigvec(:,:):       calculated eigenvectors, row: basis dimension, column: different energies
    ! skip_diag:         skipping Lanczos steps between two convergence checks, > 0
    ! max_steps:         maximum number of Lanczos steps
    ! ghost_state_level: minimum level for the coupling constant in the Lanczos matrix before subspace is diagonalized
    ! startvec(:):       start vector for Lanczos method

    integer :: i1,lanczos_step,skip_diag_intern,max_steps_intern,last_diag,seed
    complex(kind=kind(0.0d0)), allocatable :: x_n_1(:),x_n_2(:)
    double precision, allocatable :: mat_diag(:),mat_coupling(:),eigval_lanczos(:),eigvec_lanczos(:,:)
    double precision :: ghost_state_level_intern

    if(present(skip_diag)) then
       skip_diag_intern=skip_diag
    else
       skip_diag_intern=4
    end if
    if(present(max_steps)) then
       max_steps_intern=max_steps
    else
       max_steps_intern=20
    end if
    if(present(ghost_state_level)) then
       ghost_state_level_intern=ghost_state_level
    else
       ghost_state_level_intern=1.0d-9
    end if
    error=-1

    print '(/,a)','calculating Lanczos eigenvalues'
    allocate(x_n_1(n_mat))
    if(present(startvec)) then
       x_n_1=startvec
    else
       seed=34123451
       call randinit(seed)
       call get_lanczos_vec_x_0(x_n_1,n_mat)
    end if
    allocate(x_n_2(n_mat))
    allocate(mat_diag(max_steps_intern+1))
    allocate(mat_coupling(max_steps_intern+1))
    mat_coupling(1)=0.0d0
    call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
    if(mat_coupling(2).lt.ghost_state_level_intern) error=2
    lanczos_step=1
    last_diag=lanczos_step/skip_diag_intern
    allocate(eigval_lanczos(n_eigval))
    allocate(eigvec_lanczos(max_steps_intern+1,n_eigval))
    eigval_lanczos=0.0d0
    do while((lanczos_step.lt.max_steps_intern).and.(error.lt.0))
       lanczos_step=lanczos_step+1
       call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
            &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
       if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
          error=1
          exit
       end if
       if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
          call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
               &,eigval_lanczos,eigvec_lanczos,error,acc)
          if(error.ge.0) exit
          last_diag=lanczos_step/skip_diag_intern
       end if
       if(lanczos_step.ge.max_steps_intern) exit
       lanczos_step=lanczos_step+1
       call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
            &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
       if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
          error=1
          exit
       end if
       if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
          call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
               &,eigval_lanczos,eigvec_lanczos,error,acc)
          if(error.ge.0) exit
          last_diag=lanczos_step/skip_diag_intern
       end if
    end do
    if(error.eq.0) then
       print '(/,a)','lanczos converged'
    else if(error.eq.1) then
       print '(/,a)','WARNING in lanczos: possible ghost state appeared'
       call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step-1,n_eigval&
            &,eigval_lanczos,eigvec_lanczos,error,acc)
       error=1
       if(lanczos_step.lt.n_eigval) then
          eigval_lanczos(lanczos_step+1:n_eigval)=getnan()
          if(present(eigvec)) then
             eigvec(:,lanczos_step+1:n_eigval)=getnan()
          end if
          n_eigval=lanczos_step
       end if
       lanczos_step=lanczos_step-1
    else if(error.eq.2) then
       print '(/,a)','WARNING in lanczos: possible ghost state appeared in the first step'
       eigval_lanczos(1)=mat_diag(1)
       if(n_eigval.gt.1) then
          eigval_lanczos(2:n_eigval)=getnan()
          if(present(eigvec)) then
             eigvec(:,2:n_eigval)=getnan()
          end if
          n_eigval=1
       end if
       if(present(eigvec)) then
          eigvec(:,1)=x_n_1
       end if
    else
       print '(/,a,i0,a)','ERROR in lanczos: no convergence after ',max_steps_intern,' iterations'
    end if
    eigval=eigval_lanczos
    deallocate(eigval_lanczos)

    if(present(eigvec).and.((error.eq.0).or.(error.eq.1))) then
       print '(/,a)','calculating Lanczos eigenvectors'
       if(present(startvec)) then
          x_n_1=startvec
       else
          call randinit_const(seed)
          call get_lanczos_vec_x_0(x_n_1,n_mat)
       end if
       do i1=1,n_eigval
          eigvec(:,i1)=x_n_1*eigvec_lanczos(1,i1)
       end do
       call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
       do i1=1,n_eigval
          eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(2,i1)
       end do
       max_steps_intern=lanczos_step
       lanczos_step=1
       do while(lanczos_step.lt.max_steps_intern)
          lanczos_step=lanczos_step+1
          call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
               &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
          do i1=1,n_eigval
             eigvec(:,i1)=eigvec(:,i1)+x_n_1*eigvec_lanczos(lanczos_step+1,i1)
          end do
          if(lanczos_step.ge.max_steps_intern) exit
          lanczos_step=lanczos_step+1
          call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
               &,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
          do i1=1,n_eigval
             eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(lanczos_step+1,i1)
          end do
       end do
    end if

    deallocate(eigvec_lanczos)
    deallocate(x_n_1)
    deallocate(x_n_2)
    deallocate(mat_diag)
    deallocate(mat_coupling)

  contains

    subroutine get_lanczos_vec_x_0(x_0,n)
      complex(kind=kind(0.0d0)), intent(out) :: x_0(:)
      integer, intent(in) :: n
      ! initialize starting vector for Lanczos method with random numbers
      ! x_0(:): starting vector
      ! n:      size of the vector

      integer :: i1

      do i1=1,n
         x_0(i1)=cmplx(rand(),rand())
      end do
      call norm_vector_complex(x_0,n)

    end subroutine get_lanczos_vec_x_0

    subroutine get_lanczos_vec_x_1(x_0,sparse_mat,index_mat,n_mat,n_tot,x_1,eigval_0,eigval_1,coupling_1)
      complex(kind=kind(0.0d0)), intent(in) :: x_0(:),sparse_mat(:)
      integer, intent(in) :: index_mat(:,:),n_mat,n_tot
      complex(kind=kind(0.0d0)), intent(out) :: x_1(:)
      double precision, intent(out) :: eigval_0,eigval_1,coupling_1
      ! perform the first Lanczos step
      ! x_0(:):         starting vector
      ! sparse_mat(:):  matrix of the eigenvalue problem
      ! index_mat(:,1): row indices of the sparse matrix
      ! index_mat(:,2): column indices of the sparse matrix
      ! n_mat:          number of rows/columns of the matrix
      ! n_tot:          total number of non vanishing elements in the sparse matrix
      ! x_1(:):         calculated first Lanczos vector
      ! eigval_0:       first diagonal term in the Lanczos matrix
      ! eigval_0:       second diagonal term in the Lanczos matrix
      ! coupling_1:     first offdiagonal term in the Lanczos matrix

      integer :: i1

      call get_expval(sparse_mat,index_mat,n_tot,x_0,eigval_0)
      x_1=-eigval_0*x_0
      do i1=1,n_tot
         if(index_mat(i1,1).eq.index_mat(i1,2)) then
            x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,1))
         else
            x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,2))
            x_1(index_mat(i1,2))=x_1(index_mat(i1,2))+conjg(sparse_mat(i1))*x_0(index_mat(i1,1))
         end if
      end do
      call norm_vector_complex(x_1,n_mat,coupling_1)
      call get_expval(sparse_mat,index_mat,n_tot,x_1,eigval_1)

    end subroutine get_lanczos_vec_x_1

    subroutine get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n,x_n1,eigval_n,coupling_n&
         &,eigval_n1,coupling_n1)
      complex(kind=kind(0.0d0)), intent(in) :: sparse_mat(:)
      double precision, intent(in) :: eigval_n,coupling_n
      integer, intent(in) :: index_mat(:,:),n_mat,n_tot
      complex(kind=kind(0.0d0)), intent(inout) :: x_n(:),x_n1(:)
      double precision, intent(out) :: eigval_n1,coupling_n1
      ! perform a general Lanczos step
      ! sparse_mat(:):  matrix of the eigenvalue problem
      ! index_mat(:,1): row indices of the sparse matrix
      ! index_mat(:,2): column indices of the sparse matrix
      ! n_mat:          number of rows/columns of the matrix
      ! n_tot:          total number of non vanishing elements in the sparse matrix
      ! x_n(:):         vector for the n-th interation
      ! x_n1(:):        vector for the (n-1)-th iteration as input and (n+1)-th interation as output
      ! eigval_n:       n-th diagonal term in the Lanczos matrix
      ! coupling_n:     n-th offdiagonal term in the Lanczos matrix
      ! eigval_n1:      (n+1)-th diagonal term in the Lanczos matrix
      ! coupling_n1:    (n+1)-th offdiagonal term in the Lanczos matrix

      integer :: i1

      x_n1=-eigval_n*x_n-coupling_n*x_n1
      do i1=1,n_tot
         if(index_mat(i1,1).eq.index_mat(i1,2)) then
            x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,1))
         else
            x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,2))
            x_n1(index_mat(i1,2))=x_n1(index_mat(i1,2))+conjg(sparse_mat(i1))*x_n(index_mat(i1,1))
         end if
      end do
      call norm_vector_complex(x_n1,n_mat,coupling_n1)
      call get_expval(sparse_mat,index_mat,n_tot,x_n1,eigval_n1)

    end subroutine get_lanczos_vec_x_n_plus_1

    subroutine get_expval(sparse_mat,index_mat,n_tot,vector,expval)
      complex(kind=kind(0.0d0)), intent(in) :: sparse_mat(:),vector(:)
      integer, intent(in) :: index_mat(:,:),n_tot
      double precision, intent(out) :: expval
      ! calculate the expectation value of a vector and the sparse matrix (hermiticity used)
      ! sparse_mat(:):  matrix to calculate the expectation value
      ! index_mat(:,1): row indices of the sparse matrix
      ! index_mat(:,2): column indices of the sparse matrix
      ! n_tot:          total number of non vanishing elements in the sparse matrix
      ! vector(:):      vector to calculate the expectation value
      ! expval:         expectation value
      
      integer :: i1

      expval=0.0d0
      do i1=1,n_tot
         if(index_mat(i1,1).eq.index_mat(i1,2)) then
            expval=expval+real(sparse_mat(i1))*abs(vector(index_mat(i1,1)))**2
         else
            expval=expval+2.0d0*real(conjg(vector(index_mat(i1,1)))*sparse_mat(i1)*vector(index_mat(i1,2)))
         end if
      end do

    end subroutine get_expval

    subroutine  diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step&
         &,n_eigval,eigval_lanczos,eigvec_lanczos,error,acc)
      double precision, intent(in) :: mat_diag(:),mat_coupling(:),acc
      integer, intent(in) :: lanczos_step,n_eigval
      double precision, intent(inout) :: eigval_lanczos(:)
      double precision, intent(out) :: eigvec_lanczos(:,:)
      integer, intent(out) :: error
      ! calculate the lowest eigenvalues in the spaned subspace and check if they converged
      ! mat_diag(:):         main diagonal of the Lanczos matrix
      ! mat_coupling(:):     first minor diagonal of the Lanczos matrix
      ! lanczos_step:        current number of executed Lanczos steps
      ! n_eigval:            number of eigenvalues that should be calculated
      ! eigval_lanczos(:):   calculated lowest eigenvalues of the Lanczos matrix
      ! eigvec_lanczos(:,:): eigenvectors of the calculated lowest eigenvalues of the Lanczos matrix
      ! error:               unchanged: no convergence for this check
      !                      0: convergence reached
      ! acc:                 desired accuracy

      double precision, allocatable :: lanczos_diag(:),lanczos_offdiag(:),eigvec(:,:)
      double precision :: error_eigval
      integer :: i1,error_diag

      allocate(eigvec(lanczos_step+1,lanczos_step+1))
      allocate(lanczos_diag(lanczos_step+1))
      allocate(lanczos_offdiag(lanczos_step+1))
      eigvec=0.0d0
      do i1=1,lanczos_step+1
         eigvec(i1,i1)=1.0d0
      end do
      lanczos_diag=mat_diag
      lanczos_offdiag=mat_coupling
      call tql2(lanczos_step+1,lanczos_diag,lanczos_offdiag,eigvec,error_diag)
      deallocate(lanczos_offdiag)
      error_eigval=maxval(abs((eigval_lanczos-lanczos_diag(1:n_eigval))/eigval_lanczos))
      eigval_lanczos=lanczos_diag(1:n_eigval)
      eigvec_lanczos=eigvec(:,1:n_eigval)
      if(error_eigval.lt.acc) error=0
      print '(a,i3,a,f22.15)','Iteration ',lanczos_step,', Lanczos lowest eigenvalue: ',eigval_lanczos(1)
      deallocate(lanczos_diag)
      deallocate(eigvec)

    end subroutine diagonalize_mat_and_check_finished

  end subroutine lanczos_complex

  subroutine gram_schmidt_onb(vectors,n_vec,n_basis)
    double precision, intent(inout) :: vectors(:,:)
    integer, intent(in) :: n_vec,n_basis
    ! build an orthonormal basis (ONB) out of a set of linear independent vectors
    ! vectors(:,:): linear independent vectors as input and ONB as output
    !               one vector is stored as a column vector v_n=vectors(:,n)
    ! n_vec:        number of vectors that are handed over
    ! n_basis:      number of components of one vector

    double precision :: overlap_v
    integer :: i1,i2,i3

    call norm_vector(vectors(:,1),n_basis)
    do i1=2,n_vec
       do i2=1,i1-1
          overlap_v=0.0d0
          do i3=1,n_basis
             overlap_v=overlap_v+vectors(i3,i1)*vectors(i3,i2)
          end do
          vectors(:,i1)=vectors(:,i1)-overlap_v*vectors(:,i2)
       end do
       call norm_vector(vectors(:,i1),n_basis)
    end do

  end subroutine gram_schmidt_onb


  subroutine gram_schmidt_on_single_vector(onb,vector,n_vec,n_basis)
    double precision, intent(inout) :: onb(:,:),vector(:)
    integer, intent(in) :: n_vec,n_basis
    ! add one vector to an orthonormal basis (ONB)
    ! onb(:,:):  orthonormal basis with the vectors in the columns
    ! vector(:): vector to add to the ONB
    ! n_vec:     number of vectors that are handed over in the ONB
    ! n_basis:   number of components of one vector

    double precision :: overlap_v
    integer :: i1,i2

    do i1=1,n_vec
       overlap_v=0.0d0
       do i2=1,n_basis
          overlap_v=overlap_v+onb(i2,i1)*vector(i2)
       end do
       vector=vector-overlap_v*onb(:,i1)
    end do
    call norm_vector(vector,n_basis)

  end subroutine gram_schmidt_on_single_vector

end module numeric_lib