This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / bsja128.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:08  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if defined(CERNLIB_QUAD)
11 #if defined(CERNLIB_DOUBLE)
12       SUBROUTINE QBSJA(X,A,NMAX,ND,B)
13 #endif
14 #if !defined(CERNLIB_DOUBLE)
15       SUBROUTINE DBSJA(X,A,NMAX,ND,B)
16 #endif
17
18 #include "gen/imp128.inc"
19       REAL SX,D,T,Q,U,V,TC(11)
20       CHARACTER*80 ERRTXT
21       CHARACTER NAMEJ*(*),NAMEI*(*)
22 #if defined(CERNLIB_DOUBLE)
23       PARAMETER (NAMEJ = 'DBSJA/QBSJA',
24      1           NAMEI = 'DBSIA/QBSIA')
25 #endif
26 #if !defined(CERNLIB_DOUBLE)
27       PARAMETER (NAMEJ = 'DBSJA', NAMEI = 'DBSIA')
28 #endif
29       LOGICAL LJA,LIA,LEV,LER
30       DIMENSION B(0:*),BA(0:100),RR(0:100)
31
32       PARAMETER (Z1 = 1, HF = Z1/2, Z10 = 10)
33
34       DATA TC / 5.7941 E-5,-1.76148E-3, 2.08645E-2,-1.29013E-1,
35      1          8.5777 E-1, 1.0125 E+0, 7.75   E-1, 2.3026 E+0,
36      2          1.3863 E+0, 7.3576 E-1, 1.3591 E+0/
37
38       LJA=.TRUE.
39       LIA=.FALSE.
40       SGN=-1
41       GO TO 9
42
43 #if defined(CERNLIB_DOUBLE)
44       ENTRY QBSIA(X,A,NMAX,ND,B)
45 #endif
46 #if !defined(CERNLIB_DOUBLE)
47       ENTRY DBSIA(X,A,NMAX,ND,B)
48 #endif
49       LJA=.FALSE.
50       LIA=.TRUE.
51       SGN=1
52
53     9 LER=.FALSE.
54       IF(X .LE. 0) THEN
55        WRITE(ERRTXT,101) X
56        IF(LJA) CALL MTLPRT(NAMEJ,'C343.1',ERRTXT)
57        IF(LIA) CALL MTLPRT(NAMEI,'C343.1',ERRTXT)
58        LER=.TRUE.
59       ELSEIF(.NOT.(0 .LE. A .AND. A .LT. 1)) THEN
60        WRITE(ERRTXT,102) A
61        IF(LJA) CALL MTLPRT(NAMEJ,'C343.2',ERRTXT)
62        IF(LIA) CALL MTLPRT(NAMEI,'C343.2',ERRTXT)
63        LER=.TRUE.
64       ELSEIF(ABS(NMAX) .GT. 100) THEN
65        WRITE(ERRTXT,103) NMAX
66        IF(LJA) CALL MTLPRT(NAMEJ,'C343.3',ERRTXT)
67        IF(LIA) CALL MTLPRT(NAMEI,'C343.3',ERRTXT)
68        LER=.TRUE.
69       END IF
70       IF(LER) RETURN
71       EPS=HF*Z10**(-ND)
72       NMX=ABS(NMAX)
73       IF(NMAX .LE. 0) NMX=1
74       DO 5 N = 0,NMX
75       RR(N)=0
76     5 BA(N)=0
77       D=TC(8)*ND+TC(9)
78       SX=X
79       Q=0
80       IF(NMX .GT. 0) THEN
81        V=0.5*D/NMX
82        IF(V .LE. 10) THEN
83         T=TC(1)
84         DO 6 I = 2,6
85     6   T=V*T+TC(I)
86        ELSE
87         U=LOG(V)-TC(7)
88         T=V/(U*(1+(TC(7)-LOG(U))/(1+U)))
89        ENDIF
90        Q=NMX*T
91       ENDIF
92 #if defined(CERNLIB_DOUBLE)
93       F=(HF*X)**A/QGAMMA(1+A)
94 #endif
95 #if !defined(CERNLIB_DOUBLE)
96       F=(HF*X)**A/DGAMMA(1+A)
97 #endif
98       T=1
99       V=TC(10)*D/SX
100       IF(LIA) THEN
101        F=EXP(X)*F
102        V=V-TC(10)
103       ENDIF
104       IF(LJA .OR. LIA .AND. X .LT. D) THEN
105        IF(V .LE. 10) THEN
106         T=TC(1)
107         DO 7 I = 2,6
108     7   T=V*T+TC(I)
109        ELSE
110         U=LOG(V)-TC(7)
111         T=V/(U*(1+(TC(7)-LOG(U))/(1+U)))
112        ENDIF
113       ENDIF
114       NU=1+MAX(Q,TC(11)*SX*T)
115
116       MU=-1
117     2 MU=MU+1
118       AL=1
119       IF(LJA) THEN
120        DO 3 N = 1,NU/2
121        XN=N
122     3  AL=AL*(XN+A)/(XN+1)
123        R=0
124        S=0
125        LEV=.TRUE.
126        DO 4 N = 2*(NU/2),1,-1
127        XN=N
128        XA=XN+A
129        R=1/(2*XA/X-R)
130        IF(N .LE. NMX) RR(N-1)=R
131        IF(LEV) THEN
132         AL=AL*(XN+2)/(XA+A)
133         S=R*(AL*XA+S)
134        ELSE
135         S=R*S
136        ENDIF
137        LEV=.NOT.LEV
138     4  CONTINUE
139       ELSE
140        DO 23 N = 1,NU
141        XN=N
142    23  AL=AL*(XN+2*A)/(XN+1)
143        R=0
144        S=0
145        DO 24 N = NU,1,-1
146        XN=N
147        XA=XN+A
148        XA2=XA+XA
149        R=1/(XA2/X+R)
150        IF(N .LE. NMX) RR(N-1)=R
151        AL=AL*(XN+1)/(XA+A)
152        S=R*(XA2*AL+S)
153    24  CONTINUE
154       ENDIF
155       B(0)=F/(1+S)
156       DO 10 N = 0,NMX-1
157    10 B(N+1)=RR(N)*B(N)
158       DO 11 N = 0,NMX
159       IF(ABS(B(N)-BA(N)) .GT. EPS*ABS(B(N))) THEN
160        DO 12 M = 0,NMX
161    12  BA(M)=B(M)
162        NU=NU+5
163        IF(MU .LE. 50) GO TO 2
164        WRITE(ERRTXT,104) X,A
165        IF(LJA) CALL MTLPRT(NAMEJ,'C343.4',ERRTXT)
166        IF(LIA) CALL MTLPRT(NAMEI,'C343.4',ERRTXT)
167        RETURN
168       ENDIF
169    11 CONTINUE
170       IF(NMAX .LT. 0) THEN
171        AL=2/X
172        B(1)=AL*A*B(0)+SGN*B(1)
173        DO 13 N = 1,-NMAX-1
174    13  B(N+1)=AL*(A-N)*B(N)+SGN*B(N-1)
175       ENDIF
176       RETURN
177   101 FORMAT('ILLEGAL ARGUMENT X = ',1P,D15.8)
178   102 FORMAT('ILLEGAL ORDER A = ',1P,D15.8)
179   103 FORMAT('ILLEGAL NMAX = ',I5)
180   104 FORMAT('NO CONVERGENCE FOR X = ',1P,D15.8,' A = ',D15.8,
181      1       ' TRY SMALLER ND')
182       END
183 #endif