MINICERN/mathlib/gen/c/bsir464.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:03  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       FUNCTION  BSIR4(X,NU)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       FUNCTION DBSIR4(X,NU)
  15 #include "gen/imp64.inc"
  16 #endif
  17       CHARACTER*80 ERRTXT
  18       CHARACTER NAMEI*(*),NAMEK*(*),NAMEIE*(*),NAMEKE*(*)
  19 #if !defined(CERNLIB_DOUBLE)
  20       PARAMETER (NAMEI = 'BSIR4/DBSIR4', NAMEIE = 'EBSIR4/DEBIR4')
  21       PARAMETER (NAMEK = 'BSKR4/DBSKR4', NAMEKE = 'EBSKR4/DEBKR4')
  22 #endif
  23 #if defined(CERNLIB_DOUBLE)
  24       PARAMETER (NAMEI = 'BSIR4/DBSIR4', NAMEIE = 'EBSIR4/DEBIR4')
  25       PARAMETER (NAMEK = 'BSKR4/DBSKR4', NAMEKE = 'EBSKR4/DEBKR4')
  26 #endif
  27       LOGICAL LEX
  28
  29       DIMENSION BC(0:23,2),CC(0:15,2),PP(-3:3),GG(-3:3)
  30
  31       PARAMETER (EPS = 1D-14)
  32       PARAMETER (Z1 = 1, HF =Z1/2)
  33       PARAMETER (PI = 3.14159 26535 89793 24D0)
  34       PARAMETER (W2 = 1.41421 35623 73095 05D0)
  35       PARAMETER (G1 = 3.62560 99082 21908 31D0)
  36       PARAMETER (G3 = 1.22541 67024 65177 65D0)
  37       PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI = 1/PI, RW2 = 1/W2)
  38       PARAMETER (C1 = PI/(2*W2))
  39       PARAMETER (GM = 2*(1/G3-4/G1), GP = (4/G1+1/G3)/2)
  40
  41       DATA GG(-3) /0.27581 56628 30209 31D0/, PP(-3) /-0.75D0/
  42       DATA GG(-1) /0.81604 89390 98262 98D0/, PP(-1) /-0.25D0/
  43       DATA GG( 1) /1.10326 26513 20837 26D0/, PP( 1) / 0.25D0/
  44       DATA GG( 3) /1.08806 52521 31017 31D0/, PP( 3) / 0.75D0/
  45
  46       DATA BC( 0,1) / 1.00619 92270 14122 57D0/
  47       DATA BC( 1,1) / 0.00631 99620 31140 72D0/
  48       DATA BC( 2,1) / 0.00012 56131 27965 64D0/
  49       DATA BC( 3,1) / 0.00000 52052 40761 57D0/
  50       DATA BC( 4,1) / 0.00000 03591 84411 39D0/
  51       DATA BC( 5,1) / 0.00000 00355 85362 89D0/
  52       DATA BC( 6,1) / 0.00000 00036 05011 66D0/
  53       DATA BC( 7,1) /-0.00000 00001 26294 10D0/
  54       DATA BC( 8,1) /-0.00000 00002 96595 12D0/
  55       DATA BC( 9,1) /-0.00000 00001 18337 70D0/
  56       DATA BC(10,1) /-0.00000 00000 21655 68D0/
  57       DATA BC(11,1) / 0.00000 00000 03032 04D0/
  58       DATA BC(12,1) / 0.00000 00000 03041 10D0/
  59       DATA BC(13,1) / 0.00000 00000 00530 77D0/
  60       DATA BC(14,1) /-0.00000 00000 00204 53D0/
  61       DATA BC(15,1) /-0.00000 00000 00105 49D0/
  62       DATA BC(16,1) / 0.00000 00000 00005 50D0/
  63       DATA BC(17,1) / 0.00000 00000 00014 36D0/
  64       DATA BC(18,1) / 0.00000 00000 00001 14D0/
  65       DATA BC(19,1) /-0.00000 00000 00001 87D0/
  66       DATA BC(20,1) /-0.00000 00000 00000 32D0/
  67       DATA BC(21,1) / 0.00000 00000 00000 26D0/
  68       DATA BC(22,1) / 0.00000 00000 00000 06D0/
  69       DATA BC(23,1) /-0.00000 00000 00000 04D0/
  70
  71       DATA BC( 0,2) / 0.98980 19115 24008 91D0/
  72       DATA BC( 1,2) /-0.01035 09365 14827 02D0/
  73       DATA BC( 2,2) /-0.00015 85263 84973 08D0/
  74       DATA BC( 3,2) /-0.00000 60527 21962 69D0/
  75       DATA BC( 4,2) /-0.00000 04158 38597 31D0/
  76       DATA BC( 5,2) /-0.00000 00487 99346 57D0/
  77       DATA BC( 6,2) /-0.00000 00089 86835 44D0/
  78       DATA BC( 7,2) /-0.00000 00019 83283 58D0/
  79       DATA BC( 8,2) /-0.00000 00003 58969 60D0/
  80       DATA BC( 9,2) /-0.00000 00000 08766 62D0/
  81       DATA BC(10,2) / 0.00000 00000 25819 45D0/
  82       DATA BC(11,2) / 0.00000 00000 09780 24D0/
  83       DATA BC(12,2) / 0.00000 00000 00565 05D0/
  84       DATA BC(13,2) /-0.00000 00000 00851 66D0/
  85       DATA BC(14,2) /-0.00000 00000 00270 25D0/
  86       DATA BC(15,2) / 0.00000 00000 00040 96D0/
  87       DATA BC(16,2) / 0.00000 00000 00040 50D0/
  88       DATA BC(17,2) / 0.00000 00000 00001 11D0/
  89       DATA BC(18,2) /-0.00000 00000 00005 25D0/
  90       DATA BC(19,2) /-0.00000 00000 00000 70D0/
  91       DATA BC(20,2) / 0.00000 00000 00000 70D0/
  92       DATA BC(21,2) / 0.00000 00000 00000 14D0/
  93       DATA BC(22,2) /-0.00000 00000 00000 10D0/
  94       DATA BC(23,2) /-0.00000 00000 00000 02D0/
  95
  96       DATA CC( 0,1) / 0.99128 81656 75147 07D0/
  97       DATA CC( 1,1) /-0.00850 62567 20022 24D0/
  98       DATA CC( 2,1) / 0.00019 70491 57408 35D0/
  99       DATA CC( 3,1) /-0.00000 80377 10166 54D0/
 100       DATA CC( 4,1) / 0.00000 04554 01498 43D0/
 101       DATA CC( 5,1) /-0.00000 00323 27352 82D0/
 102       DATA CC( 6,1) / 0.00000 00027 16130 28D0/
 103       DATA CC( 7,1) /-0.00000 00002 60644 07D0/
 104       DATA CC( 8,1) / 0.00000 00000 27882 69D0/
 105       DATA CC( 9,1) /-0.00000 00000 03267 69D0/
 106       DATA CC(10,1) / 0.00000 00000 00414 09D0/
 107       DATA CC(11,1) /-0.00000 00000 00056 17D0/
 108       DATA CC(12,1) / 0.00000 00000 00008 09D0/
 109       DATA CC(13,1) /-0.00000 00000 00001 23D0/
 110       DATA CC(14,1) / 0.00000 00000 00000 20D0/
 111       DATA CC(15,1) /-0.00000 00000 00000 03D0/
 112
 113       DATA CC( 0,2) / 1.01476 24350 64637 87D0/
 114       DATA CC( 1,2) / 0.01449 34617 87809 66D0/
 115       DATA CC( 2,2) /-0.00025 87162 07241 80D0/
 116       DATA CC( 3,2) / 0.00000 96912 18911 49D0/
 117       DATA CC( 4,2) /-0.00000 05261 29313 99D0/
 118       DATA CC( 5,2) / 0.00000 00363 96854 29D0/
 119       DATA CC( 6,2) /-0.00000 00030 05472 76D0/
 120       DATA CC( 7,2) / 0.00000 00002 84827 80D0/
 121       DATA CC( 8,2) /-0.00000 00000 30182 91D0/
 122       DATA CC( 9,2) / 0.00000 00000 03511 10D0/
 123       DATA CC(10,2) /-0.00000 00000 00442 27D0/
 124       DATA CC(11,2) / 0.00000 00000 00059 70D0/
 125       DATA CC(12,2) /-0.00000 00000 00008 56D0/
 126       DATA CC(13,2) / 0.00000 00000 00001 30D0/
 127       DATA CC(14,2) /-0.00000 00000 00000 21D0/
 128       DATA CC(15,2) / 0.00000 00000 00000 03D0/
 129
 130       LEX=.FALSE.
 131       GO TO 8
 132
 133 #if !defined(CERNLIB_DOUBLE)
 134       ENTRY EBSIR4(X,NU)
 135 #endif
 136 #if defined(CERNLIB_DOUBLE)
 137       ENTRY DEBIR4(X,NU)
 138 #endif
 139       LEX=.TRUE.
 140
 141     8 MU=ABS(NU)
 142       IF(MU .NE. 1 .AND. MU .NE. 2 .AND. MU .NE. 3 .OR.
 143      1   NU .LT. 0 .AND. X .LE. 0 .OR. NU .GT. 0 .AND. X .LT. 0) THEN
 144        S=0
 145        WRITE(ERRTXT,101) X,NU
 146        IF(.NOT.LEX) CALL MTLPRT(NAMEI ,'C327.1',ERRTXT)
 147        IF(     LEX) CALL MTLPRT(NAMEIE,'C327.1',ERRTXT)
 148       ELSEIF(X .EQ. 0) THEN
 149        S=0
 150       ELSEIF(NU .EQ. -2) THEN
 151        IF(LEX) THEN
 152         S=HF*(1+EXP(-X-X))/SQRT(PIH*X)
 153        ELSE
 154         S=COSH(X)/SQRT(PIH*X)
 155        ENDIF
 156       ELSEIF(NU .EQ. 2) THEN
 157        IF(LEX) THEN
 158         IF(X .LT. HF) THEN
 159          S=SINH(X)*EXP(-X)/SQRT(PIH*X)
 160         ELSE
 161          S=HF*(1-EXP(-X-X))/SQRT(PIH*X)
 162         ENDIF
 163        ELSE
 164         S=SINH(X)/SQRT(PIH*X)
 165        ENDIF
 166       ELSEIF(X .LT. 8) THEN
 167        Y=(HF*X)**2
 168        XN=PP(NU)
 169        XL=XN+2
 170        A0=1
 171        A1=1+2*Y/((XL+1)*(XN+1))
 172        A2=1+Y*(4+3*Y/((XL+2)*(XN+2)))/((XL+3)*(XN+1))
 173        B0=1
 174        B1=1-Y/(XL+1)
 175        B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
 176        T1=3+XL
 177        V1=3-XL
 178        V3=XL-1
 179        V2=V3+V3
 180        C=0
 181        DO 33 N = 3,30
 182        C0=C
 183        T1=T1+2
 184        T2=T1-1
 185        T3=T2-1
 186        T4=T3-1
 187        T5=T4-1
 188        T6=T5-1
 189        V1=V1+1
 190        V2=V2+1
 191        V3=V3+1
 192        U1=N*T4
 193        E=V3/(U1*T3)
 194        U2=E*Y
 195        F1=1+Y*V1/(U1*T1)
 196        F2=(1+Y*V2/(V3*T2*T5))*U2
 197        F3=-Y*Y*U2/(T4*T5*T5*T6)
 198        A=F1*A2+F2*A1+F3*A0
 199        B=F1*B2+F2*B1+F3*B0
 200        C=A/B
 201        IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 34
 202        A0=A1
 203        A1=A2
 204        A2=A
 205        B0=B1
 206        B1=B2
 207        B2=B
 208    33  CONTINUE
 209    34  S=GG(NU)*(HF*X)**PP(NU)*C
 210        IF(LEX) S=EXP(-X)*S
 211       ELSE
 212        K=(MU+1)/2
 213        R=1/X
 214        W=SQRT(RPI*R)
 215        H=16*R-1
 216        ALFA=H+H
 217        B1=0
 218        B2=0
 219        DO 2 I = 23,0,-1
 220        B0=BC(I,K)+ALFA*B1-B2
 221        B2=B1
 222     2  B1=B0
 223        S=RW2*W*(B0-H*B2)
 224        IF(.NOT.LEX) S=EXP(X)*S
 225        T=0
 226        IF(NU .LT. 0) THEN
 227         H=10*R-1
 228         ALFA=H+H
 229         B1=0
 230         B2=0
 231         DO 3 I = 15,0,-1
 232         B0=CC(I,K)+ALFA*B1-B2
 233         B2=B1
 234     3   B1=B0
 235         R=EXP(-X)
 236         T=W*R*(B0-H*B2)
 237         IF(LEX) T=R*T
 238        ENDIF
 239        S=S+T
 240       ENDIF
 241       GO TO 99
 242
 243 #if !defined(CERNLIB_DOUBLE)
 244       ENTRY  BSKR4(X,NU)
 245 #endif
 246 #if defined(CERNLIB_DOUBLE)
 247       ENTRY DBSKR4(X,NU)
 248 #endif
 249       LEX=.FALSE.
 250       GO TO 9
 251
 252 #if !defined(CERNLIB_DOUBLE)
 253       ENTRY EBSKR4(X,NU)
 254 #endif
 255 #if defined(CERNLIB_DOUBLE)
 256       ENTRY DEBKR4(X,NU)
 257 #endif
 258       LEX=.TRUE.
 259
 260     9 MU=ABS(NU)
 261       IF(MU .NE. 1 .AND. MU .NE. 2 .AND. MU .NE. 3 .OR. X .LE. 0) THEN
 262        S=0
 263        WRITE(ERRTXT,101) X,NU
 264        IF(.NOT.LEX) CALL MTLPRT(NAMEK ,'C327.1',ERRTXT)
 265        IF(     LEX) CALL MTLPRT(NAMEKE,'C327.1',ERRTXT)
 266       ELSEIF(MU .EQ. 2) THEN
 267        S=SQRT(PIH/X)
 268        IF(.NOT.LEX) S=EXP(-X)*S
 269       ELSEIF(X .LE. 1) THEN
 270        A0=PP(-1)
 271        B=HF*X
 272        D=-LOG(B)
 273        F=A0*D
 274        E=EXP(F)
 275        G=(GM*A0+GP)*E
 276        BK=C1*(HF*GM*(E+1/E)+GP*D*SINH(F)/F)
 277        F=BK
 278        E=A0**2
 279        P=HF*C1*G
 280        Q=HF/G
 281        C=1
 282        D=B**2
 283        BK1=P
 284        DO 11 N = 1,15
 285        FN=N
 286        F=(FN*F+P+Q)/(FN**2-E)
 287        C=C*D/FN
 288        P=P/(FN-A0)
 289        Q=Q/(FN+A0)
 290        G=C*(P-FN*F)
 291        H=C*F
 292        BK=BK+H
 293        BK1=BK1+G
 294        IF(H*BK1+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
 295    11  CONTINUE
 296    12  S=BK
 297        IF(MU .EQ. 3) S=BK1/B
 298        IF(LEX) S=EXP(X)*S
 299       ELSEIF(X .LE. 5) THEN
 300        XN=4*PP(MU)**2
 301        A=9-XN
 302        B=25-XN
 303        C=768*X**2
 304        C0=48*X
 305        A0=1
 306        A1=(16*X+7+XN)/A
 307        A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
 308        B0=1
 309        B1=(16*X+9-XN)/A
 310        B2=(C+C0*B)/(A*B)+1
 311        C=0
 312        DO 24 N = 3,30
 313        C0=C
 314        FN=N
 315        FN2=FN+FN
 316        FN1=FN2-1
 317        FN3=FN1/(FN2-3)
 318        FN4=12*FN**2-(1-XN)
 319        FN5=16*FN1*X
 320        RAN=1/((FN2+1)**2-XN)
 321        F1=FN3*(FN4-20*FN)+FN5
 322        F2=28*FN-FN4-8+FN5
 323        F3=FN3*((FN2-5)**2-XN)
 324        A=(F1*A2+F2*A1+F3*A0)*RAN
 325        B=(F1*B2+F2*B1+F3*B0)*RAN
 326        C=A/B
 327        IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
 328        A0=A1
 329        A1=A2
 330        A2=A
 331        B0=B1
 332        B1=B2
 333        B2=B
 334    24  CONTINUE
 335    25  S=C/SQRT(RPIH*X)
 336        IF(.NOT.LEX) S=EXP(-X)*S
 337       ELSE
 338        K=(MU+1)/2
 339        R=1/X
 340        Y=5*R
 341        H=10*R-1
 342        ALFA=H+H
 343        B1=0
 344        B2=0
 345        DO 13 I = 15,0,-1
 346        B0=CC(I,K)+ALFA*B1-B2
 347        B2=B1
 348    13  B1=B0
 349        S=SQRT(PIH*R)*(B0-H*B2)
 350        IF(.NOT.LEX) S=EXP(-X)*S
 351       ENDIF
 352 #if !defined(CERNLIB_DOUBLE)
 353    99  BSIR4=S
 354 #endif
 355 #if defined(CERNLIB_DOUBLE)
 356    99 DBSIR4=S
 357 #endif
 358       RETURN
 359   101 FORMAT('INCORRECT ARGUMENT OR INDEX, X = ',1P,E15.6,5X,'NU = ',I5)
 360       END