@HDG,P    BESIS
@FOR,IS   BESIS
      SUBROUTINE BESIS (*,X,V,N1,N2,A)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     COMPUTE EXP(-X)*I(NU,X) WHERE I(NU,X) IS THE MODIFIED BESSEL
C     FUNCTION OF THE FIRST KIND FOR ORDER NU AND ARGUMENT X.
C     NU VARIES FROM V+MIN(N1,N2) THROUGH V+MAX(N1,N2).
C     STORE THE RESULT INTO A(1) THROUGH A(ABS(N1-N2)+1).
C
C     TAKE THE ABNORMAL EXIT IF V>=1 OR
C                               V<0 OR
C                              N1<0 OR
C                              N2<0 OR
C                               X<0.
C
      REAL A(1), V, X
      PARAMETER NN=6
      DATA PIM2/6.28318530717958647D0/
      INTEGER FTST(2)
      EQUIVALENCE (FMIN,FTST)
      DATA FTST/O160040000000,0/ @ SMALLEST S.P. NUMBER
      INTEGER PTST(2)
      EQUIVALENCE (PTST,P)
      DATA PTST/O174640000000,0/ @ 2**(-27)
C
      IF (V.GE.1) GO TO 90
      IF (V.LT.0) GO TO 90
      IF (X.LT.0) GO TO 90
      IF (N1.LT.0) GO TO 90
      IF (N2.LT.0) GO TO 90
C
C     I-TYPE BESSEL FUNCTIONS FOLLOW THE RECURRENCE RELATION
C     F(V-1,X)=(2*V/X)*F(V,X)+F(V+1,X).  TO AVOID OVERFLOW PROBLEMS,
C     LET R(V,X)=F(V+1,X)/F(V,X).  THEN R(V-1,X)=X/(2*V+X*R(V,X))
C     AND F(V-1,X)=F(V,X)/R(V-1,X)
C
C     COMPUTE I(V,X) AND I(V+1,X) FROM THE ASCENDING SEQUENCE
C     AND RECUR BACKWARD USING THE ABOVE RELATIONS.
C
      NMIN=MIN(N1,N2)
      NMAX=MAX(N1,N2)
      IF (X) ,60
      TWODX=2./X
      IF (NMAX.GT.30) GO TO 5
      IF (X.GT.MAX(NMAX*5,10)) GO TO 100
5     MU=NMAX
      VX=V+MU+1
      VXM=VX
      FM=EXP((1-VX)*LOG(TWODX)-DLGAMA(VX)-X)
      SUM=0
      ETA=0
      XSQ4=X*X*.25
      XN=1
      R=FM
      T=1
      IF ((FM-1)*(XSQ4-1).GT.0) GO TO 10
      R=1
      T=FM
10    SUM=SUM+T                        @ I(V,X) SEQUENCE
      T=T/VX
      ETA=ETA+T
      T=T*(XSQ4/XN)
      IF (T-P*SUM) 20
      XN=XN+1
      VX=VX+1
      GO TO 10
20    F=SUM*R                          @ I(V,X)
      R=ETA/(SUM*TWODX)                @ I(V+1,X) / I(V,X)
25    IF (F-FMIN) 30
      A(MU-NMIN+1)=F
      GO TO 40
30    A(MU-NMIN+1)=0
40    IF (MU.LE.NMIN) GO TO 80
      MU=MU-1
      VXM=VXM-1
      R=X/(VXM+VXM+X*R)
      F=F/R
      GO TO 25
C
C     SPECIAL FOR X = 0
C
60    N=NMAX-NMIN+1
      DO 70 I=1,N
70    A(I)=0
      IF (V) 80,,80
      IF (NMIN.EQ.0) A(1)=1.0          @ I(0,0)=1
80    RETURN
90    RETURN 1
C
C     IN THE REGION OF THE PLANE WHERE NMAX < 30 AND
C     X > MAX(NMAX*5,10), WE CAN APPROXIMATE I(V,X)
C     AND I(V+1,X) AND RECUR FORWARD.
C
100   R=4*V*V
      F=8*X
      VX=SQRT(PIM2*X)
      SUM=1
      DO 110 I=1,NN
      N=NN-I+1
      M=N+N-1
