VSH_POL_INVERSION Subroutine

public subroutine VSH_POL_INVERSION(LMAX, NTH, OUTUNIT, STATUS)

Verifies the orthogonal-rotation identity linking the standard and polar VSH bases can be inverted exactly: i.e. rotating VSH_POL_UP/VSH_POL_DN back by the same angle used to build them from PVSH_POL/PVSH_RAD recovers the originals exactly.

theta phi max_absdiff_pol max_absdiff_rad`, each the max over the 3 vector components). ).

Arguments

Type IntentOptional Attributes Name
integer(kind=i4), intent(in) :: LMAX

Maximum degree tested, .
integer(kind=i4), intent(in) :: NTH

Number of colatitude points swept over .
integer(kind=i4), intent(in) :: OUTUNIT

Fortran unit number to write results to (columns: `L M
integer(kind=i4), intent(out) :: STATUS

Output: 0 = pass, 1 = fail (either max difference

Calls

proc~~vsh_pol_inversion~~CallsGraph proc~vsh_pol_inversion VSH_POL_INVERSION interface~pvsh_pol_all PVSH_POL_ALL proc~vsh_pol_inversion->interface~pvsh_pol_all interface~pvsh_rad_all PVSH_RAD_ALL proc~vsh_pol_inversion->interface~pvsh_rad_all interface~vsh_pol_dn_all VSH_POL_DN_ALL proc~vsh_pol_inversion->interface~vsh_pol_dn_all interface~vsh_pol_up_all VSH_POL_UP_ALL proc~vsh_pol_inversion->interface~vsh_pol_up_all proc~ylm_index YLM_INDEX proc~vsh_pol_inversion->proc~ylm_index 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~vsh_pol_dn_all_dp VSH_POL_DN_ALL_DP interface~vsh_pol_dn_all->proc~vsh_pol_dn_all_dp proc~vsh_pol_up_all_dp VSH_POL_UP_ALL_DP interface~vsh_pol_up_all->proc~vsh_pol_up_all_dp proc~pvsh_pol_all_dp->proc~ylm_index interface~vsh_core VSH_CORE proc~pvsh_pol_all_dp->interface~vsh_core proc~pvsh_rad_all_dp->interface~vsh_core proc~vsh_pol_dn_all_dp->proc~ylm_index proc~vsh_pol_dn_all_dp->interface~vsh_core proc~vsh_pol_up_all_dp->proc~ylm_index proc~vsh_pol_up_all_dp->interface~vsh_core proc~vsh_core_dp VSH_CORE_DP interface~vsh_core->proc~vsh_core_dp proc~vsh_core_dp->proc~ylm_index dcos dcos proc~vsh_core_dp->dcos 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~~vsh_pol_inversion~~CalledByGraph proc~vsh_pol_inversion VSH_POL_INVERSION program~main MAIN program~main->proc~vsh_pol_inversion

Source Code

  SUBROUTINE VSH_POL_INVERSION(LMAX, NTH, OUTUNIT, STATUS)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(IN) :: LMAX
  INTEGER(KIND=i4), INTENT(IN) :: NTH
  INTEGER(KIND=i4), INTENT(IN) :: OUTUNIT
  INTEGER(KIND=i4), INTENT(OUT) :: STATUS
  INTEGER(KIND=i4) :: L, M, C, ITH, NYLM, IDX
  REAL(KIND=dp) :: THETA, PHI, DTH, SC_POL, SC_RAD
  REAL(KIND=dp) :: DIFF_POL, DIFF_RAD
  COMPLEX(KIND=dp), ALLOCATABLE :: UP(:,:), DN(:,:)
  COMPLEX(KIND=dp), ALLOCATABLE :: POL(:,:), RAD(:,:)
  NYLM = (LMAX+1)**2
  ALLOCATE(UP(3,NYLM), DN(3,NYLM), POL(3,NYLM), RAD(3,NYLM))
  STATUS = 0
  DTH = pi / (NTH + 1)
  PHI = pi / 4.d0
  WRITE(OUTUNIT,'(A)') &
      '# L  M  theta  phi  max_absdiff_pol  max_absdiff_rad'
  DO ITH = 1, NTH
    THETA = ITH * DTH
    CALL VSH_POL_UP_ALL(UP,  LMAX, THETA, PHI)
    CALL VSH_POL_DN_ALL(DN,  LMAX, THETA, PHI)
    CALL PVSH_POL_ALL(POL, LMAX, THETA, PHI)
    CALL PVSH_RAD_ALL(RAD, LMAX, THETA, PHI)
    DO L = 0, LMAX
      SC_POL = DSQRT(DBLE(L)/(2*L+1.d0))
      SC_RAD = DSQRT((L+1.d0)/(2*L+1.d0))
      DO M = -L, L
        IDX = YLM_INDEX(L, M)
        DIFF_POL = 0.d0
        DIFF_RAD = 0.d0
        DO C = 1, 3
          DIFF_POL = MAX(DIFF_POL, ABS( &
              SC_POL*UP(C,IDX) + SC_RAD*DN(C,IDX) - POL(C,IDX)))
          DIFF_RAD = MAX(DIFF_RAD, ABS( &
             -SC_RAD*UP(C,IDX) + SC_POL*DN(C,IDX) - RAD(C,IDX)))
        END DO
        WRITE(OUTUNIT, *) L, M, THETA, PHI, DIFF_POL, DIFF_RAD
        IF (DIFF_POL > 1.0E-12_dp .OR. DIFF_RAD > 1.0E-12_dp) STATUS = 1
      END DO
    END DO
  END DO
  DEALLOCATE(UP, DN, POL, RAD)
  END SUBROUTINE VSH_POL_INVERSION