ASSOC_LEGENDRE_DP Function

private function ASSOC_LEGENDRE_DP(L, K, X) result(ASSOC_LEGENDRE)

Associated Legendre function , Condon-Shortley phase convention, computed by the classic Bonnet-type upward recurrence starting from . For negative order, . Unnormalized (unlike ASSOC_LEGENDRE_NORM_ALL); loses numerical stability well before the intermediate values literally overflow, since the later upward-in- recurrence involves a subtraction of two comparably huge terms (catastrophic cancellation).

Warning

Numerical stability: verified accurate (relative error < 3e-14) for at any -- the worst case, , occurs near the equator ( degrees) at ; safety improves toward either pole. This supersedes an earlier "" claim, verified too optimistic by direct comparison against an arbitrary-precision reference. Full methodology and per-degree table: STABILITY_FINDINGS.md at the repo root, Part 2 (regenerable via the opt-in VSH_BUILD_STABILITY CMake target's vsh_unnorm_stability executable).

Returns: .

Arguments

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

Degree, .
integer(kind=i4), intent(in) :: K

Order, .
real(kind=dp), intent(in) :: X

Argument, (typically ).

Return Value real(kind=dp)


Calls

proc~~assoc_legendre_dp~~CallsGraph proc~assoc_legendre_dp ASSOC_LEGENDRE_DP interface~log_fact LOG_FACT proc~assoc_legendre_dp->interface~log_fact proc~log_fact_dp LOG_FACT_DP interface~log_fact->proc~log_fact_dp

Called by

proc~~assoc_legendre_dp~~CalledByGraph proc~assoc_legendre_dp ASSOC_LEGENDRE_DP interface~assoc_legendre ASSOC_LEGENDRE interface~assoc_legendre->proc~assoc_legendre_dp proc~batch_alm_cons BATCH_ALM_CONS proc~batch_alm_cons->interface~assoc_legendre proc~bench_legendre BENCH_LEGENDRE proc~bench_legendre->interface~assoc_legendre proc~ssh_dp SSH_DP proc~ssh_dp->interface~assoc_legendre program~convergence_main CONVERGENCE_MAIN program~convergence_main->interface~assoc_legendre program~convergence_main->proc~batch_alm_cons interface~ssh SSH program~convergence_main->interface~ssh proc~batch_ssh_cons BATCH_SSH_CONS program~convergence_main->proc~batch_ssh_cons program~unnorm_stability_main UNNORM_STABILITY_MAIN program~unnorm_stability_main->interface~assoc_legendre interface~ssh->proc~ssh_dp program~bench_main BENCH_MAIN program~bench_main->proc~bench_legendre proc~bench_ssh BENCH_SSH program~bench_main->proc~bench_ssh program~main MAIN program~main->proc~batch_alm_cons program~main->proc~batch_ssh_cons proc~test1_gw TEST1_GW program~main->proc~test1_gw proc~test2_gw TEST2_GW program~main->proc~test2_gw proc~batch_pvsh_rad_cons BATCH_PVSH_RAD_CONS program~main->proc~batch_pvsh_rad_cons proc~test2_vsh TEST2_VSH program~main->proc~test2_vsh proc~batch_ssh_cons->interface~ssh proc~bench_ssh->interface~ssh proc~pvsh_rad_dp PVSH_RAD_DP proc~pvsh_rad_dp->interface~ssh proc~run_gwi_sweep RUN_GWI_SWEEP proc~run_gwi_sweep->interface~ssh proc~run_gwj_sweep RUN_GWJ_SWEEP proc~run_gwj_sweep->interface~ssh proc~test1_gw->interface~ssh proc~test2_gw->interface~ssh program~dipole_synthesis DIPOLE_SYNTHESIS program~dipole_synthesis->interface~ssh interface~pvsh_rad PVSH_RAD program~dipole_synthesis->interface~pvsh_rad interface~pvsh_rad->proc~pvsh_rad_dp program~gwi_gwj_sweep GWI_GWJ_SWEEP program~gwi_gwj_sweep->proc~run_gwi_sweep program~gwi_gwj_sweep->proc~run_gwj_sweep proc~batch_pvsh_rad_cons->interface~pvsh_rad proc~field_known FIELD_KNOWN proc~field_known->interface~pvsh_rad proc~test2_vsh->interface~pvsh_rad program~dipole_flux_expulsion DIPOLE_FLUX_EXPULSION program~dipole_flux_expulsion->interface~pvsh_rad program~dipole_uniform_sphere DIPOLE_UNIFORM_SPHERE program~dipole_uniform_sphere->interface~pvsh_rad

Source Code

  FUNCTION ASSOC_LEGENDRE_DP(L, K, X) RESULT(ASSOC_LEGENDRE)
  IMPLICIT NONE
  INTEGER(KIND=i4), INTENT(IN) :: L
  INTEGER(KIND=i4), INTENT(IN) :: K
  REAL(KIND=dp), INTENT(IN) :: X
  INTEGER(KIND=i4) :: ABSK, I
  REAL(KIND=dp) :: ASSOC_LEGENDRE
  REAL(KIND=dp) :: Pmm, Pmm1, Plm, somx2

  ABSK = ABS(K)

  IF (ABSK > L .OR. ABS(X) > 1.0_dp) THEN
    ASSOC_LEGENDRE = 0.0_dp
    RETURN
  END IF

  Pmm = 1.0_dp
  IF (ABSK > 0) THEN
    somx2 = sqrt((1.0_dp - X) * (1.0_dp + X))
    DO I = 1, ABSK
      Pmm = -Pmm * (2*I - 1) * somx2
    END DO
  END IF

  IF (L == ABSK) THEN
    Plm = Pmm
  ELSE
    Pmm1 = X * (2 * ABSK + 1) * Pmm
    IF (L == ABSK + 1) THEN
      Plm = Pmm1
    ELSE
      DO I = ABSK + 2, L
        Plm = (X*(2*I-1)*Pmm1-(I+ABSK-1)*Pmm) / (I-ABSK)
        Pmm = Pmm1
        Pmm1 = Plm
      END DO
    END IF
  END IF

  IF (K < 0) THEN
    Plm = ((-1.0_dp)**ABSK)* &
      EXP(LOG_FACT(L-ABSK)-LOG_FACT(L+ABSK))*Plm
  END IF

  ASSOC_LEGENDRE = Plm
  END FUNCTION ASSOC_LEGENDRE_DP