SSH_ORTHO Subroutine

public subroutine SSH_ORTHO(LMAX, NQUAD, OUTUNIT, STATUS)

Verifies scalar spherical harmonic orthonormality, , for every . The integral is done analytically ( for matching , by construction 0 otherwise, so only same- pairs are checked); the integral (in ) uses a plain NQUAD-point midpoint rule via ASSOC_LEGENDRE_ALL, not SHGLQ -- sufficient here since only a scalar pass/fail diagonal check is needed, not a high-order transform.

. L2 inner_product expected(0_or_1)`, upper triangle only). off-diagonal off from 0, by more than ).

Arguments

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

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

Number of midpoint-rule quadrature points over
integer(kind=i4), intent(in) :: OUTUNIT

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

Output: 0 = pass, 1 = fail (any diagonal off from 1, or

Calls

proc~~ssh_ortho~~CallsGraph proc~ssh_ortho SSH_ORTHO interface~assoc_legendre_all ASSOC_LEGENDRE_ALL proc~ssh_ortho->interface~assoc_legendre_all interface~log_fact LOG_FACT proc~ssh_ortho->interface~log_fact proc~plm_index PLM_INDEX proc~ssh_ortho->proc~plm_index proc~assoc_legendre_all_dp ASSOC_LEGENDRE_ALL_DP interface~assoc_legendre_all->proc~assoc_legendre_all_dp proc~log_fact_dp LOG_FACT_DP interface~log_fact->proc~log_fact_dp proc~assoc_legendre_all_dp->proc~plm_index

Called by

proc~~ssh_ortho~~CalledByGraph proc~ssh_ortho SSH_ORTHO program~main MAIN program~main->proc~ssh_ortho

Source Code

  SUBROUTINE SSH_ORTHO(LMAX, NQUAD, OUTUNIT, STATUS)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(IN) :: LMAX
  INTEGER(KIND=i4), INTENT(IN) :: NQUAD
  INTEGER(KIND=i4), INTENT(IN) :: OUTUNIT
  INTEGER(KIND=i4), INTENT(OUT) :: STATUS
  INTEGER(KIND=i4) :: L1, L2, M, IQ, PSIZE, EXPECTED
  REAL(KIND=dp) :: U, DU, N_L1M, N_L2M, FULL_INT
  REAL(KIND=dp), ALLOCATABLE :: P(:), ORTHO(:,:,:)

  PSIZE = (LMAX+1)*(LMAX+2)/2
  ALLOCATE(P(PSIZE))
  ALLOCATE(ORTHO(0:LMAX, 0:LMAX, 0:LMAX))
  STATUS = 0
  ORTHO = 0.0_dp
  DU = 2.0_dp / NQUAD

  !     Accumulate int_{-1}^{1} P_l1^m(u) * P_l2^m(u) du via midpoint rule
  DO IQ = 1, NQUAD
    U = -1.0_dp + (IQ - 0.5_dp) * DU
    CALL ASSOC_LEGENDRE_ALL(P, LMAX, U)
    DO M = 0, LMAX
      DO L1 = M, LMAX
        DO L2 = L1, LMAX
          ORTHO(L1, L2, M) = ORTHO(L1, L2, M) + &
              P(PLM_INDEX(L1, M)) * P(PLM_INDEX(L2, M))
        END DO
      END DO
    END DO
  END DO
  ORTHO = ORTHO * DU

  WRITE(OUTUNIT, '(A)') &
      '# L1  M  L2  inner_product  expected(0_or_1)'
  DO M = 0, LMAX
    DO L1 = M, LMAX
      DO L2 = L1, LMAX
        N_L1M = DSQRT((2.d0*L1+1.d0)/(4.d0*pi)) * &
            EXP(0.5_dp*(LOG_FACT(L1-M) - LOG_FACT(L1+M)))
        N_L2M = DSQRT((2.d0*L2+1.d0)/(4.d0*pi)) * &
            EXP(0.5_dp*(LOG_FACT(L2-M) - LOG_FACT(L2+M)))
        FULL_INT = ORTHO(L1, L2, M) * 2.d0*pi * N_L1M * N_L2M
        IF (L1 .EQ. L2) THEN
          EXPECTED = 1
        ELSE
          EXPECTED = 0
        END IF
        WRITE(OUTUNIT, *) L1, M, L2, FULL_INT, EXPECTED
        IF (L1 .EQ. L2) THEN
          IF (ABS(FULL_INT - 1.0_dp) > 1.0E-3_dp) STATUS = 1
        ELSE
          IF (ABS(FULL_INT) > 1.0E-3_dp) STATUS = 1
        END IF
      END DO
    END DO
  END DO

  DEALLOCATE(P, ORTHO)
  END SUBROUTINE SSH_ORTHO