Commit 2e8dc74b authored by Luca Naterop's avatar Luca Naterop

some cleanup

parent 281ece3b
......@@ -7,18 +7,6 @@ MODULE mpl_module
CONTAINS
FUNCTION polylog_series(m,x,n_passed) result(res)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m
integer, optional :: n_passed
complex(kind=prec) :: x, res
integer :: i,n
integer, allocatable :: j(:)
n = merge(n_passed,GPLInfinity,present(n_passed))
j = (/(i, i=1,n,1)/)
res = sum(x**j / j**m)
END FUNCTION polylog_series
FUNCTION MPL_converges(m,x)
! checks if the MPL converges
complex(kind=prec) :: x(:)
......@@ -45,30 +33,24 @@ CONTAINS
if(size(m) /= size(x)) then
print*, 'Error: m and x must have the same length'
end if
res = 0
if(size(m) == 1) then
! base case
res = polylog_series(m(1),x(1), n)
! base case, we get a polylog
do i = 1, n
if(i**m(1) < 0) return ! roll over
if(abs(x(1)**i) < 1.e-250) return
res = res + x(1)**i / i**m(1)
end do
else
! recursion step
res = 0
do i = 1, n
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
! i = (/(j, j=1,n, 1)/)
! res = sum( x(1)**i / i**m(1) * MPL(m(2:), x(2:), i(1) - 1) )
! we could get around this problem by rewriting MPL to operate on each i and returning an array
end if
END FUNCTION MPL
END MODULE mpl_module
! PROGRAM test
! use mpl_module
! logical :: result
! result = MPL_converges( dcmplx((/1.0d0,.7d0,.3d0/)), (/ 1,2,1 /) )
! print*, result
! end PROGRAM test
......@@ -19,7 +19,6 @@ PROGRAM TEST
call do_GPL_tests()
! call do_shuffle_tests() ! put this somewhere else
if(tests_successful) then
print*, 'All tests passed. '
else
......@@ -133,8 +132,8 @@ CONTAINS
ref = cmplx((0.4925755847450199,2.6389214054743295))
call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
! ref = cmplx((-0.0015317713178859967,-0.00045003911367000565))
! call test_one_flat(cmplx([0.,0.,(1-sqrt(1-z**2))/(xchen*z), 1./(xchen*z),1.]),ref,'3.7')
ref = cmplx((-0.0015317713178859967,-0.00045003911367000565))
call test_one_flat(cmplx([0.,0.,(1-sqrt(1-z**2))/(xchen*z), 1./(xchen*z),1.]),ref,'3.7')
! here the chen integral tests start
......@@ -171,7 +170,7 @@ CONTAINS
call test_one_flat(cmplx([0.01, -1.0, 0.01, 1.]),ref,'4.10')
ref = (-0.012539108315054982, -0.015414250168437678)
call test_one_flat(cmplx([0.01, 199.99499987499374, 0.01, 1.]),ref,'4.10')
call test_one_flat(cmplx([0.01, 199.99499987499374, 0.01, 1.]),ref,'4.11')
......
......@@ -176,36 +176,3 @@ CONTAINS
! end subroutine warn
END MODULE utils
! PROGRAM test
! use globals
! use utils
! implicit none
! complex(kind=prec) :: z(4)
! integer :: m_prime(4), condensed_size
! z = cmplx((/0.0,1.7,0.0,0.0/))
! ! transform to condensed notation
! m_prime = get_condensed_m(z)
! print*, abs(z)
! m_prime = get_condensed_m(z)
! print*, abs(z)
! if(find_first_zero(m_prime) == -1) then
! condensed_size = size(m_prime)
! else
! condensed_size = find_first_zero(m_prime)-1
! end if
! print*, condensed_size
! END PROGRAM test
! Contains some functions that might be useful later
! Write your own print function with ability to suppress print
! Muss immer alle prints und warnings ausschalten können
! Test Programm schreiben mit exit codes -> gfortran 'test.f90' und dann 'echo $?'
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