28 use,
intrinsic :: iso_c_binding, only: c_double, c_int, c_f_pointer
29 use,
intrinsic :: iso_fortran_env, only: error_unit, input_unit, int32, int64, real64, output_unit
31 use blas_lapack_gbl,
only: blasint
32 use phys_const_gbl,
only: pi, imu
33 use precisn_gbl,
only: wp
39 type,
bind(C) :: gsl_sf_result
42 end type gsl_sf_result
46 function gsl_set_error_handler_off () bind(C, name='gsl_set_error_handler_off')
47 use,
intrinsic :: iso_c_binding, only: c_ptr
48 type(c_ptr) :: gsl_set_error_handler_off
49 end function gsl_set_error_handler_off
54 function gsl_sf_lngamma_complex_e (zr, zi, lnr, arg) bind(C, name='gsl_sf_lngamma_complex_e')
55 use,
intrinsic :: iso_c_binding, only: c_int, c_double
57 real(c_double),
value :: zi, zr
58 type(gsl_sf_result) :: lnr, arg
59 integer(c_int) :: gsl_sf_lngamma_complex_e
60 end function gsl_sf_lngamma_complex_e
65 function gsl_sf_coulomb_wave_fg_e (eta, rho, l, k, F, Fp, G, Gp, expF, expG) bind(C, name='gsl_sf_coulomb_wave_FG_e')
66 use,
intrinsic :: iso_c_binding, only: c_int, c_double
68 real(c_double),
value :: eta, rho, l
69 integer(c_int),
value :: k
70 type(gsl_sf_result) :: F, Fp, G, Gp
71 real(c_double) :: expF, expG
72 integer(c_int) :: gsl_sf_coulomb_wave_FG_e
73 end function gsl_sf_coulomb_wave_fg_e
78 pure function gsl_sf_hyperg_u (a, b, x) bind(C, name='gsl_sf_hyperg_U')
79 use,
intrinsic :: iso_c_binding, only: c_double
80 real(c_double),
value :: a, b, x
81 real(c_double) :: gsl_sf_hyperg_U
82 end function gsl_sf_hyperg_u
91 real(wp) :: eta,
cphaz
97 subroutine coulfg (xx, eta1, xlmin, xlmax, fc, gc, fcp, gcp, mode1, kfn, ifail)
99 integer :: mode1, kfn, ifail
100 real(wp) :: xx, eta1, xlmin, xlmax, fc(*), gc(*), fcp(*), gcp(*)
106 subroutine couln (l, Z, ERy, r, U, Up, acc, efx, ierr, nlast, fnorm)
108 integer :: l, ierr, nlast
109 real(wp) :: Z, Ery, r, U, Up, acc, efx, fnorm
115 subroutine decay (k, l, r, U, Up, iwrite)
118 real(wp) :: k, r, U, Up
124 subroutine dgemm (transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
125 import real64, blasint
126 character(len=1) :: transa, transb
127 integer(blasint) :: m, n, k, lda, ldb, ldc
128 real(real64) :: A(lda, *), B(ldb, *), C(ldc, *), alpha, beta
134 subroutine dgemv (trans, m, n, alpha, A, lda, X, incx, beta, Y, incy)
135 import real64, blasint
136 character(len=1) :: trans
137 integer(blasint) :: m, n, incx, incy, lda
138 real(real64) :: A(lda, *), X(*), Y(*), alpha, beta
144 subroutine dsyev (jobz, uplo, n, A, lda, eigs, work, lwork, info)
145 import real64, blasint
146 character(len=1) :: jobz, uplo
147 integer(blasint) :: n, lda, lwork, info
148 real(real64) :: A(lda, *), eigs(*), work(*)
154 subroutine dgetrf (m, n, A, lda, ipiv, info)
155 import real64, blasint
156 integer(blasint) :: m, n, lda, ipiv(*), info
157 real(real64) :: A(lda, *)
163 subroutine dgetrs (trans, n, nrhs, A, lda, ipiv, B, ldb, info)
164 import real64, blasint
165 character(len=1) :: trans
166 integer(blasint) :: n, nrhs, lda, ldb, info, ipiv(*)
167 real(real64) :: A(lda, *), B(ldb, *)
173 subroutine zgemm (transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
174 import real64, blasint
175 character(len=1) :: transa, transb
176 integer(blasint) :: m, n, k, lda, ldb, ldc
177 complex(real64) :: A(lda, *), B(ldb, *), C(ldc, *), alpha, beta
183 subroutine zgemv (trans, m, n, alpha, A, lda, X, incx, beta, Y, incy)
184 import real64, blasint
185 character(len=1) :: trans
186 integer(blasint) :: m, n, incx, incy, lda
187 complex(real64) :: A(lda, *), X(*), Y(*), alpha, beta
193 subroutine zgetrf (m, n, A, lda, ipiv, info)
194 import real64, blasint
195 integer(blasint) :: m, n, lda, ipiv(*), info
196 complex(real64) :: A(lda, *)
202 subroutine zgetrs (trans, n, nrhs, A, lda, ipiv, B, ldb, info)
203 import real64, blasint
204 character(len=1) :: trans
205 integer(blasint) :: n, nrhs, lda, ldb, info, ipiv(*)
206 complex(real64) :: A(lda, *), B(ldb, *)
212 subroutine zgetri (n, A, lda, ipiv, work, lwork, info)
213 import real64, blasint
214 integer(blasint) :: n, lda, ipiv(*), lwork, info
215 complex(real64) :: A(lda, *), work(*)
239 integer,
intent(in) :: nchan, nstat
240 real(wp),
intent(in) :: wmat(:, :), eig(:), Etot
241 real(wp),
intent(inout) :: Rmat(:, :)
242 real(wp),
allocatable :: E_wmat(:, :)
243 integer(blasint) :: lda, ldc, n, m, istat
245 allocate (e_wmat(nchan, nstat))
255 e_wmat(1:nchan, istat) = wmat(1:nchan, istat) / (eig(istat) - etot)
258 call dgemm(
'N',
'T', n, n, m, 0.5_wp, wmat, lda, e_wmat, n, 0.0_wp, rmat, ldc)
277 real(wp),
intent(in) :: K(:, :)
278 complex(wp),
intent(inout) :: T(:, :)
279 integer,
intent(in) :: nopen
281 integer(blasint) :: i, j, n, one = 1, ldk, info
282 integer(blasint),
allocatable :: ipiv(:)
283 real(wp),
allocatable :: A(:, :), ReT(:, :), ImT(:, :)
285 allocate (a(nopen, nopen), ipiv(nopen), ret(nopen, nopen), imt(nopen, nopen))
293 imt(i, j) = 2*k(i, j)
297 call dgemm(
'T',
'N', n, n, n,
rone, k, ldk, k, ldk,
rzero, a, n)
301 a(i, i) = 1 + a(i, i)
304 call dgetrf(n, n, a, n, ipiv, info)
307 print
'(a)',
'ERROR: Failed to factorize 1 + K^2'
311 call dgetrs(
'N', n, n, a, n, ipiv, imt, n, info)
314 print
'(a)',
'ERROR: Failed to back-substitute (1 + K^2) X = 2K'
318 call dgemm(
'T',
'N', n, n, n, -
rone, k, ldk, imt, n,
rzero, ret, n)
323 t(i, j) = cmplx(ret(i, j), imt(i, j), wp)
343 complex(wp),
intent(inout) :: S(:, :)
344 complex(wp),
intent(in) :: T(:, :)
345 integer,
intent(in) :: nopen
351 s(1:nopen, j) = t(1:nopen, j)
352 s(j, j) = 1 + s(j, j)
373 subroutine scatak (rmatr, eig, etot, w, Ek, l, Km, Ak)
377 complex(wp),
intent(inout) :: Ak(:,:)
378 real(wp),
intent(in) :: rmatr, eig(:), etot, Ek(:), w(:,:), Km(:,:)
379 integer,
intent(in) :: l(:)
381 complex(wp),
allocatable :: T(:,:)
382 real(wp),
allocatable :: V(:,:), TT(:,:,:), Sp(:,:), Cp(:,:), P(:), XX(:,:,:), wFp(:,:,:), k(:)
384 integer(blasint) :: i, j, nchan, nstat, nopen, nochf, nochf2, ldwmp
385 real(wp) :: F, Fp, G, Gp, Z = 1
388 nopen = count(ek > 0)
394 nopen = min(nopen, nchan)
395 nochf = min(nopen, nochf)
398 allocate (v(nchan, nchan), t(nopen, nopen), sp(nchan, nchan), cp(nchan, nchan), tt(nchan, 2, nopen), &
399 xx(nchan, 2, nochf), wfp(nstat, 2, nochf), p(nstat))
401 t = 0; v = 0; sp = 0; cp = 0
407 t = t + imu*km(1:nopen, 1:nopen)
409 tt(1:nopen, 1, 1:nochf) = real(t(1:nopen, 1:nochf), wp)
410 tt(1:nopen, 2, 1:nochf) = aimag(t(1:nopen, 1:nochf))
414 call coul(z, l(i), ek(i), rmatr, f, fp, g, gp)
415 v(i,i) = sqrt(2/(pi*k(i)))
421 call dgemm(
'N',
'N', nchan, nopen, nchan,
rone, cp, nchan, km, nchan,
rone, sp, nchan)
422 call dgemm(
'N',
'N', nchan, nochf2, nopen,
rone, sp, nchan, tt, nchan,
rzero, xx, nchan)
423 call dgemm(
'N',
'N', nchan, nochf2, nchan,
rone, v, nchan, xx, nchan,
rzero, tt, nchan)
424 call dgemm(
'T',
'N', nstat, nochf2, nchan,
rone, w, ldwmp, tt, nchan,
rzero, wfp, nstat)
425 p = 0.5 / (eig - etot)
427 ak(:, i) = p * cmplx(wfp(:, 1, i), wfp(:, 2, i), wp)
440 real(wp) function
cphase (z, l, k) result (sigma)
442 integer,
intent(in) :: l
443 real(wp),
intent(in) :: z, k
445 type(gsl_sf_result) :: lnr, arg
446 integer(c_int) :: err
447 real(c_double) :: zr, zi
451 err = gsl_sf_lngamma_complex_e(zr, zi, lnr, arg)
454 sigma =
cphaz(l, -z/k, 0)
468 real(wp),
intent(in) :: x
469 integer,
intent(in) :: l
471 real(wp) :: absval, phase, z = 1
485 g = imu * (-1)**l * pi / (sinh(pi*x) *
complex_gamma(-l - 1, -x))
492 absval = absval * (k**2 + x**2)
494 absval = absval * pi*x/sinh(pi*x)
495 absval = sqrt(absval)
498 phase =
cphase(z, l, -1/x)
501 g = absval * cmplx(cos(phase), sin(phase), wp)
518 subroutine coul (Z, l, Ek, r, F, Fp, G, Gp)
520 integer,
intent(in) :: l
521 real(wp),
intent(in) :: Z, Ek, r
522 real(wp),
intent(inout) :: F, Fp, G, Gp
524 call coul_gsl(z, l, ek, r, f, fp, g, gp)
537 complex(wp) function coulh (Z, s, l, Ek, r)
result(H)
539 integer,
intent(in) :: s, l
540 real(wp),
intent(in) :: z, ek, r
542 real(wp) :: f, fp, g, gp
544 call coul(z, l, ek, r, f, fp, g, gp)
546 h = cmplx(g, s*f, wp)
561 complex(wp) function coulh_asy (Z, s, l, Ek, r)
result(H)
565 integer,
intent(in) :: s, l
566 real(wp),
intent(in) :: z, ek, r
569 complex(wp) :: a, b, term
587 term = term * (a + n - 1) * (b + n - 1) / (2*s*imu*n*k*r)
589 term = term * (a + n - 1) * (b + n - 1) / (-2*s*n*k*r)
595 h = h * exp(imu*s*(k*r + z*log(2*k*r)/k - pi*l/2 +
cphase(z, l, k)))
597 h = h * exp(s*(-k*r + z*log(2*k*r)/k))
613 integer,
intent(in) :: l
614 real(wp),
intent(in) :: Z, Ek, r
615 real(wp),
intent(inout) :: F, Fp, G, Gp
617 type(gsl_sf_result) :: resF, resFp, resG, resGp
618 real(c_double) :: expF, expG, eta, rho, lc, k, a, b, x, U, Up
619 integer(c_int) :: err, zero = 0
627 err = gsl_sf_coulomb_wave_fg_e(eta, rho, lc, zero, resf, resfp, resg, resgp, expf, expg)
638 u = gsl_sf_hyperg_u(a, b, x)
639 up = -2*k*a * gsl_sf_hyperg_u(a + 1, b + 1, x)
642 g = exp(-x/2) * x**(l + 1) * u
643 gp = exp(-x/2) * x**(l + 1) * (-k * u + (l + 1)/r * u + up)
662 integer,
intent(in) :: l
663 real(wp),
intent(in) :: Z, Ek, r
664 real(wp),
intent(inout) :: F, Fp, G, Gp
665 integer :: ifail, mode = 1, kfn = 0, nterms, iwrite = 0
666 real(wp) :: fc(l + 1), gc(l + 1), fcp(l + 1), gcp(l + 1), lc, efx, fnorm, acc = 1e-10_wp, k
672 call coulfg(k*r, -z/k, lc, lc, fc, gc, fcp, gcp, mode, kfn, ifail)
680 call couln(l, z, 2*ek, r, g, gp, acc, efx, ifail, nterms, fnorm)
682 call decay(2*ek, l, r, g, gp, iwrite)
705 real(wp),
intent(in) :: M(:, :)
706 real(wp),
intent(inout) :: eigs(:), eigv(:, :)
707 real(wp),
intent(in),
optional :: check
709 integer(blasint) :: i, n, info, lwork
710 real(wp),
allocatable :: work(:), tmp(:, :)
711 real(wp) :: diag, offdiag
716 call dsyev(
'V',
'U', n, eigv, n, eigs, work, -1_blasint, info)
719 write (error_unit,
'(a,i0,a)')
'Error ', info,
' when preparing matrix diagonalization'
725 allocate (work(lwork))
726 call dsyev(
'V',
'U', n, eigv, n, eigs, work, lwork, info)
729 write (error_unit,
'(a,i0,a)')
'Error ', info,
' when diagonalizing matrix'
735 if (eigv(maxloc(abs(eigv(:, i)), 1), i) < 0)
then
736 eigv(:, i) = -eigv(:, i)
741 if (
present(check))
then
744 call dgemm(
'N',
'T', n, n, n, 1._wp, eigv, n, eigv, n, 0._wp, tmp, n)
747 offdiag = sqrt(sum(abs(tmp(:, i))**2) - diag**2)
748 if (abs(diag - 1) > check .or. abs(offdiag) > check)
then
749 write (error_unit,
'(a,i0,2e25.15)')
'WARNING: Orthogonality error in symmetric diagonalization ', &
768 complex(wp),
allocatable,
intent(inout) :: T(:, :)
770 complex(wp),
allocatable :: work(:)
771 complex(wp) :: nwork(1)
773 integer(blasint),
allocatable :: ipiv(:)
774 integer(blasint) :: m, info, lwork
780 call zgetrf(m, m, t, m, ipiv, info)
784 call zgetri(m, t, m, ipiv, nwork, lwork, info)
785 lwork = int(nwork(1))
786 allocate (work(lwork))
789 call zgetri(m, t, m, ipiv, work, lwork, info)
803 complex(wp),
allocatable,
intent(inout) :: A(:, :)
804 real(wp),
allocatable,
intent(inout) :: X(:, :), Y(:, :)
805 integer,
intent(in) :: n
807 integer(blasint) :: m, info, nrhs, lda
808 integer(blasint),
allocatable :: ipiv(:)
809 complex(wp),
allocatable :: XX(:)
816 call zgetrf(m, m, a, lda, ipiv, info)
819 xx = cmplx(x(1:n, 1), x(1:n, 2), wp)
823 call zgetrs(
'N', m, nrhs, a, lda, ipiv, xx, m, info)
825 y(1:n, 1) = real(xx, wp)
826 y(1:n, 2) = aimag(xx)
845 integer,
intent(inout) :: p(:)
846 integer :: i, j, k, n
876 if (p(i) < p(k))
then
878 do while (p(i) >= p(j))
881 call swap(p(i), p(j))
885 else if (i == 1)
then
902 integer,
intent(inout) :: a, b
920 integer,
intent(inout) :: a(:)
926 call swap(a(i), a(n + 1 - i))
948 complex(wp),
intent(inout) :: X, err
949 complex(wp),
intent(in) :: dX
992 use dipelm_special_functions,
only: threej
995 integer,
intent(in) :: j, n, p(n), li, mi, qi(n), lj, mj, qj(n)
996 integer :: mu, l(n + 1), lp(n + 1), a(n + 1), ap(n + 1), b(n + 1), bp(n + 1)
1000 lp(1) = lj; lp(2:) = 1
1002 b(1) = -mi; b(2:) = qi
1003 bp(1) = -mj; bp(2:) = qj
1009 a(1) = -mu; a(2:) = p
1010 ap(1) = -mu; ap(2:) = p
1013 q = (-1)**mu * (2*j + 1) * sqrt((2*li +
rone)*(2*lj +
rone)) &
1014 * threej(2*li, 0, 2*lj, 0, 2*j, 0) * threej(2*li, 2*mu, 2*lj, -2*mu, 2*j, 0)
1054 use dipelm_special_functions,
only: threej
1057 integer,
intent(in) :: n
1058 integer,
intent(in) :: l(n), a(n), b(n), lp(n), ap(n), bp(n)
1060 integer :: u(n - 1), v(n - 1), up(n - 1), vp(n - 1), j(n - 1), jp(n - 1), j1, jp1
1061 real(wp) :: wa, wb, wap, wbp, w
1067 if (l(1) == lp(1) .and. a(1) == ap(1) .and. b(1) == bp(1)) k = 1 / (2*l(1) +
rone)
1071 u(1) = a(1) + a(2); u(2:) = a(3:)
1072 v(1) = b(1) + b(2); v(2:) = b(3:)
1074 up(1) = ap(1) + ap(2); up(2:) = ap(3:)
1075 vp(1) = bp(1) + bp(2); vp(2:) = bp(3:)
1081 do j1 = max(abs(u(1)), max(abs(v(1)), abs(l(1) - l(2)))), l(1) + l(2)
1082 do jp1 = max(abs(up(1)), max(abs(vp(1)), abs(lp(1) - lp(2)))), lp(1) + lp(2)
1087 wa = threej(2*l(1), 2*a(1), 2*l(2), 2*a(2), 2*j(1), -2*u(1))
1088 wb = threej(2*l(1), 2*b(1), 2*l(2), 2*b(2), 2*j(1), -2*v(1))
1089 wap = threej(2*lp(1), 2*ap(1), 2*lp(2), 2*ap(2), 2*jp(1), -2*up(1))
1090 wbp = threej(2*lp(1), 2*bp(1), 2*lp(2), 2*bp(2), 2*jp(1), -2*vp(1))
1092 k = k + (-1)**(u(1) - v(1) + up(1) - vp(1)) * (2*j1 + 1) * (2*jp1 + 1) * wa * wap * wb * wbp * w
1122 real(wp) function
beta_2p_arb_pol_sph_harm (lbig, mbig, n, pb, lp, mp, qb, pa, l, m, qa) result (t)
1124 use dipelm_special_functions,
only: threej
1125 use dipelm_defs,
only: pi
1128 integer,
intent(in) :: lbig, mbig, n, pa(n), pb(n), lp, mp, qa(n), l, m, qb(n)
1129 integer :: k1, k2, m1, m2, mp1, mp2, mu
1130 real(wp) :: i, i_el, fac
1133 print *,
'Not implemented for orders other than 2', n
1137 mp1 = -(-pa(1) + pb(1))
1138 mp2 = -(-pa(2) + pb(2))
1142 if ( mbig /= -(mp1 + mp2) )
return
1149 fac = sqrt((2*l +
rone)*(2*lp +
rone)*(2*lbig +
rone) / (4*pi))
1154 i_el = (-1)**(mp) * fac * threej(2*l, 2*m, 2*lp, -2*mp, 2*lbig, 2*mu) * threej(2*l, 0, 2*lp, 0, 2*lbig, 0)
1160 i = (-1)**(pa(1) + pa(2)) * threej(2*1, 2*(-pa(1)), 2*1, 2*pb(1), 2*k1, 2*mp1) * &
1161 &threej(2*1, 2*(-pa(2)), 2*1, 2*pb(2), 2*k2, 2*mp2)
1164 m1 = -(-qa(1) + qb(1))
1165 m2 = -(-qa(2) + qb(2))
1167 i = i * (2*k1 + 1) * (2*k2 + 1) * (-1)**(-qa(1)-qa(2)) * threej(2*1, 2*(-qa(1)), 2*1, 2*qb(1), 2*k1, 2*m1) * &
1168 threej(2*1, 2*(-qa(2)), 2*1, 2*qb(2), 2*k2, 2*m2)
1171 i = i * i_el * threej(2*k1, 2*mp1, 2*k2, 2*mp2, 2*lbig, 2*mbig) * threej(2*k1, 2*m1, 2*k2, 2*m2, 2*lbig, 2*mu)
1190 real(wp) function
beta_2p_demekhin (l, p1, p2, li, mi, q1i, q2i, lj, mj, q1j, q2j) result (b)
1192 use dipelm_special_functions,
only: threej
1195 integer,
intent(in) :: l, p1, p2, li, mi, lj, mj, q1i, q1j, q2i, q2j
1197 integer :: m_l, j, m_j, k, m_k
1207 w(1) = threej(2*1, 2*p1, 2*1, -2*p1, 2*j, 0)
1208 w(2) = threej(2*1, 2*q1i, 2*1, -2*q1j, 2*j, -2*m_j)
1209 w(3) = threej(2*1, 2*q2i, 2*1, -2*q2j, 2*k, -2*m_k)
1210 w(4) = threej(2*1, 2*p2, 2*1, -2*p2, 2*k, 0)
1211 w(5) = threej(2*j, 2*m_j, 2*k, 2*m_k, 2*l, 2*m_l)
1212 w(6) = threej(2*j, 0, 2*k, 0, 2*l, 0)
1213 b = b + (-1)**(q1j + q2j) * (2*j + 1) * (2*k + 1) * product(w(1:6))
1219 w(1) = threej(2*li, 0, 2*lj, 0, 2*l, 0)
1220 w(2) = threej(2*li, -2*mi, 2*lj, 2*mj, 2*l, -2*m_l)
1222 b = b * (-1)**mi * (2*l + 1) * sqrt((2*li +
rone)*(2*lj +
rone)) * product(w(1:2))
1250 integer,
intent(in) :: q(:)
1251 integer :: i, n, iq(size(q))
1263 else if (mod(n, 2) == 1)
then
1270 if (q(1) == q(i))
then
Hard-coded parameters of MULTIDIP.
complex(wp), parameter cone
real(wp), parameter rzero
integer, parameter ntermsasy
Number of terms in expansion of the Coulomb-Hankel functions.
complex(wp), parameter czero
Special functions and objects used by MULTIDIP.
subroutine diagonalize_symm_matrix(M, eigs, eigv, check)
Diagonalize a real symmetric matrix.
real(wp) function beta_2p_arb_pol_sph_harm(Lbig, Mbig, n, pB, lp, mp, qB, pA, l, m, qA)
Two-photon asymmetry parameter for the arbitrary polarized case.
recursive real(wp) function wigner_d_multiint(n, l, a, b, lp, ap, bp)
Orientation averaging.
subroutine solve_complex_system(n, A, X, Y)
Solve system of equations with complex matrix.
subroutine kahan_add(X, dX, err)
Compensated summation.
recursive complex(wp) function complex_gamma(l, x)
Complex Gamma function.
subroutine coul_ukrmol(Z, l, Ek, r, F, Fp, G, Gp)
Coulomb functions (UKRmol)
complex(wp) function coulh(Z, s, l, Ek, r)
Coulomb-Hankel function.
subroutine coul(Z, l, Ek, r, F, Fp, G, Gp)
Coulomb functions.
complex(wp) function coulh_asy(Z, s, l, Ek, r)
Coulomb-Hankel function (asymptotic form)
subroutine calculate_t_matrix(K, T, nopen)
Calculate T-matrix from K-matrix.
subroutine scatak(rmatr, eig, etot, w, Ek, l, Km, Ak)
Photoionization coefficient.
subroutine calculate_s_matrix(T, S, nopen)
Calculate S-matrix from T-matrix.
real(wp) recursive function cartesian_vector_component_product_average(q)
Return Cartesian invariant.
subroutine invert_matrix(T)
Invert a complex matrix.
subroutine coul_gsl(Z, l, Ek, r, F, Fp, G, Gp)
Coulomb functions (GSL)
real(wp) function beta_contraction_tensor(J, n, p, li, mi, qi, lj, mj, qj)
Angular part of the beta parameters.
subroutine swap(a, b)
Exchange value of two integers.
subroutine calculate_r_matrix(nchan, nstat, wmat, eig, Etot, Rmat)
Calculate scattering R-matrix.
subroutine reverse(a)
Reverse order of elements in array.
real(wp) function cphase(Z, l, k)
Coulomb phase.
real(wp) function beta_2p_demekhin(L, p1, p2, li, mi, q1i, q2i, lj, mj, q1j, q2j)
Two-photon asymmetry parameter.
logical function next_permutation(p)
Construct or advance permutation.