Commit 84f14c42 authored by Luca's avatar Luca

refactored, cleaned up

parent ac315e3e
MODULE gpl_module
use mpl_module
implicit none
CONTAINS
FUNCTION GPL()
! integer :: m(:)
! complex(kind=prec) :: y(:)
complex(kind=prec) GPL
print*, 'hello, GPL'
END FUNCTION GPL
END MODULE gpl_module
......@@ -7,7 +7,7 @@ FFLAGS=-fdefault-real-8 -cpp
LD=gfortran
OBJ=functions.o
OBJ=mpl_module.o gpl_module.o
# Rules for main fortran files
......
MODULE functions
MODULE mpl_module
implicit none
integer, parameter :: prec = selected_real_kind(15,32)
......@@ -29,18 +29,18 @@ CONTAINS
polylog = sum(x**j / j**m)
END FUNCTION polylog
FUNCTION multiple_polylog_converges(x)
FUNCTION MPL_converges(x)
! checks if the MPL converges
complex(kind=prec) :: x(:)
logical :: multiple_polylog_converges
logical :: MPL_converges
if(abs(product(x)) < 1) then
multiple_polylog_converges = .true.
MPL_converges = .true.
else
multiple_polylog_converges = .false.
MPL_converges = .false.
end if
END FUNCTION multiple_polylog_converges
END FUNCTION MPL_converges
recursive FUNCTION multiple_polylog(m, x, n_passed) result(res)
recursive FUNCTION MPL(m, x, n_passed) result(res)
! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
integer :: m(:)
complex(kind=prec) :: x(:)
......@@ -61,7 +61,7 @@ CONTAINS
! recursion step
res = 0
do i = 1, n
res = res + x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i - 1)
res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1)
end do
! a nicer way to do it would be but problem is i
......@@ -69,14 +69,13 @@ CONTAINS
! res = sum( x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i(1) - 1) )
end if
END FUNCTION multiple_polylog
END FUNCTION MPL
END MODULE functions
END MODULE mpl_module
! PROGRAM test
! use functions
! use mpl_module
! logical :: result
! result = multiple_polylog_converges( dcmplx((/10.1d0,.7d0,.3d0/)) )
! result = MPL_converges( dcmplx((/10.1d0,.7d0,.3d0/)) )
! print*, result
! end PROGRAM test
......@@ -3,15 +3,20 @@
! These tests assume that GPLInfinity = 30
PROGRAM TEST
use functions
use mpl_module
use gpl_module
implicit none
complex(kind=prec) :: res
real, parameter :: tol = 1.0e-14
call do_multiple_polylog_tests()
call do_MPL_tests()
res = GPL()
CONTAINS
subroutine do_multiple_polylog_tests()
subroutine do_MPL_tests()
integer :: m2(2), m3(3)
integer :: bla
complex(kind=prec) :: x2(2), x3(3)
......@@ -22,27 +27,26 @@ CONTAINS
m2 = (/ 1,1 /)
x2 = cmplx((/ 0.3156498673740053, 0.3431255827785649 /))
ref = cmplx(0.022696600480693277651633)
call check_multiple_polylog(m2,x2,ref)
call check_MPL(m2,x2,ref)
m2 = (/ 1,1 /)
x2 = cmplx((/ 0.03, 0.5012562893380046 /))
ref = cmplx(0.00023134615630308335448329926098409)
call check_multiple_polylog(m2,x2,ref)
call check_MPL(m2,x2,ref)
m3 = (/ 2,1,2 /)
x3 = cmplx((/ 0.03, 0.5012562893380046, 55.3832 /))
ref = cmplx(0.000023446106415452030937059124671151)
call check_multiple_polylog(m3,x3,ref)
end subroutine do_multiple_polylog_tests
call check_MPL(m3,x3,ref)
end subroutine do_MPL_tests
subroutine check_multiple_polylog(m,x,ref)
! checks multiple_polylog(m,x) against reference value ref
subroutine check_MPL(m,x,ref)
! checks MPL(m,x) against reference value ref
integer :: m(:)
complex(kind=prec) :: x(:)
complex(kind=prec) :: res, ref
real(kind=prec) :: delta
res = multiple_polylog(m,x)
res = MPL(m,x)
delta = abs((res-ref)/ref)
if(delta < tol) then
print*, 'passed with delta = ', delta
......@@ -52,9 +56,7 @@ CONTAINS
print*, 'res=',res,'ref=',ref
stop 1
end if
end subroutine check_multiple_polylog
end subroutine check_MPL
END PROGRAM TEST
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment