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