Accumulate potential, induced by multipole spherical harmonics
This function relies on a user-provided temporary workspace
Computes the following sum: \f[ v = \beta v + \alpha \sum_{\ell=0}^p \frac{4\pi}{\sqrt{2\ell+1}} \left( \frac{r}{|c|} \right)^{\ell+1} \sum_{m=-\ell}^\ell M_\ell^m Y_\ell^m \left( \frac{c}{|c|} \right), \f] where \f$ M \f$ is a vector of coefficients of input harmonics of a degree up to \f$ p \f$ inclusively with a convergence radius \f$ r \f$ located at the origin, \f$ \alpha \f$ and \f$ \beta \f$ are scaling factors and \f$ c \f$ is a location of a particle.
Based on normalized real spherical harmonics \f$ Y_\ell^m \f$, scaled by \f$ r^{-\ell} \f$. It means corresponding coefficients are simply scaled by an additional factor \f$ r^\ell \f$.
@param[in] c: Coordinates of a particle (relative to center of harmonics)
@param[in] r: Radius of spherical harmonics
@param[in] p: Maximal degree of multipole basis functions
@param[in] vscales_rel: Relative normalization constants for
\f$ Y_\ell^m \f$. Dimension is (p+1)**2
@param[in] alpha: Scalar multiplier for src_m
@param[in] src_m: Multipole coefficients. Dimension is (p+1)**2
@param[in] beta: Scalar multiplier for v
@param[inout] v: Value of induced potential
@param[out] work: Temporary workspace of size (p+1)
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=rp), | intent(in) | :: | c(3) | |||
| real(kind=rp), | intent(in) | :: | src_r | |||
| integer, | intent(in) | :: | p | |||
| real(kind=rp), | intent(in) | :: | vscales_rel((p+1)*(p+1)) | |||
| real(kind=rp), | intent(in) | :: | alpha | |||
| real(kind=rp), | intent(in) | :: | src_m((p+1)*(p+1)) | |||
| real(kind=rp), | intent(in) | :: | beta | |||
| real(kind=rp), | intent(inout) | :: | dst_v | |||
| real(kind=rp), | intent(out), | target | :: | work(p+1) |
subroutine fmm_m2p_work(c, src_r, p, vscales_rel, alpha, src_m, &
& beta, dst_v, work)
! Inputs
integer, intent(in) :: p
real(rp), intent(in) :: c(3), src_r, vscales_rel((p+1)*(p+1)), alpha, &
& src_m((p+1)*(p+1)), beta
! Output
real(rp), intent(inout) :: dst_v
! Workspace
real(rp), intent(out), target :: work(p+1)
! Local variables
real(rp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
& tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
integer :: l, m, indl
! Scale output
if (abs(beta) < eps_rp) then
dst_v = 0.0_rp
else
dst_v = beta * dst_v
end if
! In case of 0.0_rp alpha nothing else is required no matter what is the
! value of the induced potential
if (abs(alpha) < eps_rp) then
return
end if
! Get spherical coordinates
if (abs(c(1)) < eps_rp) then
max12 = abs(c(2))
ssq12 = 1.0_rp
else if (abs(c(2)) .gt. abs(c(1))) then
max12 = abs(c(2))
ssq12 = 1.0_rp + (c(1)/c(2))**2
else
max12 = abs(c(1))
ssq12 = 1.0_rp + (c(2)/c(1))**2
end if
! Then we compute rho
if (abs(c(3)) < eps_rp) then
rho = max12 * sqrt(ssq12)
else if (abs(c(3)) .gt. max12) then
rho = 1.0_rp + ssq12 *(max12/c(3))**2
rho = abs(c(3)) * sqrt(rho)
else
rho = ssq12 + (c(3)/max12)**2
rho = max12 * sqrt(rho)
end if
! In case of a singularity (rho=0.0_rp) induced potential is infinite and is
! not taken into account.
if (abs(rho) < eps_rp) then
return
end if
! Compute the actual induced potential
rcoef = src_r / rho
! Length of a vector x(1:2)
stheta = max12 * sqrt(ssq12)
! Case x(1:2) != 0
if (abs(stheta) > eps_rp) then
! Normalize cphi and sphi
cphi = c(1) / stheta
sphi = c(2) / stheta
! Normalize ctheta and stheta
ctheta = c(3) / rho
stheta = stheta / rho
! Treat easy cases
select case(p)
case (0)
! l = 0
dst_v = dst_v + alpha*rcoef*vscales_rel(1)*src_m(1)
return
case (1)
! l = 0
tmp = src_m(1) * vscales_rel(1)
! l = 1
tmp2 = ctheta * src_m(3) * vscales_rel(3)
tmp3 = vscales_rel(4) * stheta
tmp2 = tmp2 - tmp3*sphi*src_m(2)
tmp2 = tmp2 - tmp3*cphi*src_m(4)
dst_v = dst_v + alpha*rcoef*(tmp+rcoef*tmp2)
return
end select
! Now p>1
! Precompute alpha*rcoef^{l+1}
work(1) = alpha * rcoef
do l = 1, p
work(l+1) = rcoef * work(l)
end do
! Case m = 0
! P_{l-2}^m which is P_0^0 now
pl2m = 1.0_rp
dst_v = dst_v + work(1)*src_m(1)*vscales_rel(1)
! Update P_m^m for the next iteration
pmm = -stheta
! P_{l-1}^m which is P_{m+1}^m now
pl1m = ctheta
ylm = pl1m * vscales_rel(3)
dst_v = dst_v + work(2)*src_m(3)*ylm
! P_l^m for l>m+1
do l = 2, p
plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
plm = plm / dble(l)
ylm = plm * vscales_rel(l*l+l+1)
dst_v = dst_v + work(l+1)*src_m(l*l+l+1)*ylm
! Update P_{l-2}^m and P_{l-1}^m for the next iteration
pl2m = pl1m
pl1m = plm
end do
! Prepare cos(m*phi) and sin(m*phi) for m=1
cmphi = cphi
smphi = sphi
! Case 0<m<p
do m = 1, p-1
! P_{l-2}^m which is P_m^m now
pl2m = pmm
ylm = pmm * vscales_rel((m+1)**2)
tmp1 = cmphi*src_m((m+1)**2) + smphi*src_m(m*m+1)
dst_v = dst_v + work(m+1)*ylm*tmp1
! Temporary to reduce number of operations
tmp1 = dble(2*m+1) * pmm
! Update P_m^m for the next iteration
pmm = -stheta * tmp1
! P_{l-1}^m which is P_{m+1}^m now
pl1m = ctheta * tmp1
ylm = pl1m * vscales_rel((m+1)*(m+3))
tmp1 = cmphi*src_m((m+1)*(m+3)) + smphi*src_m((m+1)*(m+2)+1-m)
dst_v = dst_v + work(m+2)*ylm*tmp1
! P_l^m for l>m+1
do l = m+2, p
plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
plm = plm / dble(l-m)
ylm = plm * vscales_rel(l*l+l+1+m)
tmp1 = cmphi*src_m(l*l+l+1+m) + smphi*src_m(l*l+l+1-m)
dst_v = dst_v + work(l+1)*ylm*tmp1
! Update P_{l-2}^m and P_{l-1}^m for the next iteration
pl2m = pl1m
pl1m = plm
end do
! Update cos(m*phi) and sin(m*phi) for the next iteration
tmp1 = cmphi
cmphi = cmphi*cphi - smphi*sphi
smphi = tmp1*sphi + smphi*cphi
end do
! Case m=p requires only to use P_m^m
ylm = pmm * vscales_rel((p+1)**2)
tmp1 = cmphi*src_m((p+1)**2) + smphi*src_m(p*p+1)
dst_v = dst_v + work(p+1)*ylm*tmp1
! Case of x(1:2) = 0 and x(3) != 0
else
! In this case Y_l^m = 0 for m != 0, so only case m = 0 is taken into
! account. Y_l^0 = ctheta^l in this case where ctheta is either +1 or
! -1. So, we copy sign(ctheta) into rcoef. But before that we
! initialize alpha/r factor for l=0 case without taking possible sign
! into account.
t = alpha * rcoef
if (c(3) .lt. 0.0_rp) then
rcoef = -rcoef
end if
! Proceed with accumulation of a potential
indl = 1
do l = 0, p
! Index of Y_l^0
indl = indl + 2*l
! Add 4*pi/(2*l+1)*rcoef^{l+1}*Y_l^0 contribution
dst_v = dst_v + t*src_m(indl)*vscales_rel(indl)
! Update t
t = t * rcoef
end do
end if
end subroutine fmm_m2p_work