110   SUM=1.0-SUM*(R-M*M)/(F*N)
      FM=SUM/VX
      IF (NMAX.NE.0) GO TO 120
      A(1)=FM
      GO TO 80
120   R=4*(V+1)*(V+1)
      SUM=1
      DO 125 I=1,NN
      N=NN-I+1
      M=N+N-1
125   SUM=1.0-SUM*(R-M*M)/(F*N)
      F=SUM/VX
      IF (NMAX.GT.1) GO TO 140
      IF (NMAX.NE.NMIN) GO TO 160
130   A(1)=F
      GO TO 80
140   VX=TWODX*V
      N=NMAX-1
      DO 150 I=1,N
      VX=VX+TWODX
      R=FM-VX*F
      IF (I.GT.NMIN) A(I-NMIN)=FM
      FM=F
150   F=R
      IF (NMAX.EQ.NMIN) GO TO 130
160   A(NMAX-NMIN+1)=F
      A(NMAX-NMIN)=FM
      GO TO 80
      END
@HDG,P    DCI
@FOR,IS   DCI
      DOUBLE PRECISION FUNCTION DCI (X)
C
C     COMPUTE THE COSINE INTEGRAL OF X =
C     GAMMA + LOG(ABS(X)) + INTEGRAL FROM 0 TO X OF ((COS(T)-1)/T DT),
C     WHERE GAMMA IS EULER'S CONSTANT.
C
C     FOR ABS(X)<16, USE A CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE
C     Z=X/16 TO EVALUATE (CI(X)-LOG(ABS(X))-GAMMA)/(Z*Z), THEN MULTIPLY
C     THE RESULT BY Z*Z AND ADD LOG(ABS(X)) AND GAMMA.  LOSS OF ABOUT
C     TWO DIGITS OCCURS NEAR ABS(X)=16 DUE TO CANCELLATION.
C
C     FOR ABS(X).GE.16, USE CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE
C     WHERE Z=16/X ARE USED TO COMPUTE F(X)/Z AND G(X)/(Z*Z).  THEN
C     CI(X)=F(X)*SIN(X)-G(X)*COS(X).
C
C     THIS ALGORITHM YIELDS 15 DIGITS OF PRECISION.
C
      DOUBLE PRECISION X
      DOUBLE PRECISION GAMMA,Z,ZW,FZ,GZ,CHBPOL
      DATA GAMMA/0.57721 56649 01533 D0/
      DOUBLE PRECISION C(21)
      DATA C/
     1 +14.99258 93678 13409 D0, -19.38612 40966 07770 D0,
     2 +12.74187 08697 58071 D0, - 8.10790 39705 62531 D0,
     3 + 4.86202 23485 00627 D0, - 2.49750 50885 39025 D0,
     4 + 1.00866 07873 58110 D0, - 0.31208 09248 25428 D0,
     5 + 0.07467 82552 94576 D0, - 0.01411 08652 53535 D0,
     6 + 0.00215 20467 52074 D0, - 0.00027 02123 31184 D0,
     7 + 0.00002 84169 45498 D0, - 0.00000 25401 25611 D0,
     8 + 0.00000 01954 37144 D0, - 0.00000 00130 84020 D0,
     9 + 0.00000 00007 69379 D0, - 0.00000 00000 40066 D0,
     A + 0.00000 00000 01861 D0, - 0.00000 00000 00078 D0,
     B + 0.00000 00000 00003 D0/
      DOUBLE PRECISION F(11),G(11)
      DATA F/
     1 + 0.06226 37290 28927 D0, - 0.00023 37560 41393 D0,
     2 + 0.00000 24537 55677 D0, - 0.00000 00586 70317 D0,
     3 + 0.00000 00023 56196 D0, - 0.00000 00001 36096 D0,
     4 + 0.00000 00000 10308 D0, - 0.00000 00000 00964 D0,
     5 + 0.00000 00000 00107 D0, - 0.00000 00000 00014 D0,
     6 + 0.00000 00000 00002 D0/
      DATA G/
     1 + 0.00386 28560 96703 D0, - 0.00004 26441 82622 D0,
     2 + 0.00000 07249 95950 D0, - 0.00000 00234 68225 D0,
     3 + 0.00000 00011 69202 D0, - 0.00000 00000 79604 D0,
     4 + 0.00000 00000 06875 D0, - 0.00000 00000 00717 D0,
     5 + 0.00000 00000 00087 D0, - 0.00000 00000 00012 D0,
     6 + 0.00000 00000 00002 D0/
