diff options
| author | stainer_t <thomas.stainer@oecd-nea.org> | 2025-09-08 13:48:49 +0200 |
|---|---|---|
| committer | stainer_t <thomas.stainer@oecd-nea.org> | 2025-09-08 13:48:49 +0200 |
| commit | 7dfcc480ba1e19bd3232349fc733caef94034292 (patch) | |
| tree | 03ee104eb8846d5cc1a981d267687a729185d3f3 /Utilib/src/ALPLGN.f | |
Initial commit from Polytechnique Montreal
Diffstat (limited to 'Utilib/src/ALPLGN.f')
| -rw-r--r-- | Utilib/src/ALPLGN.f | 63 |
1 files changed, 63 insertions, 0 deletions
diff --git a/Utilib/src/ALPLGN.f b/Utilib/src/ALPLGN.f new file mode 100644 index 0000000..8c1c030 --- /dev/null +++ b/Utilib/src/ALPLGN.f @@ -0,0 +1,63 @@ +*DECK ALPLGN + DOUBLE PRECISION FUNCTION ALPLGN(L,M,X) +* +*----------------------------------------------------------------------- +* +*Purpose: +* return the Ferrer definition of the associated Legendre function. +* +*Copyright: +* Copyright (C) 2021 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): A. Hebert +* +*Parameters: input +* L main index +* M secondary index +* X direction cosine +* +*----------------------------------------------------------------------- +* + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + INTEGER L,M + DOUBLE PRECISION X + INTEGER, ALLOCATABLE, DIMENSION(:) :: IPX +* + IF(M.LT.0) THEN + CALL XABORT('ALPLGN: BAD ARGUMENTS (1)') + ELSE IF(M.GT.L) THEN + CALL XABORT('ALPLGN: BAD ARGUMENTS (2)') + ELSE IF(ABS(X).GT.1.0D0) THEN + CALL XABORT('ALPLGN: BAD ARGUMENTS (3)') + ENDIF + PMM=1.0D0 + IF(M.GT.0) THEN + ALLOCATE(IPX(2*M)) + DO I=1,2*M + IPX(I)=I + ENDDO + PMM=PRODUCT(IPX,MASK=MOD(IPX,2)==1)*SQRT((1.0D0-X)*(1.0D0+X))**M + DEALLOCATE(IPX) + ENDIF + IF(L.EQ.M) THEN + ALPLGN=PMM + ELSE + PMMP1=(2*M+1)*X*PMM + IF(L.EQ.M+1) THEN + ALPLGN=PMMP1 + ELSE + PLL=0.0D0 + DO LL=M+2,L + PLL=((2*LL-1)*X*PMMP1-(LL+M-1)*PMM)/(LL-M) + PMM=PMMP1 + PMMP1=PLL + ENDDO + ALPLGN=PLL + ENDIF + ENDIF + RETURN + END |