SHGLQ_DP Subroutine

private subroutine SHGLQ_DP(ZERO, W, LMAX)

Gauss-Legendre quadrature nodes and weights, exact for integrating any polynomial of degree using only points -- the natural quadrature for spherical harmonic transforms, since is a polynomial in times , and the transform integral separates exactly into this -quadrature times a uniform sum over (itself exact for finite Fourier content). Nodes are the zeros of , found by Newton iteration (converges to machine precision in under 10 steps); weights are .

ascending from -1 to +1. \sum_i w_i f(x_i) ), size LMAX+1. uses LMAX+1 quadrature points.

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(out) :: ZERO(LMAX+1)
real(kind=dp), intent(out) :: W(LMAX+1)
integer(kind=i4), intent(in) :: LMAX

Requested exactness is degree ;

Calls

proc~~shglq_dp~~CallsGraph proc~shglq_dp SHGLQ_DP dcos dcos proc~shglq_dp->dcos

Called by

proc~~shglq_dp~~CalledByGraph proc~shglq_dp SHGLQ_DP interface~shglq SHGLQ interface~shglq->proc~shglq_dp proc~decompose DECOMPOSE proc~decompose->interface~shglq proc~decompose~2 DECOMPOSE proc~decompose~2->interface~shglq proc~run_gwi_sweep RUN_GWI_SWEEP proc~run_gwi_sweep->interface~shglq proc~run_part_a RUN_PART_A proc~run_part_a->proc~decompose proc~run_part_a~2 RUN_PART_A proc~run_part_a~2->proc~decompose~2 proc~run_part_b RUN_PART_B proc~run_part_b->proc~decompose proc~run_part_b~2 RUN_PART_B proc~run_part_b~2->proc~decompose~2 program~gwi_gwj_sweep GWI_GWJ_SWEEP program~gwi_gwj_sweep->proc~run_gwi_sweep program~vsh_decomposition VSH_DECOMPOSITION program~vsh_decomposition->proc~run_part_a program~vsh_decomposition->proc~run_part_b 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 SHGLQ_DP(ZERO, W, LMAX)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(IN)  :: LMAX
  REAL(KIND=dp),    INTENT(OUT) :: ZERO(LMAX+1)
  REAL(KIND=dp),    INTENT(OUT) :: W(LMAX+1)
  INTEGER(KIND=i4) :: N, I, J, K, NHALF
  REAL(KIND=dp)    :: X, P0, P1, P2, DPN, DX, TOL
  TOL = 1.0e-15_dp

  N     = LMAX + 1
  NHALF = (N + 1) / 2

  DO I = 1, NHALF
    X = DCOS(pi * (I - 0.25_dp) / (N + 0.5_dp))
    DO J = 1, 100
      P0 = 1.0_dp
      P1 = X
      DO K = 2, N
        P2  = ((2*K - 1) * X * P1 - (K-1) * P0) / K
        P0  = P1
        P1  = P2
      END DO
      DPN = N * (P0 - X * P1) / (1.0_dp - X**2)
      DX  = P1 / DPN
      X   = X - DX
      IF (ABS(DX) .LE. TOL) EXIT
    END DO
    ZERO(I)         = -X
    ZERO(N + 1 - I) =  X
    W(I)            = 2.0_dp / ((1.0_dp - X**2) * DPN**2)
    W(N + 1 - I)    = W(I)
  END DO
  END SUBROUTINE SHGLQ_DP