]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |