Added the entire code.

bin/Makefile

@@ -0,0 +1,32 @@
+#$Id: Makefile $
+
+.SUFFIXES: .f90 .o
+
+# Compilers
+
+# Laptop:
+FC = gfortran
+FFLAGS = -fstack-check -g -O2
+LIBS = -L/usr/local/lib -llapack
+#LIBS = -L/usr/local/lib -llapack -lrefblas
+
+MAIN1 = beta_and_derivatives.o
+OBJS1 = secant.o quadpack_double.o
+MODS1 = module.o
+EXEC1 = beta_and_derivatives
+
+MAIN2 = schmidt_decomposition.o
+OBJS2 = diasym.o SVD.o inverse_matrice.o
+EXEC2 = schmidt_decomposition
+
+##################### End of Configurable Options ###############
+all: beta_and_derivatives schmidt_decomposition
+
+beta_and_derivatives: ${MODS1} ${OBJS1} ${MAIN1} Makefile 
+	${FC} ${FFLAGS} -o ${EXEC1} ${MAIN1} ${OBJS1} ${MODS1} ${LIBS}
+
+schmidt_decomposition: ${OBJS2} ${MAIN2} Makefile 
+	${FC} ${FFLAGS} -o ${EXEC2} ${MAIN2} ${OBJS2} ${LIBS}
+
+.f90.o: ${MODS1} Makefile
+	${FC} ${FFLAGS} -c $<

bin/SVD.f90

@@ -0,0 +1,38 @@
+      SUBROUTINE SVD(A,U,S,S_matrix,V,M,N)
+      DOUBLE PRECISION:: A(M,N),U(M,M),VT(N,N),S(N),V(N,N),S_matrix(M,M)
+!
+! Program computes the matrix singular value decomposition. 
+! Using Lapack library.
+!
+! Programmed by sukhbinder Singh
+! 14th January 2011
+!      
+
+
+      DOUBLE PRECISION,ALLOCATABLE :: WORK(:)
+      INTEGER LDA,M,N,LWORK,LDVT,INFO
+      CHARACTER  JOBU, JOBVT
+
+      JOBU='A'
+      JOBVT='A'
+      LDA=M
+      LDU=M
+      LDVT=N
+
+      LWORK=MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N))
+
+      ALLOCATE(work(lwork))
+
+      CALL DGESVD(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,WORK, &
+                  LWORK, INFO )
+
+       DO I=1,LDVT
+         DO J=1,LDVT
+           V(I,J)=VT(J,I)
+         END DO
+       END DO
+       S_matrix=0
+       DO I=1,LDU
+         S_matrix(I,I)=S(I)
+       ENDDO
+      END SUBROUTINE SVD

bin/beta_and_derivatives.f90

@@ -0,0 +1,57 @@
+! This program is free software: you can redistribute it and/or modify
+!    it under the terms of the GNU General Public License as published by
+!    the Free Software Foundation, either version 3 of the License, or
+!    (at your option) any later version.
+!
+!    This program is distributed in the hope that it will be useful,
+!    but WITHOUT ANY WARRANTY; without even the implied warranty of
+!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+!    GNU General Public License for more details.
+!
+!    You should have received a copy of the GNU General Public License
+!    along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+program beta_and_derivatives
+
+!---------------!
+  !
+  ! ==== Global Data ====
+  !
+  USE beta_module
+  !
+  ! ==== Local Data ====
+  !
+
+real(dp) :: a
+real(dp) :: b
+real(dp) :: U,t
+real(dp) :: x,beta,dbetadt,dbetadU
+real(dp),parameter :: tol=1.0e-9_dp
+integer,parameter :: maxiter=300
+integer :: iter
+integer :: ier
+
+read(*,*) U,t
+
+a = 1.0_dp
+b = 2.0_dp
+
+U_tmp=U
+t_tmp=t
+call secant(f,a,b,x,tol,maxiter,iter,ier)
+beta = x
+
+a = 1.0_dp
+b = 2.0_dp
+U_tmp=U+0.000001_dp
+t_tmp=t
+call secant(f,a,b,x,tol,maxiter,iter,ier)
+dbetadU=(x-beta)/0.000001_dp
+U_tmp=U
+
+open (97, file='beta_dbetadU.dat',access='sequential')
+ write(97,15) beta,dbetadU
+close(97)
+
+    15 format(f15.10,f15.10,f15.10)
+end program

bin/diasym.f90