C
      IF (DABS(X).GE.16.0D0) GO TO 10
      Z=X/16.0D0
      ZW=Z*Z
      GZ=ZW+ZW-1.0D0
      Z=ZW*CHBPOL(C,21,GZ)
      Z=DLOG(DABS(X))-Z+GAMMA
      GO TO 20
C
10    Z=16.0D0/X
      ZW=Z*Z
      GZ=ZW+ZW-1.0D0
      FZ=Z*CHBPOL(F,11,GZ)
      GZ=ZW*CHBPOL(G,11,GZ)
      Z=FZ*DSIN(X)-GZ*DCOS(X)
C
20    DCI=Z
      RETURN
      END
@HDG,P    DSI
@FOR,IS   DSI
      DOUBLE PRECISION FUNCTION DSI (X)
C
C     COMPUTE THE SINE INTEGRAL OF X =
C     INTEGRAL FROM 0 TO X OF (SIN(T)/T DT).
C
C     FOR ABS(X)<16, USE A CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE
C     Z=X/16 TO EVALUATE SI(X)/Z, THEN MULTIPLY THE RESULT BY Z.  THIS
C     AVOIDS STORING ZERO COEFFICIENTS FOR EVEN ORDERS, AND PRESERVES
C     ACCURACY FOR SMALL Z.
C
C     FOR ABS(X).GE.16, USE CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1,
C     WHERE Z=16/X ARE USED TO CUMPUTE F(X)/Z AND G(X)/(Z*Z).  THEN
C     SI(X)=PI*SIGN(X)*0.5-F(X)*COS(X)-G(X)*SIN(X).
C
C     THIS ALGORITHM YIELDS 15 DIGITS OF PRECISION.
C
      DOUBLE PRECISION X
      DOUBLE PRECISION PI2,Z,ZW,FZ,GZ,CHBPOL
      DATA PI2/1.57079 63267 94896 62 D0/
      DOUBLE PRECISION S(21)
      DATA S/
     1 + 4.05292 64776 80623 D0, - 4.06398 08449 11986 D0,
     2 + 2.77875 63817 42663 D0, - 1.92656 50911 50656 D0,
     3 + 1.38930 87711 71888 D0, - 0.96832 22369 87086 D0,
     4 + 0.53014 88479 16522 D0, - 0.21126 37809 76555 D0,
     5 + 0.06203 36794 32003 D0, - 0.01386 74455 89417 D0,
     6 + 0.00243 62214 04749 D0, - 0.00034 54691 55569 D0,
     7 + 0.00004 04202 71419 D0, - 0.00000 39729 08746 D0,
     8 + 0.00000 03329 88589 D0, - 0.00000 00241 00076 D0,
     9 + 0.00000 00015 22370 D0, - 0.00000 00000 84710 D0,
     A + 0.00000 00000 04185 D0, - 0.00000 00000 00185 D0,
     B + 0.00000 00000 00007 D0/
      DOUBLE PRECISION F(11),G(11)
      DATA F/
     1 + 0.06226 37290 28927 D0, - 0.00023 37560 41393 D0,
     2 + 0.00000 24537 55677 D0, - 0.00000 00586 70317 D0,
     3 + 0.00000 00023 56196 D0, - 0.00000 00001 36096 D0,
     4 + 0.00000 00000 10308 D0, - 0.00000 00000 00964 D0,
     5 + 0.00000 00000 00107 D0, - 0.00000 00000 00014 D0,
     6 + 0.00000 00000 00002 D0/
      DATA G/
     1 + 0.00386 28560 96703 D0, - 0.00004 26441 82622 D0,
     2 + 0.00000 07249 95950 D0, - 0.00000 00234 68225 D0,
     3 + 0.00000 00011 69202 D0, - 0.00000 00000 79604 D0,
     4 + 0.00000 00000 06875 D0, - 0.00000 00000 00717 D0,
     5 + 0.00000 00000 00087 D0, - 0.00000 00000 00012 D0,
     6 + 0.00000 00000 00002 D0/
