Rotate spherical harmonics around OZ axis in an opposite direction
Compute the following matrix-vector product:
\f[
\mathrm{dst} = \beta \mathrm{dst} + \alpha R \mathrm{src},
\f]
where \f$ \mathrm{dst} \f$ is a vector of coefficients of output spherical
harmonics corresponding to a new cartesion system of coordinates, \f$
\mathrm{src} \f$ is a vector of coefficients of input spherical harmonics
corresponding to the standard cartesian system of coordinates and \f$
R \f$ is a matrix of rotation of coordinates around OZ axis on angle \f$
\phi \f$, presented by \f$ \cos(m \phi) \f$ and \f$ \sin(m \phi) \f$.
@param[in] p: Maximal order of spherical harmonics
@param[in] vcos: Vector \f$ { \cos(m \phi) }{m=0}^p \f$
@param[in] vsin: Vector \f$ { \sin(m \phi) }^p \f$
@param[in] alpha: Scalar multiplier for src
@param[in] src: Coefficients of initial spherical harmonics
@param[in] beta: Scalar multipler for dst
@param[inout] dst: Coefficients of rotated spherical harmonics
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed
arrows point from an interface to procedures which implement that interface.
This could include the module procedures in a generic interface or the
implementation in a submodule of an interface in a parent module.
subroutine fmm_sph_rotate_oz_adj_work(p,vcos,vsin,alpha,src,beta,dst)! Inputsinteger,intent(in)::preal(rp),intent(in)::vcos(p+1),vsin(p+1),alpha,src((p+1)*(p+1)),beta! Outputreal(rp),intent(inout)::dst((p+1)*(p+1))! Local variablesinteger::l,m,indreal(rp)::v1,v2,v3,v4! In case alpha is 0.0 just scale outputif(abs(alpha)<eps_rp)then! Set output to 0.0 if beta is also 0.0if(abs(beta)<eps_rp)thendst=0.0elsedst=beta*dstend if! Exit subroutinereturn end if! Now alpha is non-0.0! In case beta is 0.0 output is just overwritten without being readif(abs(beta)<eps_rp)then! l = 0dst(1)=alpha*src(1)! l > 0do l=1,pind=l*l+l+1! m = 0dst(ind)=alpha*src(ind)! m != 0do m=1,lv1=src(ind+m)v2=src(ind-m)v3=vcos(1+m)v4=vsin(1+m)! m > 0dst(ind+m)=alpha*(v1*v3+v2*v4)! m < 0dst(ind-m)=alpha*(v2*v3-v1*v4)end do end do else! l = 0dst(1)=beta*dst(1)+alpha*src(1)! l > 0do l=1,pind=l*l+l+1! m = 0dst(ind)=beta*dst(ind)+alpha*src(ind)! m != 0do m=1,lv1=src(ind+m)v2=src(ind-m)v3=vcos(1+m)v4=vsin(1+m)! m > 0dst(ind+m)=beta*dst(ind+m)+alpha*(v1*v3+v2*v4)! m < 0dst(ind-m)=beta*dst(ind-m)+alpha*(v2*v3-v1*v4)end do end do end if end subroutine fmm_sph_rotate_oz_adj_work