RUN_GWI_SWEEP Subroutine

subroutine RUN_GWI_SWEEP(N_FAIL)

GWI vs. the brute-force Gaunt-coefficient integral (via SHGLQ theta-quadrature times a uniform phi sum), swept over JMAX, all valid . Writes examples/gwi_gwj_sweep/gwi_check.dat.

across the sweep exceeds TOL.

Arguments

Type IntentOptional Attributes Name
integer(kind=i4), intent(inout) :: N_FAIL

Running failure counter, incremented if the max error

Calls

proc~~run_gwi_sweep~~CallsGraph proc~run_gwi_sweep RUN_GWI_SWEEP dacos dacos proc~run_gwi_sweep->dacos interface~gwi GWI proc~run_gwi_sweep->interface~gwi interface~shglq SHGLQ proc~run_gwi_sweep->interface~shglq interface~ssh SSH proc~run_gwi_sweep->interface~ssh proc~gwi_dp GWI_DP interface~gwi->proc~gwi_dp proc~shglq_dp SHGLQ_DP interface~shglq->proc~shglq_dp proc~ssh_dp SSH_DP interface~ssh->proc~ssh_dp interface~cgcoeff CGCOEFF proc~gwi_dp->interface~cgcoeff dcos dcos proc~shglq_dp->dcos proc~ssh_dp->dcos interface~assoc_legendre ASSOC_LEGENDRE proc~ssh_dp->interface~assoc_legendre interface~log_fact LOG_FACT proc~ssh_dp->interface~log_fact proc~assoc_legendre_dp ASSOC_LEGENDRE_DP interface~assoc_legendre->proc~assoc_legendre_dp proc~cgcoeff_dp CGCOEFF_DP interface~cgcoeff->proc~cgcoeff_dp proc~log_fact_dp LOG_FACT_DP interface~log_fact->proc~log_fact_dp proc~assoc_legendre_dp->interface~log_fact proc~cgcoeff_dp->interface~log_fact

Called by

proc~~run_gwi_sweep~~CalledByGraph proc~run_gwi_sweep RUN_GWI_SWEEP program~gwi_gwj_sweep GWI_GWJ_SWEEP program~gwi_gwj_sweep->proc~run_gwi_sweep

Source Code

  SUBROUTINE RUN_GWI_SWEEP(N_FAIL)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(INOUT) :: N_FAIL
  REAL(KIND=dp) :: ZERO(LMAX_QUAD+1), WTHETA(LMAX_QUAD+1)
  REAL(KIND=dp) :: DPHI, THETA, PHI, MAX_ERR, ABSDIFF
  COMPLEX(KIND=dp) :: CLOSED, BRUTE
  INTEGER(KIND=i4) :: J1, J2, M1, M2, L, M, IQ, IP

  CALL SHGLQ(ZERO, WTHETA, LMAX_QUAD)
  DPHI = 2.d0*pi/NPHI
  MAX_ERR = 0.0_dp

  OPEN(UNIT=81, FILE="./examples/gwi_gwj_sweep/gwi_check.dat")
  WRITE(81,'(A)') '# J1  M1  J2  M2  L  M  '// &
                  'Re(closed)  Im(closed)  Re(brute)  Im(brute)  absdiff'

  DO J1 = 0, JMAX
    DO J2 = 0, JMAX
      DO M1 = -J1, J1
        DO M2 = -J2, J2
          DO L = ABS(J1-J2), J1+J2
            M = M1+M2
            IF (ABS(M) > L) CYCLE

            CLOSED = GWI(J1,M1,J2,M2,L,M)

            BRUTE = DCMPLX(0.d0, 0.d0)
            DO IQ = 1, LMAX_QUAD+1
              THETA = DACOS(ZERO(IQ))
              DO IP = 0, NPHI-1
                PHI = IP*DPHI
                BRUTE = BRUTE + WTHETA(IQ)*DPHI * SSH(J1,M1,THETA,PHI) * &
                        SSH(J2,M2,THETA,PHI) * CONJG(SSH(L,M,THETA,PHI))
              ENDDO
            ENDDO

            ABSDIFF = ABS(CLOSED-BRUTE)
            MAX_ERR = MAX(MAX_ERR, ABSDIFF)
            WRITE(81,*) J1, M1, J2, M2, L, M, DBLE(CLOSED), AIMAG(CLOSED), &
                        DBLE(BRUTE), AIMAG(BRUTE), ABSDIFF
          ENDDO
        ENDDO
      ENDDO
    ENDDO
  ENDDO
  CLOSE(81)

  IF (MAX_ERR > TOL) THEN
    N_FAIL = N_FAIL + 1
    WRITE(*,'(A,ES10.3)') "FAIL  gwi_gaunt_sweep  max_err=", MAX_ERR
  ELSE
    WRITE(*,'(A,ES10.3)') "PASS  gwi_gaunt_sweep  max_err=", MAX_ERR
  END IF
  END SUBROUTINE RUN_GWI_SWEEP