C
      IF (DABS(X).GE.16.0D0) GO TO 10
      Z=X/16.0D0
      ZW=2.0D0*Z*Z-1.0D0
      Z=Z*CHBPOL(S,21,ZW)
      GO TO 20
C
10    Z=16.0D0/X
      ZW=Z*Z
      GZ=ZW+ZW-1.0D0
      FZ=Z*CHBPOL(F,11,GZ)
      GZ=ZW*CHBPOL(G,11,GZ)
      Z=DSIGN(PI2,X)-FZ*DCOS(X)-GZ*DSIN(X)
C
20    DSI=Z
      RETURN
      END
@HDG,P    PSI
@FOR,IS   PSI
      FUNCTION PSI(Z)
C
C     COMPUTE THE FUNCTION  PSI(Z) = GAMMA'(Z)/GAMMA(Z)
C     RATIONAL APPROXIMATIONS DUE TO CODY ET AL ARE USED.
C
      DOUBLE PRECISION ZD,ZF,ZW,PSIW,PI,P(3),Q(3)
      DATA PI /3.14159265358979324D0/
      DATA P /-2.002880963995D-9, -1.528276172927D0, -1.399162842582D0/
      DATA Q /1.833932086804D1, 1.862345253239D1, 1.D0/
C
C     SEPARATE Z INTO FRACTIONAL AND INTEGER PARTS.
C
      IZ=Z
      ZD=Z
      ZF=ABS(ZD-IZ)
C
C     RETURN  -INFINITY FOR Z=0  OR  Z=ANY NEGATIVE INTEGER.
C
      IF (IZ.GT.0) GO TO 10
      IF (ABS(ZF).GT.1.7D-38) GO TO 10
      PSI=-1.7D38
      RETURN
C
C     COMPUTE PSI(Z).  IF  0 < Z < 3  COMPUTE PSI(3+ZF) AND RECUR
C     BACKWARD TO PSI(Z) USING  PSI(Z) = PSI(Z+1) - 1/Z.  IF  -2 < Z < 0
C     USE THE ABOVE RECURRENCE FORMULA STARTING WITH PSI(3+ZF) TO
C     COMPUTE PSI(1+ABS(Z)).  THEN IF Z IS NEGATIVE USE THE REFLECTION
C     FORMULA  PSI(-Z) = PSI(1+Z) + PI*COT(PI*ZF).
C
10    ZW=ABS(ZD)
      NREC=0
      IF (IZ.GT.2) GO TO 20
      IF (IZ.LT.-1) GO TO 20
      ZW=3+ZF
      NREC=3-ABS(IZ)
      IF (Z.LT.0) NREC=NREC-1
      GO TO 30
20    IF (Z.LT.0) ZW=1+ZW
30    PSIW=LOG(ZW)-1/(ZW+ZW)
      ZD=ZW*ZW
      PSIW=PSIW+(P(1)+(P(2)+P(3)/ZD)/ZD)/(Q(1)+(Q(2)+Q(3)/ZD)/ZD)
40    IF (NREC.EQ.0) GO TO 50
      ZW=ZW-1
      PSIW=PSIW-1/ZW
      NREC=NREC-1
      GO TO 40
50    IF (Z.GT.0) GO TO 60
      PSIW=PSIW+PI*DCOTAN(PI*ZF)
60    PSI=PSIW
      RETURN
      END
