MODULE constants INTEGER, PARAMETER :: dp = KIND(1.0D0) INTEGER, PARAMETER :: dpc = KIND((1.0D0,1.0D0)) REAL(DP), PUBLIC, PARAMETER :: pi = 3.141592741012573_dp REAL(DP), PUBLIC, PARAMETER :: pi_2 = 1.570796370506287_dp REAL(DP), PUBLIC, PARAMETER :: pi_4 = 0.7853981852531433_dp END MODULE constants ! ! Routines for various wave functions are also included ! MODULE wave_functions USE constants CONTAINS ! ! H.O. functions using Kummers function ! REAL(DP) FUNCTION rnl(n,l,z) IMPLICIT NONE INTEGER :: lll, nn INTEGER, INTENT(IN) :: l, n REAL(DP) :: y, dl, gamfaa, dfll, gamfab, dfnn REAL(DP), INTENT(IN) :: z rnl=0. ; y=0.5_dp*z*z IF(y > 60.0_dp) RETURN dl = l IF((ABS(z) < 1.0d-6) .AND. (l == 0)) rnl = 1.0_dp IF( ABS(z) > 1.0d-6) rnl = (z**l) * EXP(-y) * hypkum(n,dl+1.5_dp,z*z) gamfaa = 0.5_dp * SQRT(pi) IF(l /= 0) THEN DO lll = 1, l dfll = lll - 1 gamfaa = gamfaa * (dfll + 1.5_dp) ENDDO ENDIF gamfab = gamfaa IF(n /= 0) THEN dfll = dl + 0.5_dp DO nn = 1, n dfnn = nn gamfab = gamfab * ((dfnn + dfll) / dfnn) ENDDO ENDIF rnl = rnl * (SQRT(2.0_dp * gamfab) / gamfaa) END FUNCTION rnl ! ! Kummers function, Abramowitz & Stegun ! exp. 13.1.2. a(there) equals (-n) ! REAL(DP) FUNCTION hypkum(n,b,z) IMPLICIT NONE INTEGER :: nmax, nf INTEGER, INTENT(IN) :: n REAL(DP) :: af, bf, zf, term, dfnf, xadd, sum REAL(DP), INTENT(IN) :: b, z IF(n < 0) WRITE (6,*)' error exit in hypkum ', n,b,z hypkum = 1.0_dp IF(n == 0) RETURN nmax = n ; af = - n ; bf = b ; zf = z ; sum = 1.0 ; term = 1.0_dp DO nf = 1, nmax dfnf = nf xadd = dfnf - 1.0_dp term = term * ((af + xadd) / (bf + xadd)) * (zf / dfnf) IF(ABS(term) < 1.0d-12) EXIT sum = sum + term ENDDO hypkum = sum END FUNCTION hypkum ! This function sets up the recursive relation ! for the associated Legendre polynomials REAL(DP) FUNCTION legendre_polynomials(l, m, x) IMPLICIT NONE REAL(DP) :: fact,pll,pmm,pmmp1,somx2 REAL(DP), INTENT(IN) :: x INTEGER :: i,ll INTEGER, INTENT(IN) :: l, m ! check whether m, l and x are ok IF((M < 0).OR.(M > L).OR.(ABS(X) > 1.)) THEN WRITE(6,*) 'bad arguments', m, l, x; RETURN ENDIF ! calculate now pmm as starting point for iterations pmm=1.0 IF (m > 0) THEN somx2=SQRT((1.0-x)*(1.0+x)) fact=1.0_dp; DO i=1, m pmm = -fact*somx2*pmm fact = fact+2.0_dp ENDDO ENDIF ! if l == m we do not need to use recursion relation IF (l == m) THEN legendre_polynomials=pmm ! recursive relation for associated Legendre polynomials ELSE pmmp1=x*(2*m+1)*pmm ! analytical formula for the case l == m+1 IF (l == (m+1)) THEN legendre_polynomials=pmmp1 ELSE DO ll=m+2, l pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m) pmm=pmmp1 pmmp1=pll ENDDO legendre_polynomials= pll ENDIF ENDIF END FUNCTION legendre_polynomials ! ! ! This routine calculates gauss-legendre mesh points and weights ! input: ! x1 : lower limit of the integration interval ! x2 : upper limit ---------- "" ------------- ! n : the desired number of mesh points ! output : ! x : gauss-legendre mesh points on the interval (x1,x2) ! w : the corresponding weights ! From : Numerical recipes ! F90 version : M. Hjorth-Jensen ! SUBROUTINE gauss_legendre(x1,x2,x,w,n) IMPLICIT NONE INTEGER, INTENT(IN) :: n INTEGER :: i, j, m REAL(DP), INTENT(IN) :: x1, x2 REAL(DP), INTENT(INOUT) :: x, w REAL(DP) :: eps DIMENSION :: x(n), w(n) PARAMETER (eps=3.D-14) REAL(DP) :: p1,p2,p3,pp,xl,xm,z,z1 m=(n+1)/2 xm=0.5_dp*(x2+x1) xl=0.5_dp*(x2-x1) DO i=1,m z1=0. z=COS(pi*(i-.25_dp)/(n+.5_dp)) DO WHILE ( ABS(z-z1) > EPS) p1=1.0_dp p2=0.0_dp DO j=1,n p3=p2 p2=p1 p1=((2.0_dp*j-1.)*z*p2-(j-1.0_dp)*p3)/j ENDDO pp=n*(z*p1-p2)/(z*z-1.) z1=z z=z-p1/pp ENDDO x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.0_dp*xl/((1.0_dp-z*z)*pp*pp) w(n+1-i)=w(i) ENDDO END SUBROUTINE gauss_legendre SUBROUTINE laguerre_general( n, alpha, x, cx ) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp ) :: alpha REAL ( dp ) :: cx(0:n) INTEGER :: i REAL ( dp ), INTENT(IN) :: x IF ( alpha <= -1.0D+00 ) THEN WRITE ( *, '(a)' ) ' ' WRITE ( *, '(a)' ) 'LAGUERRE_GENERAL - Fatal error!' WRITE ( *, '(a,g14.6)' ) ' The input value of ALPHA is ', alpha WRITE ( *, '(a)' ) ' but ALPHA must be greater than -1.' STOP END IF IF ( n < 0 ) THEN RETURN END IF cx(0) = 1.0D+00 IF ( n == 0 ) THEN RETURN END IF cx(1) = 1.0D+00 + alpha - x DO i = 2, n cx(i) = ( ( REAL ( 2 * i - 1, kind = 8 ) + alpha - x ) * cx(i-1) & + ( REAL ( - i + 1, kind = 8 ) - alpha ) * cx(i-2) ) & / REAL ( i, kind = 8 ) END DO END SUBROUTINE laguerre_general SUBROUTINE laguerre_complex( n, alpha, x, cx ) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL (dp ) :: alpha COMPLEX ( DPC ) :: cx(0:n) INTEGER :: i COMPLEX ( DPC ), INTENT(IN) :: x IF ( alpha <= -1.0D+00 ) THEN WRITE ( *, '(a)' ) ' ' WRITE ( *, '(a)' ) 'LAGUERRE_GENERAL - Fatal error!' WRITE ( *, '(a,g14.6)' ) ' The input value of ALPHA is ', alpha WRITE ( *, '(a)' ) ' but ALPHA must be greater than -1.' STOP END IF IF ( n < 0 ) THEN RETURN END IF cx(0) = 1.0D+00 IF ( n == 0 ) THEN RETURN END IF cx(1) = 1.0D+00 + alpha - x DO i = 2, n cx(i) = ( ( REAL ( 2 * i - 1, kind = 8 ) + alpha - x ) * cx(i-1) & + ( REAL ( - i + 1, kind = 8 ) - alpha ) * cx(i-2) ) & / REAL ( i, kind = 8 ) END DO END SUBROUTINE laguerre_complex DOUBLE PRECISION FUNCTION fac(m) IMPLICIT NONE INTEGER, INTENT(IN) :: m INTEGER :: i fac = 0.0D0 IF(m == 0) RETURN DO i=1,m fac=fac+LOG(FLOAT(i)) ENDDO END FUNCTION fac DOUBLE PRECISION FUNCTION dfac(m) IMPLICIT NONE INTEGER, INTENT(IN) :: m INTEGER :: i IF (MOD(m,2).NE.1) STOP 'wrong argument to dfac' dfac = 0.0D0 IF (m == 1)RETURN DO i=3,m,2 dfac=dfac+LOG(FLOAT(i)) ENDDO END FUNCTION dfac END MODULE wave_functions