18 #include "../base/base_uses.f90"
22 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'whittaker'
52 INTEGER,
INTENT(IN) :: n, l2, l1
53 REAL(kind=
dp),
INTENT(IN) :: alpha
54 REAL(kind=
dp),
DIMENSION(n),
INTENT(IN) :: erfa, expa, r
55 REAL(kind=
dp),
DIMENSION(n),
INTENT(OUT) :: wc
62 IF (mod(l, 2) /= 0)
THEN
63 cpabort(
"Total angular momentum has to be even")
75 wc(i) = x**l1*(x**2/(3._dp + y) - alpha*x**4/(5._dp + y) + &
76 alpha**2*x**6/(14._dp + 2._dp*y) - &
77 alpha**3*x**8/(54._dp + 6._dp*y) + &
78 alpha**4*x**10/(256._dp + 24._dp*y) - &
79 alpha**5*x**12/120._dp/(13._dp + y))
83 wc(i) = wc(i) + expa(i)*x**(2*k + 1)*t1**(2*k + 3)* &
86 wc(i) = -wc(i)/2._dp**(l/2 + 2)/t1**(l + 5)/x**(l2 + 1)
111 SUBROUTINE whittaker_c0(wc, r, expa, erfa, alpha, l, n)
112 INTEGER,
INTENT(IN) :: n, l
113 REAL(kind=
dp),
INTENT(IN) :: alpha
114 REAL(kind=
dp),
DIMENSION(n),
INTENT(IN) :: erfa, expa, r
115 REAL(kind=
dp),
DIMENSION(n),
INTENT(OUT) :: wc
118 REAL(
dp) :: t1, t10, t11, t12, t13, t14, t16, t17, t18, t19, t2, t21, t22, t23, t25, t28, &
119 t3, t30, t31, t34, t36, t39, t4, t41, t44, t45, t46, t5, t50, t51, t52, t56, t6, t61, t7, &
122 IF (mod(l, 2) /= 0)
THEN
123 cpabort(
"Angular momentum has to be even")
137 wc(i) = wc(i) + expa(i)*x**(2*k + 1)*t1**(2*k + 3)* &
138 dfac(l + 1)/
dfac(2*k + 1)*2**(k + 1)
140 wc(i) = -wc(i)/2._dp**(l/2 + 2)/t1**(l + 5)
153 t18 = -1._dp/t2/t1*(2._dp*x*t7*t1 - t11*t13)/4._dp
168 t25 = -1._dp/t3/t1*(4._dp*t6*x*t9*t2*t1 + 6._dp*x*t9*t1 - 3*t17*t19)/8._dp
185 t31 = -1._dp/t4/t3*(8._dp*t8*x*t11*t4*t1 + 20._dp*t7*x*t11*t3 + 30._dp*x*t11*t1 - &
186 15._dp*t23*t25)/16._dp
205 t39 = -(16._dp*t3*t2*t6*t11*t8 + 56._dp*t3*x*t6*t10*t17 + 140._dp*t2*t6*t10*t8 + &
206 210._dp*x*t6*t17 - 105._dp*t28*t30*alpha)/t11/t17/32._dp
227 t45 = -(32._dp*t3*x*t6*t12*t10 + 144._dp*t2*t16*t6*t12*t8 + 504._dp*t2*x*t6*t11*t10 + &
228 1260._dp*t16*t6*t28 + 1890._dp*x*t6*t10 - 945._dp*t34*t36*alpha)/t12/t28/64._dp
250 t50 = -(64._dp*t4*t2*t7*t13*t12 + 352._dp*t4*x*t7*t13*t19 + &
251 1584._dp*t3*t2*t7*t13*t9 + 5544._dp*t3*x*t7*t30 + &
252 13860._dp*t2*t7*t12 + 20790._dp*x*t7*t19 - 10395._dp*t39*t41*alpha)/ &
277 t56 = -(128._dp*t4*t3*t7*t13*t14 + 832._dp*t4*t18*t7*t14*t21 + &
278 4576._dp*t4*x*t7*t14*t11 + 20592._dp*t2*t18*t7*t14*t9 + 72072._dp*t3*t7*t13 + &
279 180180._dp*t18*t7*t21 + 270270._dp*x*t7*t11 - 135135._dp*t44*t46*alpha)/ &
305 t61 = -(256._dp*t5*t4*t8*t14*t10 + 1920._dp*t5*t18*t8*t13*t22 + &
306 12480._dp*t5*t2*t8*t13*t28 + 68640._dp*t5*x*t8*t13*t21 + &
307 308880._dp*t4*t8*t13*t10 + 1081080._dp*t18*t8*t22 + &
308 2702700._dp*t2*t8*t28 + 4054050._dp*x*t8*t21 - &
309 2027025._dp*t50*t52*alpha)/t14/t21/512._dp
315 END SUBROUTINE whittaker_c0
340 INTEGER,
INTENT(IN) :: n, l
341 REAL(kind=
dp),
INTENT(IN) :: alpha
342 REAL(kind=
dp),
DIMENSION(n),
INTENT(IN) :: expa, r
343 REAL(kind=
dp),
DIMENSION(n),
INTENT(OUT) :: wc
346 REAL(
dp) :: t1, t10, t13, t14, t17, t18, t2, t21, &
347 t25, t29, t3, t30, t33, t4, t5, t6, &
350 IF (mod(l, 2) /= 0)
THEN
351 cpabort(
"Angular momentum has to be even")
364 wc(i) = wc(i) + alpha**k*x**(2*k)*
fac(l/2)/
fac(k)
366 wc(i) = 0.5_dp*wc(i)/alpha**(l/2 + 1)*expa(i)
375 t6 = 1._dp/alpha*t4/2._dp
387 t9 = t3*(t2 + 1)/t6/2._dp
400 t13 = t3*(t4*t5 + 2._dp*t2 + 2._dp)/t5/alpha/2._dp
414 t17 = t3*(t4*t1*t6*alpha + 3._dp*t4*t6 + 6._dp*t2 + 6._dp)/t14/2._dp
429 t21 = t3*(t5*t7 + 4._dp*t4*t1*t6*alpha + 12._dp*t4*t6 + 24._dp*t2 + 24._dp)/t7/alpha/2._dp
444 t25 = t3*(t5*t1*t8*alpha + 5._dp*t5*t8 + 20._dp*t4*t1*t7*alpha + 60._dp*t4*t7 + &
445 120._dp*t2 + 120._dp)/t8/t7/2._dp
461 t29 = t3*(t5*t4*t8*t7 + 6._dp*t5*t1*t8*alpha + 30._dp*t5*t8 + 120._dp*t4*t1*t18 + &
462 360._dp*t4*t7 + 720._dp*t2 + 720._dp)/t8/t18/2._dp
480 t33 = t3*(t6*t5*t10*t9 + 7*t6*t4*t10*t8 + 42._dp*t6*t1*t10*alpha + &
481 210._dp*t6*t10 + 840._dp*t5*t9 + 2520._dp*t4*t8 + 5040._dp*t2 + 5040._dp)/t30/2._dp
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), parameter, public rootpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Calculates special integrals.
real(kind=dp), parameter epsilon
subroutine, public whittaker_c0a(wc, r, expa, erfa, alpha, l1, l2, n)
int(y^(2+l1+l2) * exp(-alpha*y*y),y=0..x) / x^(l2+1); wc(:) :: output r(:) :: coordinate expa(:) :: e...
subroutine, public whittaker_ci(wc, r, expa, alpha, l, n)
int(y^(l+1) * exp(-alpha*y*y),y=x..infinity);