Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/besk.ssp
There are 2 other files named besk.ssp in the archive. Click here to see a list.
C                                                                       BESK  10
C     ..................................................................BESK  20
C                                                                       BESK  30
C        SUBROUTINE BESK                                                BESK  40
C                                                                       BESK  50
C           COMPUTE THE K BESSEL FUNCTION FOR A GIVEN ARGUMENT AND ORDERBESK  60
C                                                                       BESK  70
C        USAGE                                                          BESK  80
C           CALL BESK(X,N,BK,IER)                                       BESK  90
C                                                                       BESK 100
C        DESCRIPTION OF PARAMETERS                                      BESK 110
C           X  -THE ARGUMENT OF THE K BESSEL FUNCTION DESIRED           BESK 120
C           N  -THE ORDER OF THE K BESSEL FUNCTION DESIRED              BESK 130
C           BK -THE RESULTANT K BESSEL FUNCTION                         BESK 140
C           IER-RESULTANT ERROR CODE WHERE                              BESK 150
C              IER=0  NO ERROR                                          BESK 160
C              IER=1  N IS NEGATIVE                                     BESK 170
C              IER=2  X IS ZERO OR NEGATIVE                             BESK 180
C              IER=3  X .GT. 170, MACHINE RANGE EXCEEDED                BESK 190
C              IER=4  BK .GT. 10**70                                    BESK 200
C                                                                       BESK 210
C        REMARKS                                                        BESK 220
C           N MUST BE GREATER THAN OR EQUAL TO ZERO                     BESK 230
C                                                                       BESK 240
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  BESK 250
C           NONE                                                        BESK 260
C                                                                       BESK 270
C        METHOD                                                         BESK 280
C           COMPUTES ZERO ORDER AND FIRST ORDER BESSEL FUNCTIONS USING  BESK 290
C           SERIES APPROXIMATIONS AND THEN COMPUTES N TH ORDER FUNCTION BESK 300
C           USING RECURRENCE RELATION.                                  BESK 310
C           RECURRENCE RELATION AND POLYNOMIAL APPROXIMATION TECHNIQUE  BESK 320
C           AS DESCRIBED BY A.J.M.HITCHCOCK,'POLYNOMIAL APPROXIMATIONS  BESK 330
C           TO BESSEL FUNCTIONS OF ORDER ZERO AND ONE AND TO RELATED    BESK 340
C           FUNCTIONS', M.T.A.C., V.11,1957,PP.86-88, AND G.N. WATSON,  BESK 350
C           'A TREATISE ON THE THEORY OF BESSEL FUNCTIONS', CAMBRIDGE   BESK 360
C           UNIVERSITY PRESS, 1958, P. 62                               BESK 370
C                                                                       BESK 380
C     ..................................................................BESK 390
C                                                                       BESK 400
      SUBROUTINE BESK(X,N,BK,IER)                                       BESK 410
      DIMENSION T(12)                                                   BESK 420
      BK=.0                                                             BESK 430
      IF(N)10,11,11                                                     BESK 440
   10 IER=1                                                             BESK 450
      RETURN                                                            BESK 460
   11 IF(X)12,12,20                                                     BESK 470
   12 IER=2                                                             BESK 480
      RETURN                                                            BESK 490
   20 IF(X-170.0)22,22,21                                               BESK 500
   21 IER=3                                                             BESK 510
      RETURN                                                            BESK 520
   22 IER=0                                                             BESK 530
      IF(X-1.)36,36,25                                                  BESK 540
   25 A=EXP(-X)                                                         BESK 550
      B=1./X                                                            BESK 560
      C=SQRT(B)                                                         BESK 570
      T(1)=B                                                            BESK 580
      DO 26 L=2,12                                                      BESK 590
   26 T(L)=T(L-1)*B                                                     BESK 600
      IF(N-1)27,29,27                                                   BESK 610
