source: trunk/fsource/texturesubs/qlmn.for @ 305

Last change on this file since 305 was 305, checked in by vondreele, 10 years ago

remove pypowder & powsubs - no longer needed
add texturesubs & pytexture

File size: 2.0 KB
Line 
1      SUBROUTINE QLMN(L,MM,NN,Q)
2
3!PURPOSE: Compute Ql,m,n for spherical harmonics from lookup table
4
5      INCLUDE       '../INCLDS/COPYRIGT.FOR' 
6
7!CALLING ARGUMENTS:         
8
9      INTEGER*4     L,MM,NN             !order & subindices (m may be <0)
10      REAL*4        Q                   !Output value
11
12!INCLUDE STATEMENTS:
13
14      REAL*4        QT                 
15      COMMON /QLMNVAL/QT(2109)
16
17
18!LOCAL VARIABLES:
19
20      REAL*8        SUM,TEMP,TEMP1     
21      INTEGER*4     LMN,I,J,M,N         
22
23!FUNCTION DEFINITIONS:                     
24
25      REAL*8        FACTLN              !Compute ln-factorial & binominal coeffs.
26
27!DATA STATEMENTS:
28
29!CODE: 
30                         
31      M = ABS(MM)           
32      N = ABS(NN)
33      IF ( N.GT.M ) THEN
34        I = M
35        M = N
36        N = I
37      END IF
38      IF ( MOD(N,2).EQ.0 .AND. MOD(L,2).EQ.0 ) THEN         !Even L,N - do lookup
39        J = 0
40        DO I=2,L,2
41          J = J+(I/2)**2      !points to last term in L-2 block
42        END DO
43        J = J+1
44        DO I=0,M-1
45          J = J+(I+2)/2       !points to 1st term in M block
46        END DO
47        J = J+N/2             !offset to N term
48!       PRINT '(A,I4,F12.8)',' J,Q ',J,QT(J)
49        Q = QT(J)                         
50      ELSE                       !Odd L or N - calculate Q
51        LMN = L-M-N
52        TEMP = 0.5D0*(FACTLN(L+N)+FACTLN(L+M)+
53     1    FACTLN(L-M)+FACTLN(L-N))
54        SUM = 0.0D0
55        DO I=0,L-N
56          IF ( (L-M-I).GE.0 .AND. (M+N+I).GE.0 ) THEN
57            TEMP1 = TEMP-FACTLN(I)-FACTLN(L-M-I)-
58     1        FACTLN(L-N-I)-FACTLN(M+N+I)
59            TEMP1 = DEXP(TEMP1)
60            IF ( MOD(I,2).NE.0 ) TEMP1 = -TEMP1
61            SUM = SUM+TEMP1
62          END IF
63        END DO
64        Q = SUM/2.0**L
65        IF ( MOD(LMN,2).NE.0 ) Q = -Q
66!       PRINT '(A,3I4,F12.8)',' l,m,n,Q(lmn)',L,M,N,Q     
67      END IF
68      IF ( MM.LT.0 .AND. MOD(L+NN,2).NE.0 ) Q = -Q
69      IF ( NN.LT.0 .AND. MOD(L+MM,2).NE.0 ) Q = -Q
70      RETURN
71      END
72
Note: See TracBrowser for help on using the repository browser.