25 use,
intrinsic :: iso_fortran_env
129 real(real64),
intent(in) :: x, y
130 real(real64),
optional,
intent(in) :: rtol
131 real(real64),
optional,
intent(in) :: atol
132 real(real64) :: atol_, rtol_
142 complex(real64),
intent(in) :: x, y
143 real(real64),
optional,
intent(in) :: rtol
144 real(real64),
optional,
intent(in) :: atol
145 real(real64) :: atol_, rtol_
158 integer,
intent(in) :: dim
159 integer,
intent(in) :: diag
160 integer :: matrix(dim, dim)
166 matrix(ii, ii) = diag
172 recursive function hermite(n, x)
result (h)
173 integer,
intent(in) :: n
174 real(real64),
intent(in) :: x
192 recursive function factorial (n)
RESULT (fac)
193 integer,
intent(in) :: n
209 real(real64),
intent(in) :: xx(3)
210 integer,
intent(in) :: li, mi
211 complex(real64),
intent(out) :: ylm
214 real(real64) :: dx, dy, dz, r, plm, cosm, sinm, cosmm1, sinmm1, cosphi, sinphi
215 real(real64),
parameter :: tiny = 1.e-15_real64
225 r =
sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
242 r =
sqrt(dx*dx + dy*dy)
256 cosm = cosmm1*cosphi - sinmm1*sinphi
257 sinm = cosmm1*sinphi + sinmm1*cosphi
261 ylm = plm*cmplx(cosm, sinm, real64)
280 real(real64),
intent(in) :: xx(3)
281 integer,
intent(in) :: li, mi
282 real(real64),
intent(out) :: ylm
284 integer,
parameter :: lmaxd = 20
285 real(real64),
parameter :: tiny = 1.e-15_real64
286 integer :: i, ilm0, l, m, mabs
287 integer,
save :: lmax = -1
289 real(real64) :: cmi, cosm, cosmm1, cosphi, dphase, dplg, fac, &
290 fourpi, plgndr, phase, pll, pmm, pmmp1, sinm, &
291 sinmm1, sinphi, rsize, rx, ry, rz, xysize
292 real(real64),
save :: c(0:(lmaxd+1)*(lmaxd+1))
303 do i = l - m + 1, l + m
306 c(ilm0 + m) =
sqrt(fac)
308 if (m /= 0) c(ilm0 + m) = c(ilm0 + m)*
sqrt(
m_two)
309 c(ilm0 - m) = c(ilm0 + m)
322 rsize =
sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
323 if (rsize < tiny)
then
348 ylm = c(4)*6.0_real64*rx*ry
356 ylm = c(8)*
m_three*(rx*rx - ry*ry)
363 xysize =
sqrt(max(rx*rx + ry*ry, tiny))
371 cosm = cosmm1*cosphi - sinmm1*sinphi
372 sinm = cosmm1*sinphi + sinmm1*cosphi
380 dphase = (-mabs)*sinm
388 pmm = (-pmm)*fac*xysize
395 dplg = (-li)*rz*pmm/(xysize**2)
397 pmmp1 = rz*(2*mabs + 1)*pmm
398 if (li == mabs + 1)
then
400 dplg = -((li*rz*pmmp1 - (mabs + li)*pmm)/(xysize**2))
403 pll = (rz*(2*l - 1)*pmmp1 - (l + mabs - 1)*pmm)/(l - mabs)
408 dplg = -((li*rz*pll - (l + mabs - 1)*pmm)/(xysize**2))
414 ylm = cmi*plgndr*phase
429 subroutine weights(N, M, cc, side)
430 integer,
intent(in) :: n, m
431 real(real64),
intent(out) :: cc(0:,0:,0:)
432 integer,
optional,
intent(in) :: side
434 integer :: i, j, k, mn, side_
435 real(real64) :: c1, c2, c3, c4, c5, xi
436 real(real64),
allocatable :: x(:)
440 safe_allocate(x(0:m))
442 if (
present(side))
then
452 x(:) = (/(-i,i=0,mn)/)
456 x(:) = (/(i,i=0,mn)/)
460 x(:) = (/0,(-i,i,i=1,mn)/)
483 if (j <= n) cc(k, j - 1, j) =
m_zero
484 cc(k, j, 0) = c4*cc(k, j - 1, 0)/c3
487 cc(k, j, i) = (c4*cc(k, j - 1, i) - i*cc(k, j - 1, i - 1))/c3
490 cc(j, j, 0) = -c1*c5*cc(j - 1, j - 1, 0) / c2
494 cc(j, j, i) = c1*(i*cc(j - 1, j - 1, i - 1) - c5*cc(j - 1, j - 1, i))/c2
505 real(real64)
pure function
ddelta(i, j)
506 integer,
intent(in) :: i
507 integer,
intent(in) :: j
524 integer,
intent(in) :: n
525 integer,
allocatable,
intent(out) :: out(:)
526 integer,
intent(out) :: length
527 integer,
optional,
intent(in) :: in(:)
533 if (
present(in))
then
534 length = ubound(in, 1)
535 safe_allocate(out(1:length))
539 safe_allocate(out(1:length))
552 logical function member(n, a)
553 integer,
intent(in) :: n
554 integer,
intent(in) :: a(:)
562 do i = 1, ubound(a, 1)
574 integer,
intent(in) :: nn
575 real(real64),
intent(in) :: xa(:)
576 real(real64),
intent(in) :: xx
577 real(real64),
intent(out) :: cc(:)
587 cc(ii) = cc(ii)*(xx - xa(kk))/(xa(ii) - xa(kk))
596 logical pure function
even(n)
597 integer,
intent(in) :: n
599 even = (mod(n, 2) == 0)
606 logical pure function
odd(n)
607 integer,
intent(in) :: n
617 real(real64),
intent(in) :: x(:)
618 real(real64),
intent(out) :: u(:)
621 real(real64) :: sumx2
622 real(real64),
allocatable :: xx(:)
628 assert(
size(u) == n-1)
631 safe_allocate(xx(1:n))
633 if (abs(x(j))<1.0e-8_real64)
then
644 sumx2 = sumx2 + xx(j)**2
646 if (abs(sumx2) <= m_epsilon .and. abs(xx(k)) <= m_epsilon)
exit
650 u(n-1) =
atan2(xx(n), xx(n-1))
664 real(real64),
intent(in) :: u(:)
665 real(real64),
intent(out) :: x(:)
673 assert(
size(u) == n-1)
680 x(2) =
sin(u(1))*
cos(u(2))
681 x(3) =
sin(u(1))*
sin(u(2))
684 x(2) =
sin(u(1))*
cos(u(2))
691 x(j) = x(j) *
cos(u(j))
695 x(n) = x(n) *
sin(u(k))
697 x(n) = x(n) *
sin(u(n-1))
710 real(real64),
intent(out) :: grad_matrix(:,:)
711 real(real64),
intent(in) :: r
712 real(real64),
intent(in) :: x(:)
722 grad_matrix(1,1) = -r*
sin(x(1))
724 grad_matrix(m,1) = m_zero
731 grad_matrix(m,l) = -r*grad_matrix(m,l)*product(
sin(x(1:l)))
733 grad_matrix(m,l) = grad_matrix(m,l)*r*
cos(x(m))*
cos(x(l))*product(
sin(x(1:m-1)))*product(
sin(x(m+1:l-1)))
735 grad_matrix(m,l) = m_zero
742 grad_matrix(m,n) = r*
cos(x(m))*grad_matrix(m,n)*product(
sin(x(1:m-1)))*product(
sin(x(m+1:n-1)))
749 integer(int64) pure function
pad88(size, blk)
750 integer(int64),
intent(in) :: size
751 integer(int64),
intent(in) :: blk
759 pad88 =
size + blk - mm
764 integer(int64) pure function
pad48(size, blk)
765 integer,
intent(in) :: size
766 integer(int64),
intent(in) :: blk
772 integer(int64) pure function
pad8(size, blk)
773 integer(int64),
intent(in) :: size
774 integer,
intent(in) :: blk
780 integer pure function
pad4(size, blk)
781 integer,
intent(in) :: size
782 integer,
intent(in) :: blk
784 pad4 = int(
pad88(int(
size, int64), int(blk, int64)), int32)
793 integer pure function
pad_pow2(size)
794 integer,
intent(in) :: size
802 if (
size == 0)
return
808 if (mm - 2*mm2 /= 0)
exit
822 real(real64)
pure function
dlog2(xx)
823 real(real64),
intent(in) :: xx
830 integer pure function
ilog2(xx)
831 integer,
intent(in) :: xx
838 integer(int64) pure function
llog2(xx)
839 integer(int64),
intent(in) :: xx
841 llog2 = nint(
log2(real(xx, real64) ), kind=int64)
847 complex(real64),
intent(in) :: z
858 complex(real64) pure function
phi1(z)
859 complex(real64),
intent(in) :: z
861 real(real64),
parameter :: cut = 0.002_real64
863 if (abs(z) < cut)
then
864 phi1 = m_one + m_half*z*(m_one + z/m_three*(m_one + m_fourth*z*(m_one + z/m_five)))
877 complex(real64) pure function
phi2(z)
878 complex(real64),
intent(in) :: z
880 real(real64),
parameter :: cut = 0.002_real64
882 if (abs(z) < cut)
then
883 phi2 = m_half + z/6.0_real64*(m_one + m_fourth*z*(m_one + z/m_five*(m_one + z/6.0_real64)))
885 phi2 = (
exp(z) - z - m_z1)/z**2
892 real(real64),
intent(in) :: x
900 integer,
intent(in) :: n
907 message(1) =
"Internal error: is_prime does not take negative numbers."
908 call messages_fatal(1)
916 root = nint(
sqrt(real(n, real64) ))
918 if (mod(n,i) == 0)
then
935 real(real64),
intent(out) :: r(:,:)
936 real(real64),
intent(in) :: ff(:)
937 real(real64),
intent(in) :: tt(:)
940 real(real64) :: th, uv, uu, vv, ft
941 real(real64),
allocatable :: axis(:), u(:), v(:), f(:), t(:), p(:)
947 assert((dim < 3) .or. (dim > 2))
948 assert(
size(tt,1) == dim)
949 assert((
size(r,1) == dim) .and. (
size(r,2) == dim))
951 safe_allocate(u(1:dim))
952 safe_allocate(v(1:dim))
953 safe_allocate(f(1:dim))
954 safe_allocate(t(1:dim))
961 ft = dot_product(f,t)
963 if (abs(ft) < m_one)
then
976 safe_allocate(axis(1:dim))
979 u = f / dot_product(f,f)
980 v = t /dot_product(t,t)
983 axis = axis / norm2(axis)
985 r(1,1) =
cos(th) + axis(1)**2 * (1 -
cos(th))
986 r(1,2) = axis(1)*axis(2)*(1-
cos(th)) + axis(3)*
sin(th)
987 r(1,3) = axis(1)*axis(3)*(1-
cos(th)) - axis(2)*
sin(th)
989 r(2,1) = axis(2)*axis(1)*(1-
cos(th)) - axis(3)*
sin(th)
990 r(2,2) =
cos(th) + axis(2)**2 * (1 -
cos(th))
991 r(2,3) = axis(2)*axis(3)*(1-
cos(th)) + axis(1)*
sin(th)
993 r(3,1) = axis(3)*axis(1)*(1-
cos(th)) + axis(2)*
sin(th)
994 r(3,2) = axis(3)*axis(2)*(1-
cos(th)) - axis(1)*
sin(th)
995 r(3,3) =
cos(th) + axis(3)**2 * (1 -
cos(th))
997 safe_deallocate_a(axis)
1006 r(1,1) = u(1)**2 + (1-u(1)**2)*
cos(th)
1007 r(1,2) = u(1)*u(2)*(1-
cos(th)) - u(3)*
sin(th)
1008 r(1,3) = u(1)*u(3) + u(2)*
sin(th)
1010 r(2,1) = u(1)*u(2)*(1-
cos(th)) + u(3)*
sin(th)
1011 r(2,2) = u(2)**2 + (1-u(2)**2)*
cos(th)
1012 r(2,3) = u(2)*u(3)*(1-
cos(th)) - u(1)*
sin(th)
1014 r(3,1) = u(1)*u(3)*(1-
cos(th)) - u(2)*
sin(th)
1015 r(3,2) = u(2)*u(3)*(1-
cos(th)) + u(1)*
sin(th)
1016 r(3,3) = u(3)**2 + (1-u(3)**2)*
cos(th)
1021 safe_allocate(p(1:dim))
1023 if (abs(f(1)) <= abs(f(2)) .and. abs(f(1)) < abs(f(3)))
then
1024 p = (/m_one, m_zero, m_zero/)
1025 else if (abs(f(2)) < abs(f(1)) .and. abs(f(2)) <= abs(f(3)))
then
1026 p = (/m_zero, m_one, m_zero/)
1027 else if (abs(f(3)) <= abs(f(1)) .and. abs(f(3)) < abs(f(2)))
then
1028 p = (/m_zero, m_zero, m_one/)
1034 uu = dot_product(u,u)
1035 vv = dot_product(v,v)
1036 uv = dot_product(u,v)
1041 r(i,j) =
ddelta(i,j) - m_two * u(i)*u(j)/uu - m_two * v(i)*v(j)/vv &
1042 + 4_real64*uv * v(i)*u(j) /(uu*vv)
1047 safe_deallocate_a(p)
1063 safe_deallocate_a(u)
1064 safe_deallocate_a(v)
1080 real(real64),
intent(in) :: x, h
1081 real(real64),
intent(out) :: res, err
1084 use,
intrinsic :: iso_fortran_env
1086 real(real64),
intent(in) :: x
1087 real(real64),
intent(inout) :: fx
1091 real(real64) :: con = 1.4_real64, big = 1.0e30_real64, safe = m_two
1092 integer :: ntab = 20, i, j
1093 real(real64) :: errt, fac, hh, fx1, fx2
1094 real(real64),
allocatable :: a(:, :)
1098 if (abs(h) <= m_epsilon)
then
1099 message(1) =
"h must be nonzero in numder_ridders"
1100 call messages_fatal(1)
1103 safe_allocate(a(1:ntab, 1:ntab))
1108 a(1,1) = (fx1-fx2) / (m_two*hh)
1114 a(1,i) = (fx1-fx2) / (m_two*hh)
1117 a(j,i) = (a(j-1,i)*fac-a(j-1,i-1)) / (fac-m_one)
1119 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
1120 if (errt .le. err)
then
1125 if (abs(a(i,i)-a(i-1,i-1)) .ge. safe*err)
return
1128 safe_deallocate_a(a)
1134 real(real64),
intent(in) :: a(:)
1135 complex(real64),
intent(in) :: b(:)
1137 complex(real64) :: c(1:3)
1139 c(1) = a(2)*b(3) - a(3)*b(2)
1140 c(2) = a(3)*b(1) - a(1)*b(3)
1141 c(3) = a(1)*b(2) - a(2)*b(1)
1147 complex(real64),
intent(in) :: a(:)
1148 real(real64),
intent(in) :: b(:)
1150 complex(real64) :: c(1:3)
1152 c(1) = a(2)*b(3) - a(3)*b(2)
1153 c(2) = a(3)*b(1) - a(1)*b(3)
1154 c(3) = a(1)*b(2) - a(2)*b(1)
1260 integer,
intent(in) :: np
1261 integer,
intent(in) :: nn, mm
1262 real(real64),
intent(in) :: xx(np)
1263 real(real64),
intent(out) :: cx(np)
1266 real(real64),
allocatable :: cx_tmp(:,:)
1270 safe_allocate(cx_tmp(1:np, 0:nn))
1282 cx_tmp(1:np,0) = m_one
1283 cx_tmp(1:np,1) = real(mm + 1, real64) - xx(1:np)
1286 cx_tmp(1:np, ii) = &
1287 (( real(mm + 2*ii - 1, real64) - xx(1:np)) * cx_tmp(1:np,ii-1) &
1288 + real(-mm -ii + 1, real64) * cx_tmp(1:np,ii-2)) / real(ii, real64)
1291 cx(1:np) = cx_tmp(1:np, nn)
1293 safe_deallocate_a(cx_tmp)
1298 integer,
intent(in) :: nn
1299 real(real64),
intent(inout) :: aa(:, :)
1305 aa(jj, ii) = aa(ii, jj)
1311 integer,
intent(in) :: nn
1312 complex(real64),
intent(inout) :: aa(:, :)
1318 aa(jj, ii) = conjg(aa(ii, jj))
1320 aa(jj, jj) = real(aa(jj, jj), real64)
1325 integer,
intent(in) :: nn
1326 real(real64),
intent(inout) :: aa(:, :)
1332 aa(ii, jj) = aa(jj, ii)
1338 integer,
intent(in) :: nn
1339 complex(real64),
intent(inout) :: aa(:, :)
1345 aa(ii, jj) = conjg(aa(jj, ii))
1347 aa(jj, jj) = real(aa(jj, jj), real64)
1353 integer,
intent(in) :: nn
1354 real(real64),
intent(inout) :: aa(:, :)
1360 aa(ii, jj) = m_half * (aa(ii, jj) + aa(jj, ii))
1361 aa(jj, ii) = aa(ii, jj)
1367 integer,
intent(in) :: nn
1368 real(real64),
intent(inout) :: aa(:, :)
1378 logical pure function
is_diagonal(dim, matrix, tol)
1379 integer,
intent(in) :: dim
1380 real(real64),
intent(in) :: matrix(:, :)
1381 real(real64),
intent(in) :: tol
1384 real(real64) :: max_diag
1388 max_diag = max(abs(matrix(ii, ii)), max_diag)
1395 if (abs(matrix(ii, jj)) > tol*max_diag)
is_diagonal = .false.
1402#include "complex.F90"
1403#include "math_inc.F90"
1407#include "math_inc.F90"
This is the common interface to a simple-minded polynomical interpolation procedure (simple use of th...
double acos(double __x) __attribute__((__nothrow__
double log(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double atan2(double __y, double __x) __attribute__((__nothrow__
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_pi
some mathematical constants
complex(real64), parameter, public m_z0
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_three
This module is intended to contain "only mathematical" functions and procedures.
subroutine, public dsymmetrize_matrix(nn, aa)
complex(real64) pure function, public phi2(z)
Compute phi2(z) = (phi1(z)-1)/z = (exp(z) - z - 1)/z^2.
integer pure function ilog2(xx)
logical pure function, public even(n)
Returns if n is even.
pure real(real64) function, dimension(dim, dim) ddiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the integer versi...
subroutine, public generalized_laguerre_polynomial(np, nn, mm, xx, cx)
complex(real64) pure function, public exponential(z)
Wrapper for exponential.
subroutine dinterpolate_0(xa, ya, x, y)
subroutine, public ylmr_cmplx(xx, li, mi, ylm)
Computes spherical harmonics ylm at position (x, y, z)
subroutine zinterpolate_0(xa, ya, x, y)
subroutine dinterpolate_1(xa, ya, x, y)
real(real64) pure function dlog2(xx)
integer(int64) pure function pad88(size, blk)
real(real64) pure function, public square_root(x)
Wrapper for sqrt.
logical function, public member(n, a)
Considers a(1:ubound(a, 1)) as an integer set and checks if n is a member of it.
integer pure function pad4(size, blk)
subroutine dlower_triangular_to_hermitian(nn, aa)
pure complex(real64) function, dimension(1:3), public dzcross_product(a, b)
subroutine, public cartesian2hyperspherical(x, u)
Performs a transformation of an n-dimensional vector from Cartesian coordinates to hyperspherical coo...
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
subroutine zinterpolate_1(xa, ya, x, y)
pure real(real64) function, dimension(1:3), public dcross_product(a, b)
integer(int64) pure function pad48(size, blk)
subroutine, public ylmr_real(xx, li, mi, ylm)
This is a Numerical Recipes-based subroutine computes real spherical harmonics ylm at position (x,...
subroutine dupper_triangular_to_hermitian(nn, aa)
subroutine, public interpolation_coefficients(nn, xa, xx, cc)
subroutine, public generate_rotation_matrix(R, ff, tt)
Generates a rotation matrix R to rotate a vector f to t.
elemental logical function zis_close_scalar(x, y, rtol, atol)
Same as dis_close_scalar for complex numbers.
logical pure function, public odd(n)
Returns if n is odd.
real(real64) pure function, public ddelta(i, j)
subroutine, public make_idx_set(n, out, length, in)
Construct out(1:length) = (/1, ..., n/) if in is not present, out(1:length) = in otherwise.
logical function, public is_prime(n)
complex(real64) pure function, public phi1(z)
Compute phi1(z) = (exp(z)-1)/z.
pure complex(real64) function, dimension(1:3), public zcross_product(a, b)
subroutine, public weights(N, M, cc, side)
Compute the weights for finite-difference calculations:
subroutine zinterpolate_2(xa, ya, x, y)
pure complex(real64) function, dimension(dim, dim) zdiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the integer versi...
subroutine zlower_triangular_to_hermitian(nn, aa)
logical pure function, public is_diagonal(dim, matrix, tol)
Returns true is the matrix of size dim x dim is diagonal, given a relative tolerance.
subroutine zupper_triangular_to_hermitian(nn, aa)
pure integer function, dimension(dim, dim) idiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the real and comp...
integer(int64) pure function pad8(size, blk)
subroutine dinterpolate_2(xa, ya, x, y)
subroutine, public numder_ridders(x, h, res, err, f)
Numerical derivative (Ridder`s algorithm).
integer(int64) pure function llog2(xx)
recursive real(real64) function, public hermite(n, x)
elemental logical function dis_close_scalar(x, y, rtol, atol)
Are and equal within a tolerance.
pure complex(real64) function, dimension(1:3), public zdcross_product(a, b)
subroutine, public hyperspherical2cartesian(u, x)
Performs the inverse transformation of cartesian2hyperspherical.
subroutine, public dzero_small_elements_matrix(nn, aa, tol)
subroutine, public hypersphere_grad_matrix(grad_matrix, r, x)
Gives the hyperspherical gradient matrix, which contains the derivatives of the Cartesian coordinates...
recursive integer function, public factorial(n)