C                                                                       BESK 620
C     COMPUTE KO USING POLYNOMIAL APPROXIMATION                         BESK 630
C                                                                       BESK 640
   27 G0=A*(1.2533141-.1566642*T(1)+.08811128*T(2)-.09139095*T(3)       BESK 650
     2+.1344596*T(4)-.2299850*T(5)+.3792410*T(6)-.5247277*T(7)          BESK 660
     3+.5575368*T(8)-.4262633*T(9)+.2184518*T(10)-.06680977*T(11)       BESK 670
     4+.009189383*T(12))*C                                              BESK 680
      IF(N)20,28,29                                                     BESK 690
   28 BK=G0                                                             BESK 700
      RETURN                                                            BESK 710
C                                                                       BESK 720
C     COMPUTE K1 USING POLYNOMIAL APPROXIMATION                         BESK 730
C                                                                       BESK 740
   29 G1=A*(1.2533141+.4699927*T(1)-.1468583*T(2)+.1280427*T(3)         BESK 750
     2-.1736432*T(4)+.2847618*T(5)-.4594342*T(6)+.6283381*T(7)          BESK 760
     3-.6632295*T(8)+.5050239*T(9)-.2581304*T(10)+.07880001*T(11)       BESK 770
     4-.01082418*T(12))*C                                               BESK 780
      IF(N-1)20,30,31                                                   BESK 790
   30 BK=G1                                                             BESK 800
      RETURN                                                            BESK 810
C                                                                       BESK 820
C     FROM KO,K1 COMPUTE KN USING RECURRENCE RELATION                   BESK 830
C                                                                       BESK 840
   31 DO 35 J=2,N                                                       BESK 850
      GJ=2.*(FLOAT(J)-1.)*G1/X+G0                                       BESK 860
      IF(GJ-1.7E33)33,33,32                                             BESK 870
   32 IER=4                                                             BESK 880
      GO TO 34                                                          BESK 890
   33 G0=G1                                                             BESK 900
   35 G1=GJ                                                             BESK 910
   34 BK=GJ                                                             BESK 920
      RETURN                                                            BESK 930
   36 B=X/2.                                                            BESK 940
      A=.5772157+ALOG(B)                                                BESK 950
      C=B*B                                                             BESK 960
      IF(N-1)37,43,37                                                   BESK 970
C                                                                       BESK 980
C     COMPUTE KO USING SERIES EXPANSION                                 BESK 990
C                                                                       BESK1000
   37 G0=-A                                                             BESK1010
      X2J=1.                                                            BESK1020
      FACT=1.                                                           BESK1030
      HJ=.0                                                             BESK1040
      DO 40 J=1,6                                                       BESK1050
      RJ=1./FLOAT(J)                                                    BESK1060
      X2J=X2J*C                                                         BESK1070
      FACT=FACT*RJ*RJ                                                   BESK1080
      HJ=HJ+RJ                                                          BESK1090
   40 G0=G0+X2J*FACT*(HJ-A)                                             BESK1100
      IF(N)43,42,43                                                     BESK1110
   42 BK=G0                                                             BESK1120
      RETURN                                                            BESK1130
C                                                                       BESK1140
C     COMPUTE K1 USING SERIES EXPANSION                                 BESK1150
C                                                                       BESK1160
   43 X2J=B                                                             BESK1170
      FACT=1.                                                           BESK1180
      HJ=1.                                                             BESK1190
      G1=1./X+X2J*(.5+A-HJ)                                             BESK1200
      DO 50 J=2,8                                                       BESK1210
      X2J=X2J*C                                                         BESK1220
      RJ=1./FLOAT(J)                                                    BESK1230
      FACT=FACT*RJ*RJ                                                   BESK1240
      HJ=HJ+RJ                                                          BESK1250
   50 G1=G1+X2J*FACT*(.5+(A-HJ)*FLOAT(J))                               BESK1260
      IF(N-1)31,52,31                                                   BESK1270
   52 BK=G1                                                             BESK1280
      RETURN                                                            BESK1290
      END                                                               BESK1300