This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / wpsipg.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:01  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       FUNCTION WPSIPG(Z,K)                                                      
11                                                                                 
12       IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                       
13       COMPLEX*16 WPSIPG,Z,U,V,H,R,P,GCMPLX                                      
14       CHARACTER NAME*(*)                                                        
15       CHARACTER*80 ERRTXT                                                       
16       PARAMETER (NAME = 'CPSIPG/WPSIPG')                                        
17       DIMENSION C(6,0:4),FCT(-1:4),SGN(0:4)                                     
18                                                                                 
19       PARAMETER (DELTA = 5D-13)                                                 
20       PARAMETER (R1 = 1, HF = R1/2)                                             
21       PARAMETER (PI = 3.14159 26535 89793 24D0)                                 
22       PARAMETER (C1 = PI**2, C2 = 2*PI**3, C3 = 2*PI**4, C4 = 8*PI**5)          
23                                                                                 
24       DATA SGN /-1,+1,-1,+1,-1/, FCT /0,1,1,2,6,24/                             
25                                                                                 
26       DATA C(1,0) / 8.33333 33333 33333 33D-2/                                  
27       DATA C(2,0) /-8.33333 33333 33333 33D-3/                                  
28       DATA C(3,0) / 3.96825 39682 53968 25D-3/                                  
29       DATA C(4,0) /-4.16666 66666 66666 67D-3/                                  
30       DATA C(5,0) / 7.57575 75757 57575 76D-3/                                  
31       DATA C(6,0) /-2.10927 96092 79609 28D-2/                                  
32                                                                                 
33       DATA C(1,1) / 1.66666 66666 66666 67D-1/                                  
34       DATA C(2,1) /-3.33333 33333 33333 33D-2/                                  
35       DATA C(3,1) / 2.38095 23809 52380 95D-2/                                  
36       DATA C(4,1) /-3.33333 33333 33333 33D-2/                                  
37       DATA C(5,1) / 7.57575 75757 57575 76D-2/                                  
38       DATA C(6,1) /-2.53113 55311 35531 14D-1/                                  
39                                                                                 
40       DATA C(1,2) / 5.00000 00000 00000 00D-1/                                  
41       DATA C(2,2) /-1.66666 66666 66666 67D-1/                                  
42       DATA C(3,2) / 1.66666 66666 66666 67D-1/                                  
43       DATA C(4,2) /-3.00000 00000 00000 00D-1/                                  
44       DATA C(5,2) / 8.33333 33333 33333 33D-1/                                  
45       DATA C(6,2) /-3.29047 61904 76190 48D+0/                                  
46                                                                                 
47       DATA C(1,3) / 2.00000 00000 00000 00D+0/                                  
48       DATA C(2,3) /-1.00000 00000 00000 00D+0/                                  
49       DATA C(3,3) / 1.33333 33333 33333 33D+0/                                  
50       DATA C(4,3) /-3.00000 00000 00000 00D+0/                                  
51       DATA C(5,3) / 1.00000 00000 00000 00D+1/                                  
52       DATA C(6,3) /-4.60666 66666 66666 67D+1/                                  
53                                                                                 
54       DATA (C(I,4),I=1,6) / 10, -7, 12, -33, 130, -691/                         
55                                                                                 
56       GCMPLX(X,Y)=DCMPLX(X,Y)                                                   
57                                                                                 
58       U=Z                                                                       
59       X=U                                                                       
60       A=ABS(X)                                                                  
61       IF(K .LT. 0 .OR. K .GT. 4) THEN                                           
62        H=0                                                                      
63        WRITE(ERRTXT,101) K                                                      
64        CALL MTLPRT(NAME,'C317.1',ERRTXT)                                        
65       ELSEIF(ABS(IMAG(U)) .LT. DELTA .AND. ABS(X+NINT(A)) .LT. DELTA)           
66      1                                                        THEN              
67        H=0                                                                      
68        WRITE(ERRTXT,102) X                                                      
69        CALL MTLPRT(NAME,'C317.2',ERRTXT)                                        
70       ELSE                                                                      
71        K1=K+1                                                                   
72        IF(X .LT. 0) U=-U                                                        
73        V=U                                                                      
74        H=0                                                                      
75        IF(A .LT. 15) THEN                                                       
76         H=1/V**K1                                                               
77         DO 1 I = 1,14-INT(A)                                                    
78         V=V+1                                                                   
79     1   H=H+1/V**K1                                                             
80         V=V+1                                                                   
81        END IF                                                                   
82        R=1/V**2                                                                 
83        P=R*C(6,K)                                                               
84        DO 2 I = 5,1,-1                                                          
85     2  P=R*(C(I,K)+P)                                                           
86        H=SGN(K)*(FCT(K)*H+(V*(FCT(K-1)+P)+HF*FCT(K))/V**K1)                     
87        IF(K .EQ. 0) H=H+LOG(V)                                                  
88        IF(X .LT. 0) THEN                                                        
89         V=PI*U                                                                  
90         X=V                                                                     
91         Y=IMAG(V)                                                               
92         A=SIN(X)                                                                
93         B=COS(X)                                                                
94         T=TANH(Y)                                                               
95         P=GCMPLX(B,-A*T)/GCMPLX(A,B*T)                                          
96         IF(K .EQ. 0) THEN                                                       
97          H=H+1/U+PI*P                                                           
98         ELSEIF(K .EQ. 1) THEN                                                   
99          H=-H+1/U**2+C1*(P**2+1)                                                
100         ELSEIF(K .EQ. 2) THEN                                                   
101          H=H+2/U**3+C2*P*(P**2+1)                                               
102         ELSEIF(K .EQ. 3) THEN                                                   
103          R=P**2                                                                 
104          H=-H+6/U**4+C3*((3*R+4)*R+1)                                           
105         ELSEIF(K .EQ. 4) THEN                                                   
106          R=P**2                                                                 
107          H=H+24/U**5+C4*P*((3*R+5)*R+2)                                         
108         ENDIF                                                                   
109        ENDIF                                                                    
110       ENDIF                                                                     
111       WPSIPG=H                                                                  
112       RETURN                                                                    
113   101 FORMAT('K = ',I5,'  (< 0  OR  > 4)')                                      
114   102 FORMAT('ARGUMENT EQUALS NON-POSITIVE INTEGER = ',F8.1)                    
115       END