mpl_module.f90 1.07 KB
Newer Older
Luca's avatar
Luca committed
1 2
 
MODULE mpl_module
3 4
  use globals
  use utils
5
  use maths_functions
6
  implicit none
Luca's avatar
Luca committed
7

8 9
CONTAINS 

10
  FUNCTION MPL_converges(m,x)
Luca's avatar
Luca committed
11 12
    ! checks if the MPL converges 
    complex(kind=prec) :: x(:)
13
    integer :: m(:)
Luca's avatar
Luca committed
14
    logical :: MPL_converges
15
    MPL_converges = .false.
Luca's avatar
Luca committed
16
    if(abs(product(x)) < 1) then
17
      if(m(1) /= 1 .or. abs(x(1) - 1) < zero) then
18 19
        MPL_converges = .true.
      end if
Luca's avatar
Luca committed
20
    end if
Luca's avatar
Luca committed
21
  END FUNCTION MPL_converges
Luca's avatar
Luca committed
22

ulrich_y's avatar
ulrich_y committed
23
  FUNCTION MPL(m, x) result(res)
24
    integer :: m(:)
Luca's avatar
Luca committed
25 26
    complex(kind=prec) :: x(:)
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
27 28 29
    complex(kind=prec) :: t(size(x))
    integer :: q, j, k
    j = size(x)
30

Luca's avatar
Luca committed
31

32 33 34
    if(size(m) /= size(x)) then
      print*, 'Error: m and x must have the same length'
    end if
ulrich_y's avatar
ulrich_y committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
    res=0
    q=0
    t=0
    do while(.true.)
      res = t(1)
      q = q+1
      t(j) = t(j) + x(j)**q/q**m(j)
      do k=1,j-1
        t(j-k) = t(j-k) + t(j-k+1) * x(j-k)**(k+q)/(k+q)**m(j-k)
      enddo

      if (mod(q,2) .eq. 1) then
        if (abs(t(1)-res).lt.MPLdel) exit
      endif
    enddo
    res = t(1)

Luca Naterop's avatar
Luca Naterop committed
52

Luca's avatar
Luca committed
53
  END FUNCTION MPL
Luca's avatar
Luca committed
54

Luca's avatar
Luca committed
55
END MODULE mpl_module