MINICERN/mathlib/gen/c/bsir364.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:07  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       FUNCTION  BSIR3(X,NU)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       FUNCTION DBSIR3(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 = 'BSIR3', NAMEIE = 'EBSIR3')
  21       PARAMETER (NAMEK = 'BSKR3', NAMEKE = 'EBSKR3')
  22 #endif
  23 #if defined(CERNLIB_DOUBLE)
  24       PARAMETER (NAMEI = 'BSIR3/DBSIR3', NAMEIE = 'EBSIR3/DEBIR3')
  25       PARAMETER (NAMEK = 'BSKR3/DBSKR3', NAMEKE = 'EBSKR3/DEBKR3')
  26 #endif
  27       LOGICAL LEX
  28
  29       DIMENSION BC(0:23,2),CC(0:15,2),PP(-2:2),GG(-2:2)
  30
  31       PARAMETER (EPS = 1D-15)
  32       PARAMETER (Z1 = 1, HF = Z1/2)
  33       PARAMETER (PI = 3.14159 26535 89793 24D0)
  34       PARAMETER (W3 = 1.73205 08075 68877 29D0)
  35       PARAMETER (G1 = 2.67893 85347 07747 63D0)
  36       PARAMETER (G2 = 1.35411 79394 26400 42D0)
  37       PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI2 = 1/(2*PI))
  38       PARAMETER (C1 = 2*PI/(3*W3))
  39       PARAMETER (GM = 3*(1/G2-3/G1)/2, GP = (3/G1+1/G2)/2)
  40
  41       DATA PP(-2) /-0.66666 66666 66666 67D0/
  42       DATA PP(-1) /-0.33333 33333 33333 33D0/
  43       DATA PP( 1) / 0.33333 33333 33333 33D0/
  44       DATA PP( 2) / 0.66666 66666 66666 67D0/
  45
  46       DATA GG(-2) / 0.37328 21739 07395 23D0/
  47       DATA GG(-1) / 0.73848 81116 21648 31D0/
  48       DATA GG( 1) / 1.11984 65217 22185 68D0/
  49       DATA GG( 2) / 1.10773 21674 32472 47D0/
  50
  51       DATA BC( 0,1) / 1.00458 61710 93207 35D0/
  52       DATA BC( 1,1) / 0.00467 34791 99873 60D0/
  53       DATA BC( 2,1) / 0.00009 08034 04815 04D0/
  54       DATA BC( 3,1) / 0.00000 37262 16110 59D0/
  55       DATA BC( 4,1) / 0.00000 02520 73237 90D0/
  56       DATA BC( 5,1) / 0.00000 00227 82110 77D0/
  57       DATA BC( 6,1) / 0.00000 00012 91332 28D0/
  58       DATA BC( 7,1) /-0.00000 00006 11915 16D0/
  59       DATA BC( 8,1) /-0.00000 00003 75616 85D0/
  60       DATA BC( 9,1) /-0.00000 00001 16415 46D0/
  61       DATA BC(10,1) /-0.00000 00000 14443 25D0/
  62       DATA BC(11,1) / 0.00000 00000 05373 69D0/
  63       DATA BC(12,1) / 0.00000 00000 03074 27D0/
  64       DATA BC(13,1) / 0.00000 00000 00297 66D0/
  65       DATA BC(14,1) /-0.00000 00000 00265 20D0/
  66       DATA BC(15,1) /-0.00000 00000 00091 37D0/
  67       DATA BC(16,1) / 0.00000 00000 00015 52D0/
  68       DATA BC(17,1) / 0.00000 00000 00014 12D0/
  69       DATA BC(18,1) /-0.00000 00000 00000 23D0/
  70       DATA BC(19,1) /-0.00000 00000 00001 98D0/
  71       DATA BC(20,1) /-0.00000 00000 00000 13D0/
  72       DATA BC(21,1) / 0.00000 00000 00000 29D0/
  73       DATA BC(22,1) / 0.00000 00000 00000 03D0/
  74       DATA BC(23,1) /-0.00000 00000 00000 05D0/
  75
  76       DATA BC( 0,2) / 0.99363 49867 16925 14D0/
  77       DATA BC( 1,2) /-0.00646 71526 00616 03D0/
  78       DATA BC( 2,2) /-0.00010 60188 22351 55D0/
  79       DATA BC( 3,2) /-0.00000 41406 57716 24D0/
  80       DATA BC( 4,2) /-0.00000 02916 95418 21D0/
  81       DATA BC( 5,2) /-0.00000 00365 71574 33D0/
  82       DATA BC( 6,2) /-0.00000 00075 81590 37D0/
  83       DATA BC( 7,2) /-0.00000 00019 23008 52D0/
  84       DATA BC( 8,2) /-0.00000 00004 20438 80D0/
  85       DATA BC( 9,2) /-0.00000 00000 39372 04D0/
  86       DATA BC(10,2) / 0.00000 00000 19007 44D0/
  87       DATA BC(11,2) / 0.00000 00000 10137 64D0/
  88       DATA BC(12,2) / 0.00000 00000 01331 30D0/
  89       DATA BC(13,2) /-0.00000 00000 00676 92D0/
  90       DATA BC(14,2) /-0.00000 00000 00311 72D0/
  91       DATA BC(15,2) / 0.00000 00000 00011 87D0/
  92       DATA BC(16,2) / 0.00000 00000 00040 21D0/
  93       DATA BC(17,2) / 0.00000 00000 00004 78D0/
  94       DATA BC(18,2) /-0.00000 00000 00004 74D0/
  95       DATA BC(19,2) /-0.00000 00000 00001 16D0/
  96       DATA BC(20,2) / 0.00000 00000 00000 59D0/
  97       DATA BC(21,2) / 0.00000 00000 00000 21D0/
  98       DATA BC(22,2) /-0.00000 00000 00000 08D0/
  99       DATA BC(23,2) /-0.00000 00000 00000 03D0/
 100
 101       DATA CC( 0,1) / 0.99353 64122 76093 39D0/
 102       DATA CC( 1,1) /-0.00631 44392 60798 63D0/
 103       DATA CC( 2,1) / 0.00014 30095 80961 13D0/
 104       DATA CC( 3,1) /-0.00000 57870 60592 03D0/
 105       DATA CC( 4,1) / 0.00000 03265 50333 20D0/
 106       DATA CC( 5,1) /-0.00000 00231 23231 95D0/
 107       DATA CC( 6,1) / 0.00000 00019 39555 14D0/
 108       DATA CC( 7,1) /-0.00000 00001 85897 89D0/
 109       DATA CC( 8,1) / 0.00000 00000 19868 42D0/
 110       DATA CC( 9,1) /-0.00000 00000 02326 79D0/
 111       DATA CC(10,1) / 0.00000 00000 00294 68D0/
 112       DATA CC(11,1) /-0.00000 00000 00039 95D0/
 113       DATA CC(12,1) / 0.00000 00000 00005 75D0/
 114       DATA CC(13,1) /-0.00000 00000 00000 87D0/
 115       DATA CC(14,1) / 0.00000 00000 00000 14D0/
 116       DATA CC(15,1) /-0.00000 00000 00000 02D0/
 117
 118       DATA CC( 0,2) / 1.00914 95380 72789 40D0/
 119       DATA CC( 1,2) / 0.00897 12068 42483 60D0/
 120       DATA CC( 2,2) /-0.00017 13895 98261 54D0/
 121       DATA CC( 3,2) / 0.00000 65547 92549 82D0/
 122       DATA CC( 4,2) /-0.00000 03595 19190 49D0/
 123       DATA CC( 5,2) / 0.00000 00250 24412 19D0/
 124       DATA CC( 6,2) /-0.00000 00020 74924 13D0/
 125       DATA CC( 7,2) / 0.00000 00001 97223 67D0/
 126       DATA CC( 8,2) /-0.00000 00000 20946 47D0/
 127       DATA CC( 9,2) / 0.00000 00000 02440 93D0/
 128       DATA CC(10,2) /-0.00000 00000 00307 91D0/
 129       DATA CC(11,2) / 0.00000 00000 00041 61D0/
 130       DATA CC(12,2) /-0.00000 00000 00005 97D0/
 131       DATA CC(13,2) / 0.00000 00000 00000 91D0/
 132       DATA CC(14,2) /-0.00000 00000 00000 14D0/
 133       DATA CC(15,2) / 0.00000 00000 00000 02D0/
 134
 135       LEX=.FALSE.
 136       GO TO 8
 137
 138 #if !defined(CERNLIB_DOUBLE)
 139       ENTRY  EBSIR3(X,NU)
 140 #endif
 141 #if defined(CERNLIB_DOUBLE)
 142       ENTRY DEBIR3(X,NU)
 143 #endif
 144       LEX=.TRUE.
 145
 146     8 MU=ABS(NU)
 147       IF(MU .NE. 1 .AND. MU .NE. 2 .OR. NU .LT. 0 .AND. X .LE. 0
 148      1   .OR. NU .GT. 0 .AND. X .LT. 0) THEN
 149        S=0
 150        WRITE(ERRTXT,101) X,NU
 151        IF(.NOT.LEX) CALL MTLPRT(NAMEI ,'C340.1',ERRTXT)
 152        IF(     LEX) CALL MTLPRT(NAMEIE,'C340.1',ERRTXT)
 153       ELSEIF(X .EQ. 0) THEN
 154        S=0
 155       ELSEIF(X .LT. 8) THEN
 156        Y=(HF*X)**2
 157        XN=PP(NU)
 158        XL=XN+2
 159        A0=1
 160        A1=1+2*Y/((XL+1)*(XN+1))
 161        A2=1+Y*(4+3*Y/((XL+2)*(XN+2)))/((XL+3)*(XN+1))
 162        B0=1
 163        B1=1-Y/(XL+1)
 164        B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
 165        T1=3+XL
 166        V1=3-XL
 167        V3=XL-1
 168        V2=V3+V3
 169        C=0
 170        DO 33 N = 3,30
 171        C0=C
 172        T1=T1+2
 173        T2=T1-1
 174        T3=T2-1
 175        T4=T3-1
 176        T5=T4-1
 177        T6=T5-1
 178        V1=V1+1
 179        V2=V2+1
 180        V3=V3+1
 181        U1=N*T4
 182        E=V3/(U1*T3)
 183        U2=E*Y
 184        F1=1+Y*V1/(U1*T1)
 185        F2=(1+Y*V2/(V3*T2*T5))*U2
 186        F3=-Y*Y*U2/(T4*T5*T5*T6)
 187        A=F1*A2+F2*A1+F3*A0
 188        B=F1*B2+F2*B1+F3*B0
 189        C=A/B
 190        IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 34
 191        A0=A1
 192        A1=A2
 193        A2=A
 194        B0=B1
 195        B1=B2
 196        B2=B
 197    33  CONTINUE
 198    34  S=GG(NU)*(HF*X)**PP(NU)*C
 199        IF(LEX) S=EXP(-X)*S
 200       ELSE
 201        R=1/X
 202        W=SQRT(RPI2*R)
 203        H=16*R-1
 204        ALFA=H+H
 205        B1=0
 206        B2=0
 207        DO 2 I = 23,0,-1
 208        B0=BC(I,MU)+ALFA*B1-B2
 209        B2=B1
 210     2  B1=B0
 211        S=W*(B0-H*B2)
 212        IF(.NOT.LEX) S=EXP(X)*S
 213        T=0
 214        IF(NU .LT. 0) THEN
 215         H=10*R-1
 216         ALFA=H+H
 217         B1=0
 218         B2=0
 219         DO 3 I = 15,0,-1
 220         B0=CC(I,MU)+ALFA*B1-B2
 221         B2=B1
 222     3   B1=B0
 223         R=EXP(-X)
 224         T=W3*W*R*(B0-H*B2)
 225         IF(LEX) T=R*T
 226        END IF
 227        S=S+T
 228       END IF
 229       GO TO 99
 230
 231 #if !defined(CERNLIB_DOUBLE)
 232       ENTRY  BSKR3(X,NU)
 233 #endif
 234 #if defined(CERNLIB_DOUBLE)
 235       ENTRY DBSKR3(X,NU)
 236 #endif
 237       LEX=.FALSE.
 238       GO TO 9
 239
 240 #if !defined(CERNLIB_DOUBLE)
 241       ENTRY  EBSKR3(X,NU)
 242 #endif
 243 #if defined(CERNLIB_DOUBLE)
 244       ENTRY DEBKR3(X,NU)
 245 #endif
 246       LEX=.TRUE.
 247
 248     9 MU=ABS(NU)
 249       IF(MU .NE. 1 .AND. MU .NE. 2 .OR. X .LE. 0) THEN
 250        S=0
 251        WRITE(ERRTXT,101) X,NU
 252        IF(.NOT.LEX) CALL MTLPRT(NAMEK ,'C340.1',ERRTXT)
 253        IF(     LEX) CALL MTLPRT(NAMEKE,'C340.1',ERRTXT)
 254       ELSEIF(X .LE. 1) THEN
 255        A0=PP(-1)
 256        B=HF*X
 257        D=-LOG(B)
 258        F=A0*D
 259        E=EXP(F)
 260        G=(GM*A0+GP)*E
 261        BK=C1*(HF*GM*(E+1/E)+GP*D*SINH(F)/F)
 262        F=BK
 263        E=A0**2
 264        P=HF*C1*G
 265        Q=HF/G
 266        C=1
 267        D=B**2
 268        BK1=P
 269        DO 11 N = 1,15
 270        FN=N
 271        F=(FN*F+P+Q)/(FN**2-E)
 272        C=C*D/FN
 273        P=P/(FN-A0)
 274        Q=Q/(FN+A0)
 275        G=C*(P-FN*F)
 276        H=C*F
 277        BK=BK+H
 278        BK1=BK1+G
 279        IF(H*BK1+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
 280    11  CONTINUE
 281    12  S=BK
 282        IF(MU .EQ. 2) S=BK1/B
 283        IF(LEX) S=EXP(X)*S
 284       ELSEIF(X .LE. 5) THEN
 285        XN=4*PP(MU)**2
 286        A=9-XN
 287        B=25-XN
 288        C=768*X**2
 289        C0=48*X
 290        A0=1
 291        A1=(16*X+7+XN)/A
 292        A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
 293        B0=1
 294        B1=(16*X+9-XN)/A
 295        B2=(C+C0*B)/(A*B)+1
 296        C=0
 297        DO 24 N = 3,30
 298        C0=C
 299        FN=N
 300        FN2=FN+FN
 301        FN1=FN2-1
 302        FN3=FN1/(FN2-3)
 303        FN4=12*FN**2-(1-XN)
 304        FN5=16*FN1*X
 305        RAN=1/((FN2+1)**2-XN)
 306        F1=FN3*(FN4-20*FN)+FN5
 307        F2=28*FN-FN4-8+FN5
 308        F3=FN3*((FN2-5)**2-XN)
 309        A=(F1*A2+F2*A1+F3*A0)*RAN
 310        B=(F1*B2+F2*B1+F3*B0)*RAN
 311        C=A/B
 312        IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
 313        A0=A1
 314        A1=A2
 315        A2=A
 316        B0=B1
 317        B1=B2
 318        B2=B
 319    24  CONTINUE
 320    25  S=C/SQRT(RPIH*X)
 321        IF(.NOT.LEX) S=EXP(-X)*S
 322       ELSE
 323        R=1/X
 324        H=10*R-1
 325        ALFA=H+H
 326        B1=0
 327        B2=0
 328        DO 13 I = 15,0,-1
 329        B0=CC(I,MU)+ALFA*B1-B2
 330        B2=B1
 331    13  B1=B0
 332        S=SQRT(PIH*R)*(B0-H*B2)
 333        IF(.NOT.LEX) S=EXP(-X)*S
 334       END IF
 335 #if !defined(CERNLIB_DOUBLE)
 336    99  BSIR3=S
 337 #endif
 338 #if defined(CERNLIB_DOUBLE)
 339    99 DBSIR3=S
 340 #endif
 341       RETURN
 342   101 FORMAT('INCORRECT ARGUMENT OR INDEX, X = ',1P,E15.6,' NU = ',I5)
 343       END