@@ -0,0 +1,19 @@
+!---------------------------------------------------------!
+!Calls the LAPACK diagonalization subroutine DSYEV        !
+!input:  a(n,n) = real symmetric matrix to be diagonalized!
+!            n  = size of a                               !
+!output: a(n,n) = orthonormal eigenvectors of a           !
+!        eig(n) = eigenvalues of a in ascending order     !
+!---------------------------------------------------------!
+!--------------------------!
+ subroutine diasym(a,eig,n)
+ implicit none
+
+ integer n,l,inf
+ real*8  a(n,n),eig(n),work(n*(3+n/2))
+
+ l=n*(3+n/2)
+ call dsyev('V','U',n,a,n,eig,work,l,inf)
+
+ end subroutine diasym
+!---------------------!

bin/inverse_matrice.f90

@@ -0,0 +1,79 @@
+subroutine inverse(a,c,n)
+!============================================================
+! Inverse matrix
+! Method: Based on Doolittle LU factorization for Ax=b
+! Alex G. December 2009
+!-----------------------------------------------------------
+! input ...
+! a(n,n) - array of coefficients for matrix A
+! n      - dimension
+! output ...
+! c(n,n) - inverse matrix of A
+! comments ...
+! the original matrix a(n,n) will be destroyed 
+! during the calculation
+!===========================================================
+implicit none 
+integer n
+double precision a(n,n), c(n,n)
+double precision L(n,n), U(n,n), b(n), d(n), x(n)
+double precision coeff
+integer i, j, k
+
+! step 0: initialization for matrices L and U and b
+! Fortran 90/95 aloows such operations on matrices
+L=0.0
+U=0.0
+b=0.0
+
+! step 1: forward elimination
+do k=1, n-1
+   do i=k+1,n
+      coeff=a(i,k)/a(k,k)
+      L(i,k) = coeff
+      do j=k+1,n
+         a(i,j) = a(i,j)-coeff*a(k,j)
+      end do
+   end do
+end do
+
+! Step 2: prepare L and U matrices 
+! L matrix is a matrix of the elimination coefficient
+! + the diagonal elements are 1.0
+do i=1,n
+  L(i,i) = 1.0
+end do
+! U matrix is the upper triangular part of A
+do j=1,n
+  do i=1,j
+    U(i,j) = a(i,j)
+  end do
+end do
+
+! Step 3: compute columns of the inverse matrix C
+do k=1,n
+  b(k)=1.0
+  d(1) = b(1)
+! Step 3a: Solve Ld=b using the forward substitution
+  do i=2,n
+    d(i)=b(i)
+    do j=1,i-1
+      d(i) = d(i) - L(i,j)*d(j)
+    end do
+  end do
+! Step 3b: Solve Ux=d using the back substitution
+  x(n)=d(n)/U(n,n)
+  do i = n-1,1,-1
+    x(i) = d(i)
+    do j=n,i+1,-1
+      x(i)=x(i)-U(i,j)*x(j)
+    end do
+    x(i) = x(i)/u(i,i)
+  end do
+! Step 3c: fill the solutions x(n) into column k of C
+  do i=1,n
+    c(i,k) = x(i)
+  end do
+  b(k)=0.0
+end do
+end subroutine inverse

bin/module.f90

@@ -0,0 +1,50 @@
+module beta_module
+
+integer,parameter :: dp=kind(1.d0)
+real, parameter :: Pi=3.14159265358979323846264338327950288419716939937510
+real(dp) :: t_tmp,U_tmp
+
+contains
+function integrand(x) result(res)
+
+real(dp),intent(in) :: x
+real(dp) :: res
+real(dp) :: num,den
+
+num = bessel_j0(x)*bessel_j1(x)
+den = x*(1.0_dp + exp(x*U_tmp/2.0_dp/t_tmp))
+
+res = num/den
+
+end function
+
+function f(x) result(res)
+
+real(dp),intent(in) :: x
+real(dp) :: res
+real(dp) :: a,b
+
+integer,parameter :: limit=1000
+integer,parameter :: lenw=limit*4
+real(dp) :: abserr
+real(dp),parameter :: epsabs=0.0_dp
+real(dp),parameter :: epsrel=0.00001_dp
+integer :: ier
+integer :: iwork(limit)
+integer,parameter :: inf=1
+integer :: last
+integer :: neval
+real(dp) :: work(lenw)
+
+! double, quadrature de Gauss pour integrer de 0 a inf. (inf = 1 
+! equivaut a infini
+call dqagi(integrand,0.0_dp,inf,epsabs,epsrel,a,abserr,neval, &
+           ier,limit,lenw,last,iwork,work)
+
+b = -sin(Pi/x)*2.0_dp*t_tmp*x/Pi
+
+res = b + 4.0_dp*t_tmp*a
+
+end function
+
+end module beta_module