GAUSS.inc

1       SUBROUTINE GAUSS(A,W,K)
2 C
3 C  A; COSINE OF COLATITUDE
4 C  W; GAUSSIAN WEIGHT
5 C  K; ORDER OF LEGENDRE FUNCTIONS
6 C
7       IMPLICIT REAL*8(A-H,O-Z)
8 C
9       DIMENSION A(K),W(K)
10 C
11       ESP=1.E-14
12       C=(1.E0-(2.E0/pi)**2)*0.25E0
13       FK=K
14       KK=K/2
15       CALL BSSLZ1(A,KK)
16       DO 30 IS=1,KK
17       XZ=COS(A(IS)/SQRT((FK+0.5E0)**2+C))
18       ITER=0
19    10 PKM2=1.0
20       PKM1=XZ
21       ITER=ITER+1
22       IF(ITER.GT.10) GO TO 70
23       DO 20 N=2,K
24       FN=N
25       PK=((2.E0*FN-1.E0)*XZ*PKM1-(FN-1.0)*PKM2)/FN
26       PKM2=PKM1
27    20 PKM1=PK
28       PKM1=PKM2
29       PKMRK=(FK*(PKM1-XZ*PK))/(1.E0-XZ**2)
30       SP=PK/PKMRK
31       XZ=XZ-SP
32       AVSP=ABS(SP)
33       IF(AVSP.GT.ESP) GO TO 10
34       A(IS)=XZ
35       W(IS)=(2.E0*(1.E0-XZ**2))/(FK*PKM1)**2
36    30 CONTINUE
37       IF(K.EQ.KK*2) GO TO 50
38       A(KK+1)=0.E0
39       PK=2.E0/FK**2
40       DO 40 N=2,K,2
41       FN=N
42    40 PK=PK*FN**2/(FN-1.E0)**2
43       W(KK+1)=PK
44    50 CONTINUE
45       DO 60 N=1,KK
46       L=K+1-N
47       A(L)=-A(N)
48    60 W(L)=W(N)
49       RETURN
50    70 WRITE(96,6000)
51  6000 FORMAT(//5X,14HERROR IN GAUAW//)
52 C
53       STOP
54       END SUBROUTINE GAUSS