MINICERN/mathlib/gen/c/besi064.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:00  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       FUNCTION  BESI0(X)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       FUNCTION DBESI0(X)
  15 #include "gen/imp64.inc"
  16 #endif
  17       LOGICAL LEX
  18       CHARACTER NAME0*(*),NAME1*(*),NAME0E*(*),NAME1E*(*)
  19       CHARACTER*80 ERRTXT
  20       DIMENSION CI(0:24,0:1),CK(0:16,0:1)
  21
  22 #if !defined(CERNLIB_DOUBLE)
  23       PARAMETER (NAME0 = 'BESK0', NAME0E = 'EBESK0')
  24       PARAMETER (NAME1 = 'BESK1', NAME1E = 'EBESK1')
  25 #endif
  26 #if defined(CERNLIB_DOUBLE)
  27       PARAMETER (NAME0 = 'BESK0/DBESK0', NAME0E = 'EBESK0/DBESK0')
  28       PARAMETER (NAME1 = 'BESK1/DBESK1', NAME1E = 'EBESK1/DBESK1')
  29 #endif
  30       PARAMETER (EPS=1D-15)
  31       PARAMETER (Z1 = 1, HF = Z1/2)
  32       PARAMETER (PI = 3.14159 26535 89793 24D0)
  33       PARAMETER (CE = 0.57721 56649 01532 86D0)
  34       PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI2 = 1/(2*PI))
  35
  36       DATA CI( 0,0) /+1.00827 92054 58740 032D0/
  37       DATA CI( 1,0) /+0.00844 51226 24920 943D0/
  38       DATA CI( 2,0) /+0.00017 27006 30777 567D0/
  39       DATA CI( 3,0) /+0.00000 72475 91099 959D0/
  40       DATA CI( 4,0) /+0.00000 05135 87726 878D0/
  41       DATA CI( 5,0) /+0.00000 00568 16965 808D0/
  42       DATA CI( 6,0) /+0.00000 00085 13091 223D0/
  43       DATA CI( 7,0) /+0.00000 00012 38425 364D0/
  44       DATA CI( 8,0) /+0.00000 00000 29801 672D0/
  45       DATA CI( 9,0) /-0.00000 00000 78956 698D0/
  46       DATA CI(10,0) /-0.00000 00000 33127 128D0/
  47       DATA CI(11,0) /-0.00000 00000 04497 339D0/
  48       DATA CI(12,0) /+0.00000 00000 01799 790D0/
  49       DATA CI(13,0) /+0.00000 00000 00965 748D0/
  50       DATA CI(14,0) /+0.00000 00000 00038 604D0/
  51       DATA CI(15,0) /-0.00000 00000 00104 039D0/
  52       DATA CI(16,0) /-0.00000 00000 00023 950D0/
  53       DATA CI(17,0) /+0.00000 00000 00009 554D0/
  54       DATA CI(18,0) /+0.00000 00000 00004 443D0/
  55       DATA CI(19,0) /-0.00000 00000 00000 859D0/
  56       DATA CI(20,0) /-0.00000 00000 00000 709D0/
  57       DATA CI(21,0) /+0.00000 00000 00000 087D0/
  58       DATA CI(22,0) /+0.00000 00000 00000 112D0/
  59       DATA CI(23,0) /-0.00000 00000 00000 012D0/
  60       DATA CI(24,0) /-0.00000 00000 00000 018D0/
  61
  62       DATA CI( 0,1) /+0.97580 06023 26285 926D0/
  63       DATA CI( 1,1) /-0.02446 74429 63276 385D0/
  64       DATA CI( 2,1) /-0.00027 72053 60763 829D0/
  65       DATA CI( 3,1) /-0.00000 97321 46728 020D0/
  66       DATA CI( 4,1) /-0.00000 06297 24238 640D0/
  67       DATA CI( 5,1) /-0.00000 00659 61142 154D0/
  68       DATA CI( 6,1) /-0.00000 00096 13872 919D0/
  69       DATA CI( 7,1) /-0.00000 00014 01140 901D0/
  70       DATA CI( 8,1) /-0.00000 00000 47563 167D0/
  71       DATA CI( 9,1) /+0.00000 00000 81530 681D0/
  72       DATA CI(10,1) /+0.00000 00000 35408 148D0/
  73       DATA CI(11,1) /+0.00000 00000 05102 564D0/
  74       DATA CI(12,1) /-0.00000 00000 01804 409D0/
  75       DATA CI(13,1) /-0.00000 00000 01023 594D0/
  76       DATA CI(14,1) /-0.00000 00000 00052 678D0/
  77       DATA CI(15,1) /+0.00000 00000 00107 094D0/
  78       DATA CI(16,1) /+0.00000 00000 00026 120D0/
  79       DATA CI(17,1) /-0.00000 00000 00009 561D0/
  80       DATA CI(18,1) /-0.00000 00000 00004 713D0/
  81       DATA CI(19,1) /+0.00000 00000 00000 829D0/
  82       DATA CI(20,1) /+0.00000 00000 00000 743D0/
  83       DATA CI(21,1) /-0.00000 00000 00000 080D0/
  84       DATA CI(22,1) /-0.00000 00000 00000 117D0/
  85       DATA CI(23,1) /+0.00000 00000 00000 011D0/
  86       DATA CI(24,1) /+0.00000 00000 00000 019D0/
  87
  88       DATA CK( 0,0) /+0.98840 81742 30825 800D0/
  89       DATA CK( 1,0) /-0.01131 05046 46928 281D0/
  90       DATA CK( 2,0) /+0.00026 95326 12762 724D0/
  91       DATA CK( 3,0) /-0.00001 11066 85196 665D0/
  92       DATA CK( 4,0) /+0.00000 06325 75108 500D0/
  93       DATA CK( 5,0) /-0.00000 00450 47337 641D0/
  94       DATA CK( 6,0) /+0.00000 00037 92996 456D0/
  95       DATA CK( 7,0) /-0.00000 00003 64547 179D0/
  96       DATA CK( 8,0) /+0.00000 00000 39043 756D0/
  97       DATA CK( 9,0) /-0.00000 00000 04579 936D0/
  98       DATA CK(10,0) /+0.00000 00000 00580 811D0/
  99       DATA CK(11,0) /-0.00000 00000 00078 832D0/
 100       DATA CK(12,0) /+0.00000 00000 00011 360D0/
 101       DATA CK(13,0) /-0.00000 00000 00001 727D0/
 102       DATA CK(14,0) /+0.00000 00000 00000 275D0/
 103       DATA CK(15,0) /-0.00000 00000 00000 046D0/
 104       DATA CK(16,0) /+0.00000 00000 00000 008D0/
 105
 106       DATA CK( 0,1) /+1.03595 08587 72358 331D0/
 107       DATA CK( 1,1) /+0.03546 52912 43331 114D0/
 108       DATA CK( 2,1) /-0.00046 84750 28166 889D0/
 109       DATA CK( 3,1) /+0.00001 61850 63810 053D0/
 110       DATA CK( 4,1) /-0.00000 08451 72048 124D0/
 111       DATA CK( 5,1) /+0.00000 00571 32218 103D0/
 112       DATA CK( 6,1) /-0.00000 00046 45554 607D0/
 113       DATA CK( 7,1) /+0.00000 00004 35417 339D0/
 114       DATA CK( 8,1) /-0.00000 00000 45757 297D0/
 115       DATA CK( 9,1) /+0.00000 00000 05288 133D0/
 116       DATA CK(10,1) /-0.00000 00000 00662 613D0/
 117       DATA CK(11,1) /+0.00000 00000 00089 048D0/
 118       DATA CK(12,1) /-0.00000 00000 00012 726D0/
 119       DATA CK(13,1) /+0.00000 00000 00001 921D0/
 120       DATA CK(14,1) /-0.00000 00000 00000 305D0/
 121       DATA CK(15,1) /+0.00000 00000 00000 050D0/
 122       DATA CK(16,1) /-0.00000 00000 00000 009D0/
 123
 124       NU=0
 125       LEX=.FALSE.
 126       GO TO 6
 127
 128 #if !defined(CERNLIB_DOUBLE)
 129       ENTRY  EBESI0(X)
 130 #endif
 131 #if defined(CERNLIB_DOUBLE)
 132       ENTRY DEBSI0(X)
 133 #endif
 134       NU=0
 135       LEX=.TRUE.
 136       GO TO 6
 137
 138 #if !defined(CERNLIB_DOUBLE)
 139       ENTRY  BESI1(X)
 140 #endif
 141 #if defined(CERNLIB_DOUBLE)
 142       ENTRY DBESI1(X)
 143 #endif
 144       NU=1
 145       LEX=.FALSE.
 146       GO TO 6
 147
 148 #if !defined(CERNLIB_DOUBLE)
 149       ENTRY  EBESI1(X)
 150 #endif
 151 #if defined(CERNLIB_DOUBLE)
 152       ENTRY DEBSI1(X)
 153 #endif
 154       NU=1
 155       LEX=.TRUE.
 156
 157     6 V=ABS(X)
 158       IF(V .LT. 8) THEN
 159        Y=(HF*V)**2
 160        XL=NU+2
 161        A0=1
 162        A1=1+2*Y/((XL+1)*(XL-1))
 163        A2=1+Y*(4+3*Y/((XL+2)*XL))/((XL+3)*(XL-1))
 164        B0=1
 165        B1=1-Y/(XL+1)
 166        B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
 167        W1=3+XL
 168        V1=3-XL
 169        V3=XL-1
 170        V2=V3+V3
 171        C=0
 172        DO 3 N = 3,30
 173        C0=C
 174        FN=N
 175        W1=W1+2
 176        W2=W1-1
 177        W3=W2-1
 178        W4=W3-1
 179        W5=W4-1
 180        W6=W5-1
 181        V1=V1+1
 182        V2=V2+1
 183        V3=V3+1
 184        U1=FN*W4
 185        E=V3/(U1*W3)
 186        U2=E*Y
 187        F1=1+Y*V1/(U1*W1)
 188        F2=(1+Y*V2/(V3*W2*W5))*U2
 189        F3=-Y*Y*U2/(W4*W5*W5*W6)
 190        A=F1*A2+F2*A1+F3*A0
 191        B=F1*B2+F2*B1+F3*B0
 192        C=A/B
 193        IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 4
 194        A0=A1
 195        A1=A2
 196        A2=A
 197        B0=B1
 198        B1=B2
 199        B2=B
 200     3  CONTINUE
 201     4  H=C
 202        IF(NU .EQ. 1) H=HF*X*H
 203        IF(LEX) H=EXP(-V)*H
 204       ELSE
 205        R=1/V
 206        H=16*R-1
 207        ALFA=H+H
 208        B1=0
 209        B2=0
 210        DO 1 I = 24,0,-1
 211        B0=CI(I,NU)+ALFA*B1-B2
 212        B2=B1
 213     1  B1=B0
 214        H=SQRT(RPI2*R)*(B0-H*B2)
 215        IF(NU*X .LT. 0) H=-H
 216        IF(.NOT.LEX) H=EXP(V)*H
 217       ENDIF
 218       GO TO 9
 219
 220 #if !defined(CERNLIB_DOUBLE)
 221       ENTRY  BESK0(X)
 222 #endif
 223 #if defined(CERNLIB_DOUBLE)
 224       ENTRY DBESK0(X)
 225 #endif
 226       NU=0
 227       LEX=.FALSE.
 228       GO TO 8
 229
 230 #if !defined(CERNLIB_DOUBLE)
 231       ENTRY  EBESK0(X)
 232 #endif
 233 #if defined(CERNLIB_DOUBLE)
 234       ENTRY DEBSK0(X)
 235 #endif
 236       NU=0
 237       LEX=.TRUE.
 238       GO TO 8
 239
 240 #if !defined(CERNLIB_DOUBLE)
 241       ENTRY  BESK1(X)
 242 #endif
 243 #if defined(CERNLIB_DOUBLE)
 244       ENTRY DBESK1(X)
 245 #endif
 246       NU=1
 247       LEX=.FALSE.
 248       GO TO 8
 249
 250 #if !defined(CERNLIB_DOUBLE)
 251       ENTRY  EBESK1(X)
 252 #endif
 253 #if defined(CERNLIB_DOUBLE)
 254       ENTRY DEBSK1(X)
 255 #endif
 256       NU=1
 257       LEX=.TRUE.
 258
 259     8 IF(X .LE. 0) THEN
 260        H=0
 261        WRITE(ERRTXT,101) X
 262        IF(NU .EQ. 0 .AND. .NOT.LEX) CALL MTLPRT(NAME0 ,'C313.1',ERRTXT)
 263        IF(NU .EQ. 0 .AND.      LEX) CALL MTLPRT(NAME0E,'C313.1',ERRTXT)
 264        IF(NU .EQ. 1 .AND. .NOT.LEX) CALL MTLPRT(NAME1 ,'C313.1',ERRTXT)
 265        IF(NU .EQ. 1 .AND.      LEX) CALL MTLPRT(NAME1E,'C313.1',ERRTXT)
 266       ELSEIF(X .LT. 1) THEN
 267        B=HF*X
 268        BK=-(LOG(B)+CE)
 269        F=BK
 270        P=HF
 271        Q=HF
 272        C=1
 273        D=B**2
 274        BK1=P
 275        DO 11 N = 1,15
 276        FN=N
 277        RFN=1/FN
 278        P=P*RFN
 279        Q=Q*RFN
 280        F=(F+P+Q)*RFN
 281        C=C*D*RFN
 282        G=C*(P-FN*F)
 283        H=C*F
 284        BK=BK+H
 285        BK1=BK1+G
 286        IF(BK1*H+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
 287    11  CONTINUE
 288    12  H=BK
 289        IF(NU .EQ. 1) H=BK1/B
 290        IF(LEX) H=EXP(X)*H
 291       ELSEIF(X .LE. 5) THEN
 292        XN=4*NU**2
 293        A=9-XN
 294        B=25-XN
 295        C=768*X**2
 296        C0=48*X
 297        A0=1
 298        A1=(16*X+7+XN)/A
 299        A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
 300        B0=1
 301        B1=(16*X+9-XN)/A
 302        B2=(C+C0*B)/(A*B)+1
 303        C=0
 304        DO 24 N = 3,30
 305        C0=C
 306        FN=N
 307        FN2=FN+FN
 308        FN1=FN2-1
 309        FN3=FN1/(FN2-3)
 310        FN4=12*FN**2-(1-XN)
 311        FN5=16*FN1*X
 312        RAN=1/((FN2+1)**2-XN)
 313        F1=FN3*(FN4-20*FN)+FN5
 314        F2=28*FN-FN4-8+FN5
 315        F3=FN3*((FN2-5)**2-XN)
 316        A=(F1*A2+F2*A1+F3*A0)*RAN
 317        B=(F1*B2+F2*B1+F3*B0)*RAN
 318        C=A/B
 319        IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
 320        A0=A1
 321        A1=A2
 322        A2=A
 323        B0=B1
 324        B1=B2
 325        B2=B
 326    24  CONTINUE
 327    25  H=C/SQRT(RPIH*X)
 328        IF(.NOT.LEX) H=EXP(-X)*H
 329       ELSE
 330        R=1/X
 331        H=10*R-1
 332        ALFA=H+H
 333        B1=0
 334        B2=0
 335        DO 23 I = 16,0,-1
 336        B0=CK(I,NU)+ALFA*B1-B2
 337        B2=B1
 338    23  B1=B0
 339        H=SQRT(PIH*R)*(B0-H*B2)
 340        IF(.NOT.LEX) H=EXP(-X)*H
 341       ENDIF
 342     9 CONTINUE
 343 #if !defined(CERNLIB_DOUBLE)
 344        BESI0=H
 345 #endif
 346 #if defined(CERNLIB_DOUBLE)
 347       DBESI0=H
 348 #endif
 349       RETURN
 350   101 FORMAT(' NON-POSITIVE ARGUMENT X = ',1P,E15.6)
 351       END