DECOMPOSE Subroutine

subroutine DECOMPOSE(FIELD_FN, LMAX, LMAX_QUAD, NPHI, A_RAD, A_POL, A_TOR)

Decomposes FIELD_FN into for , via 2D Gauss-Legendre() times uniform-trapezoidal() quadrature. Identical to vsh_decomposition.f90's DECOMPOSE -- duplicated rather than shared, matching this repository's convention of fully independent worked-example programs (see CMakeLists.txt).

nodes); must scale with LMAX for an accurate result. YLM_INDEX(l,m).

Arguments

Type IntentOptional Attributes Name
function FIELD_FN(THETA, PHI) result(F)
Arguments
Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: THETA
real(kind=dp), intent(in) :: PHI
Return Value complex(kind=dp), DIMENSION(3)
integer(kind=i4), intent(in) :: LMAX

Maximum degree to decompose, .
integer(kind=i4), intent(in) :: LMAX_QUAD

Theta-quadrature order (SHGLQ uses LMAX_QUAD+1
integer(kind=i4), intent(in) :: NPHI

Number of uniform phi-quadrature points.
complex(kind=dp), intent(out) :: A_RAD((LMAX+1)**2)
complex(kind=dp), intent(out) :: A_POL((LMAX+1)**2)
complex(kind=dp), intent(out) :: A_TOR((LMAX+1)**2)

Calls

proc~~decompose~2~~CallsGraph proc~decompose~2 DECOMPOSE dacos dacos proc~decompose~2->dacos interface~dot DOT proc~decompose~2->interface~dot interface~pvsh_pol_all PVSH_POL_ALL proc~decompose~2->interface~pvsh_pol_all interface~pvsh_rad_all PVSH_RAD_ALL proc~decompose~2->interface~pvsh_rad_all interface~pvsh_tor_all PVSH_TOR_ALL proc~decompose~2->interface~pvsh_tor_all interface~shglq SHGLQ proc~decompose~2->interface~shglq proc~dot_dp DOT_DP interface~dot->proc~dot_dp proc~pvsh_pol_all_dp PVSH_POL_ALL_DP interface~pvsh_pol_all->proc~pvsh_pol_all_dp proc~pvsh_rad_all_dp PVSH_RAD_ALL_DP interface~pvsh_rad_all->proc~pvsh_rad_all_dp proc~pvsh_tor_all_dp PVSH_TOR_ALL_DP interface~pvsh_tor_all->proc~pvsh_tor_all_dp proc~shglq_dp SHGLQ_DP interface~shglq->proc~shglq_dp interface~vsh_core VSH_CORE proc~pvsh_pol_all_dp->interface~vsh_core proc~ylm_index YLM_INDEX proc~pvsh_pol_all_dp->proc~ylm_index proc~pvsh_rad_all_dp->interface~vsh_core proc~pvsh_tor_all_dp->interface~vsh_core proc~pvsh_tor_all_dp->proc~ylm_index dcos dcos proc~shglq_dp->dcos proc~vsh_core_dp VSH_CORE_DP interface~vsh_core->proc~vsh_core_dp proc~vsh_core_dp->dcos proc~vsh_core_dp->proc~ylm_index dsin dsin proc~vsh_core_dp->dsin interface~assoc_legendre_norm_all ASSOC_LEGENDRE_NORM_ALL proc~vsh_core_dp->interface~assoc_legendre_norm_all interface~ddx_assoc_legendre_norm_all DDX_ASSOC_LEGENDRE_NORM_ALL proc~vsh_core_dp->interface~ddx_assoc_legendre_norm_all proc~plm_index PLM_INDEX proc~vsh_core_dp->proc~plm_index proc~assoc_legendre_norm_all_dp ASSOC_LEGENDRE_NORM_ALL_DP interface~assoc_legendre_norm_all->proc~assoc_legendre_norm_all_dp proc~ddx_assoc_legendre_norm_all_dp DDX_ASSOC_LEGENDRE_NORM_ALL_DP interface~ddx_assoc_legendre_norm_all->proc~ddx_assoc_legendre_norm_all_dp proc~assoc_legendre_norm_all_dp->proc~plm_index proc~ddx_assoc_legendre_norm_all_dp->proc~plm_index

Called by

proc~~decompose~2~~CalledByGraph proc~decompose~2 DECOMPOSE proc~run_part_a~2 RUN_PART_A proc~run_part_a~2->proc~decompose~2 proc~run_part_b~2 RUN_PART_B proc~run_part_b~2->proc~decompose~2 program~vsh_decomposition_tor VSH_DECOMPOSITION_TOR program~vsh_decomposition_tor->proc~run_part_a~2 program~vsh_decomposition_tor->proc~run_part_b~2

Source Code

  SUBROUTINE DECOMPOSE(FIELD_FN, LMAX, LMAX_QUAD, NPHI, A_RAD, A_POL, A_TOR)
  IMPLICIT NONE
  INTERFACE
    FUNCTION FIELD_FN(THETA, PHI) RESULT(F)
      IMPORT :: dp
      REAL(KIND=dp), INTENT(IN) :: THETA, PHI
      COMPLEX(KIND=dp), DIMENSION(3) :: F
    END FUNCTION FIELD_FN
  END INTERFACE
  INTEGER(KIND=i4), INTENT(IN)  :: LMAX
  INTEGER(KIND=i4), INTENT(IN)  :: LMAX_QUAD
  INTEGER(KIND=i4), INTENT(IN)  :: NPHI
  COMPLEX(KIND=dp), INTENT(OUT) :: A_RAD((LMAX+1)**2)
  COMPLEX(KIND=dp), INTENT(OUT) :: A_POL((LMAX+1)**2)
  COMPLEX(KIND=dp), INTENT(OUT) :: A_TOR((LMAX+1)**2)
  REAL(KIND=dp), ALLOCATABLE :: ZERO(:), WTHETA(:)
  COMPLEX(KIND=dp), ALLOCATABLE :: RADALL(:,:), POLALL(:,:), TORALL(:,:)
  COMPLEX(KIND=dp), DIMENSION(3) :: FVEC
  REAL(KIND=dp) :: THETA, PHI, DPHI, WEIGHT
  INTEGER(KIND=i4) :: NYLM, IQ, IP, IDX

  NYLM = (LMAX+1)**2
  ALLOCATE(ZERO(LMAX_QUAD+1), WTHETA(LMAX_QUAD+1))
  ALLOCATE(RADALL(3,NYLM), POLALL(3,NYLM), TORALL(3,NYLM))
  CALL SHGLQ(ZERO, WTHETA, LMAX_QUAD)

  DPHI = 2.d0*pi/NPHI
  A_RAD = DCMPLX(0.d0, 0.d0)
  A_POL = DCMPLX(0.d0, 0.d0)
  A_TOR = DCMPLX(0.d0, 0.d0)

  DO IQ = 1, LMAX_QUAD+1
    THETA = DACOS(ZERO(IQ))
    DO IP = 0, NPHI-1
      PHI = IP*DPHI
      FVEC = FIELD_FN(THETA, PHI)
      CALL PVSH_RAD_ALL(RADALL, LMAX, THETA, PHI)
      CALL PVSH_POL_ALL(POLALL, LMAX, THETA, PHI)
      CALL PVSH_TOR_ALL(TORALL, LMAX, THETA, PHI)
      WEIGHT = WTHETA(IQ) * DPHI
      DO IDX = 1, NYLM
        A_RAD(IDX) = A_RAD(IDX) + WEIGHT * DOT(FVEC, CONJG(RADALL(:,IDX)))
        A_POL(IDX) = A_POL(IDX) + WEIGHT * DOT(FVEC, CONJG(POLALL(:,IDX)))
        A_TOR(IDX) = A_TOR(IDX) + WEIGHT * DOT(FVEC, CONJG(TORALL(:,IDX)))
      ENDDO
    ENDDO
  ENDDO

  DEALLOCATE(ZERO, WTHETA, RADALL, POLALL, TORALL)
  END SUBROUTINE DECOMPOSE