Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-04 - decus/20-0111/beslib.for
There is 1 other file named beslib.for in the archive. Click here to see a list.
      FUNCTION DJNP(N,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-R,T-Z)
      COMMON/CIEXP/IEXP,FEXPC
      REAL FLOAT
      M=IABS(N)
      DM=DBLE(FLOAT(M))
      IF(X-0.0D0)11,8,11
    8 IEXP=0
      IF(M-1)10,9,10
    9 DJNP=0.5D0*DBLE(FLOAT(N))
      GO TO 39
   10 DJNP=0.0D0
      GO TO 39
   11 IF(M)13,12,13
   12 DJNP=-DJN(1,X)
      GO TO 35
   13 DJX=DJN(M,X)
      IJX=IEXP
      DJTX=DJN(M-1,X)
      IJTX=IEXP
      DJNP=QFADD((-(DM*DJX)/X),DJTX,IJX,IJTX)
   35 MEO=M-2*(M/2)
      IF((N.LT.0).AND.(MEO.GT.0)) DJNP=-DJNP
   39 RETURN
      END
      FUNCTION DINP(N,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-R,T-Z)
      COMMON/CIEXP/IEXP,FEXPC
      REAL FLOAT
      M=IABS(N)
      DM=DBLE(FLOAT(M))
      IF(X-0.0D0)11,8,11
    8 IEXP=0
      IF(M-1)10,9,10
    9 DINP=0.5D0*DBLE(FLOAT(N))
      GO TO 35
   10 DINP=0.0D0
      GO TO 35
   11 IF(M)13,12,13
   12 DINP=DIN(1,X)
      GO TO 35
   13 DIX=DIN(M,X)
      IIX=IEXP
      DITX=DIN(M+1,X)
      IITX=IEXP
      DINP=QFADD(((DM*DIX)/X),DITX,IIX,IITX)
   35 RETURN
      END
      FUNCTION DNNP(N,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-R,T-Z)
      COMMON/CIEXP/IEXP,FEXPC
      REAL FLOAT
      M=IABS(N)
      IF(X-0.0D0)11,10,11
   10 DNNP=1.0D10
      IEXP=80000
      GO TO 35
   11 IF(M)13,12,13
   12 DNNP=-DNN(1,X)
      GO TO 35
   13 DM=DBLE(FLOAT(M))
      DNX=DNN(M,X)
      INX=IEXP
      DNTX=DNN(M-1,X)
      INTX=IEXP
      DNNP=QFADD((-(DM*DNX)/X),DNTX,INX,INTX)
   35 MEO=M-2*(M/2)
      IF((N.LT.0).AND.(MEO.GT.0)) DNNP=-DNNP
   39 RETURN
      END
      FUNCTION DKNP(N,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-R,T-Z)
      COMMON/CIEXP/IEXP,FEXPC
      REAL FLOAT
      M=IABS(N)
      IF(X-0.0D0)11,10,11
   10 DKNP=-1.0D10
      IEXP=80000
      GO TO 35
   11 IF(M)13,12,13
   12 DKNP=-DKN(1,X)
      GO TO 35
   13 DM=DBLE(FLOAT(M))
      DKX=DKN(M,X)
      IKX=IEXP
      DKTX=DKN(M-1,X)
      IKTX=IEXP
      DKNP=-QFADD(((DM*DKX)/X),DKTX,IKX,IKTX)
   35 RETURN
      END
      FUNCTION DJN(N,X)
      DOUBLE PRECISION DJN,X,DTN,FIEXP,FEXPC,DJR,DJR1,DJR2,FCR,DBLE,FNX,
     1FNX22,DFN,P0,P0T,XTJ,DJNF,XA,DABS
      COMMON/CIEXP/IEXP,FEXPC
      XA=DABS(X)
      IXS=1
      M=IABS(N)
      IF(X-0.0D0)9,10,11
   10 IF(M-1)12,13,13
   12 DJN=1.0D0
      IEXP=0
      GO TO 35
   13 DJN=0.0D0
      IEXP=0
      GO TO 35
    9 IXS=-1
   11 IF(XA-XTJ)8,5,8
    5 IF(M-MTJ)8,6,8
    8 IF(M-5)14,15,15
   14 IF(XA-2.0D0)25,25,16
   25 DJN=DTN(1,0,M,XA)
      P0T=P0(M,XA)
      IES=IEXP
      DJN=FIEXP(DJN*P0T)
      IEXP=IEXP+IES
      GO TO 35
   16 IF(XA-17.0D0)26,30,30
   26 DJR1=DTN(1,0,40,XA)
      IES=0
      DJR2=DTN(1,0,41,XA)
      FCR=1.0D0
      NSMAX=40
   27 FNX22=(XA**2)/4.0D0
      DO 9999 NS=NSMAX,M+1,-1
      FLNS=FLOAT(NS)
      DJR=FIEXP(DJR1-((DJR2*FNX22)/(FLNS*(FLNS+1.0D0)))*(FCR))
      FCR=FEXPC
      IES=IES+IEXP
      DJR2=DJR1
      DJR1=DJR
 9999 CONTINUE
      P0T=P0(M,XA)
      IESS=IEXP
      DJN=FIEXP(DJR*P0T)
      IEXP=IEXP+IES+IESS
      GO TO 35
   30 DJN=DTN(1,1,M,XA)
      IEXP=0
      GO TO 35
   15 FN=FLOAT(M)
      DFN=DBLE(FN)
      IF(M-15)17,17,21
   17 IF(XA-DFN/2.0D0)25,25,18
   18 IF(XA-17.0D0)26,19,19
   19 IF(M-13)30,30,51
   51 IF(XA-21.0D0)29,30,30
   21 X2=XA*XA
      LMAX=INT(X2/(2.0E0*(FN+SQRT(FN*FN+X2))))
      IF(LMAX-1)25,25,22
   22 IF(XA-DFN)28,52,52
   52 IF(M-18)53,53,54
   53 IF(XA-25.0D0)29,30,30
   54 FN2=FLOAT(M*M)
      XS=XA
      FKIN=(FN2-0.25E0)/(XS-0.5E0+SQRT(XS*(XS-1.0E0)+FN2))
      IF(FKIN-6.3E0)30,30,29
   29 DJR1=DTN(1,1,1,XA)
      DJR2=DTN(1,1,2,XA)
      IER=0
      FCR=1.0D0
      DO 9998 NS=2,M-1
      DJR=FIEXP(((2.0D0*NS*DJR2)/XA)-DJR1*FCR)
      IER=IEXP+IER
      FCR=FEXPC
      DJR1=DJR2
      DJR2=DJR
 9998 CONTINUE
      DJN=DJR
      IEXP=IER
      GO TO 35
   28 DJR1=DTN(1,0,100,XA)
      IES=0
      DJR2=DTN(1,0,101,XA)
      FCR=1.0D0
      NSMAX=100
      GO TO 27
   35 DJNF=DJN
      IJNF=IEXP
      XTJ=XA
      MTJ=M
      GO TO 36
    6 DJN=DJNF
      IEXP=IJNF
   36 MEO=M-2*(M/2)
      IF((N.LT.0).AND.(MEO.GT.0)) DJN=-DJN
      IF((IXS.LT.0).AND.(MEO.GT.0)) DJN=-DJN
   39 RETURN
      END
      FUNCTION DIN(N,X)
      DOUBLE PRECISION DIN,X,FEXPC,FIEXP,DTN,P0,QEXP,QX,DFN,DBLE,P0T,XTI
     1,XA,DABS,DINF
      COMMON/CIEXP/IEXP,FEXPC
      XA=DABS(X)
      IXS=1
      M=IABS(N)
      IF(X-0.0D0)9,10,11
   10 IF(M-1)12,13,13
   12 DIN=1.0D0
      IEXP=0
      GO TO 35
   13 DIN=0.0D0
      IEXP=0
      GO TO 35
    9 IXS=-1
   11 IF(XA-XTI)8,5,8
    5 IF(M-MTI)8,6,8
    8 IF(M-5)14,15,15
   14 IF(XA-17.0D0)40,45,45
   15 DFN=DBLE(FLOAT(M))
      IF(XA-5.0D0*DFN)40,45,45
   40 DIN=DTN(4,0,M,XA)
      IES=IEXP
      P0T=P0(M,XA)
      IESS=IEXP
      DIN=FIEXP(DIN*P0T)
      IEXP=IES+IEXP+IESS
      GO TO 35
   45 DIN=DTN(4,1,M,XA)
      QX=QEXP(XA)
      IES=IEXP
      DIN=FIEXP(DIN*QX)
      IEXP=IEXP+IES
   35 DINF=DIN
      IINF=IEXP
      XTI=XA
      MTI=M
      GO TO 36
    6 DIN=DINF
      IEXP=IINF
   36 MEO=M-2*(M/2)
      IF((IXS.LT.0).AND.(MEO.GT.0)) DIN=-DIN
   39 RETURN
      END
      FUNCTION DNN(N,X)
      DOUBLE PRECISION DNN,X,DNR,DNR1,DNR2,FCR,FIEXP,FEXPC,DTN
     1,DTEN(-2/3),XTN,DNNF,XA,DABS
      DATA DTEN/1.0D-24,1.0D-12,1.0D0,1.0D12,1.0D24,1.0D36/
      COMMON/CIEXP/IEXP,FEXPC
      XA=DABS(X)
      IXS=1
      IF(X-0.0D0)7,8,9
    8 DNN=-1.0D10
      IEXP=80000
      GO TO 35
    7 IXS=-1
    9 M=IABS(N)
      IF(XA-XTN)4,5,4
    5 IF(M-MTN)4,6,4
    4 IF(M-3)10,11,12
   10 IF(XA-15.0D0)41,41,42
   11 IF(XA-16.0D0)41,42,42
   12 IF(M-7)13,14,14
   13 IF(XA-17.0D0)44,42,42
   14 IF(M-10)15,43,16
   43 IF(XA-19.0D0)45,42,42
   15 IF(XA-18.0D0)45,42,42
   16 IF(M-14)17,18,18
   17 IF(XA-19.0D0)46,42,42
   18 IF(M-16)19,20,21
   19 IF(XA-21.0D0)47,42,42
   20 IF(XA-22.0D0)48,42,42
   21 IF(M-18)22,23,24
   22 IF(XA-24.0D0)48,42,42
   23 IF(XA-25.0D0)48,42,42
   24 FN2=FLOAT(M*M)
      XS=XA
      FKIN=(FN2-0.25E0)/(XS-0.5E0+SQRT(XS*(XS-1.0E0)+FN2))
      IF(FKIN-6.3E0)42,42,48
   41 IF(XA-6.0D0)55,60,60
   44 IF(XA-7.0D0)55,55,60
   45 IF(XA-9.0D0)55,55,61
   46 IF(XA-11.0D0)55,55,62
   47 IF(XA-13.0D0)55,55,63
   48 FN=FLOAT(M)
      XS=XA
      X2S=XS*XS
      FLMAX=X2S/(2.0E0*(FN+SQRT(FN*FN+X2S)))
      IF(FLMAX-3.0E0)55,55,64
   42 DNN=DTN(2,1,M,XA)
      GO TO 35
   55 DNN=DTN(2,0,M,XA)
      GO TO 35
   60 IF(XA-9.0D0)70,61,61
   61 IF(XA-11.0D0)71,71,62
   62 IF(XA-13.0D0)72,72,63
   63 IF(XA-15.0D0)73,73,64
   70 DNR1=DTN(2,0,8,XA)
      IER1=IEXP
      DNR2=DTN(2,0,7,XA)
      IER=IEXP
      NSS=7
      GO TO 75
   71 DNR1=DTN(2,0,12,XA)
      IER1=IEXP
      DNR2=DTN(2,0,11,XA)
      IER=IEXP
      NSS=11
      GO TO 75
   72 DNR1=DTN(2,0,15,XA)
      IER1=IEXP
      DNR2=DTN(2,0,14,XA)
      IER=IEXP
      NSS=14
      GO TO 75
   73 DNR1=DTN(2,0,17,XA)
      IER1=IEXP
      DNR2=DTN(2,0,16,XA)
      IER=IEXP
      NSS=16
      GO TO 75
   64 DNR1=DTN(2,1,0,XA)
      DNR2=DTN(2,1,1,XA)
      IER=0
      FCR=1.0D0
      NSS=1
      NST=M-1
      NSU=1
      GO TO 80
   75 NST=M+1
      NSU=-1
      FCR=DTEN((IER1-IER)/12)
      GO TO 80
   80 DO 9977 NS=NSS,NST,NSU
      DNR=FIEXP(((2.0D0*NS*DNR2)/XA)-DNR1*FCR)
      IER=IEXP+IER
      FCR=FEXPC
      DNR1=DNR2
      DNR2=DNR
 9977 CONTINUE
      DNN=DNR
      IEXP=IER
      GO TO 35
   35 DNNF=DNN
      INNF=IEXP
      XTN=XA
      MTN=M
      GO TO 36
    6 DNN=DNNF
      IEXP=INNF
   36 MEO=M-2*(M/2)
      IF((N.LT.0).AND.(MEO.GT.0)) DNN=-DNN
      IF((IXS.LT.0).AND.(MEO.GT.0)) DNN=-DNN
   39 RETURN
      END
      FUNCTION DKN(N,X)
      DOUBLE PRECISION DKN,X,FEXPC,FIEXP,DTN,QEXP,DK0NET(12/48),
     1DK1NET(12/48),DBLE,QX,XI,DJ,DKR,DKR1,DKR2,FCR,DK0T,DK1T,DK0TT,
     2DK1TT,DELTA,DELT,DEL(14),XXX,XX(14),A0(15),A1(15),XTEST,QFLOAT,
     3QFADD,XTK,XA,DABS,DKNF
      DIMENSION IEL(14)
      COMMON/CIEXP/IEXP,FEXPC
      COMMON/DKNET/DK0NET,DK1NET
      XA=DABS(X)
      IXS=1
      IF(X-0.0D0)7,8,9
    8 DKN=1.0D10
      IEXP=80000
      GO TO 35
    7 IXS=-1
    9 M=IABS(N)
      IF(XA-XTK)4,5,4
    5 IF(M-MTK)4,6,4
    4 IF(M-7)10,11,12
   10 IF(XA-3.0D0)20,20,21
   20 DKN=DTN(3,0,M,XA)
      GO TO 35
   21 IF(XA-16.0D0)22,23,23
   22 IF(XA-XTEST)54,50,54
   54 IF(XA-8.0D0)41,42,42
   41 I=IDINT(4.0D0*XA)+1
      XI=FLOAT(I)/4.0D0
      IF(XA-XI+0.25D0)44,45,44
   42 I=IDINT(2.0D0*XA)+17
      XI=(FLOAT(I)/2.0D0)-8.0D0
      IF(XA-XI+0.5D0)44,45,44
   45 DK0T=DK0NET(I-1)
      DK1T=DK1NET(I-1)
      GO TO 50
   44 DK0T=DK0NET(I)
      DK1T=DK1NET(I)
      DELTA=(XA-XI)
      DELT=1.0D0
      IELT=0
      XXX=1.0D0
      DO 8888 J=1,14
      DJ=DBLE(FLOAT(J))
      DELT=FIEXP((DELT*DELTA)/DJ)
      IELT=IELT+IEXP
      DEL(J)=DELT
      IEL(J)=IELT
      XXX=XXX*XI
      XX(J)=1.0D0/XXX
 8888 CONTINUE
      A0(1)=0.0D0
      A1(1)=-1.0D0
      A0(2)=1.0D0
      A1(2)=XX(1)
      A0(3)=-XX(1)
      A1(3)=-(1.0D0+2.0D0*XX(2))
      A0(4)=1.0D0+3.0D0*XX(2)
      A1(4)=2.0D0*XX(1)+6.0D0*XX(3)
      A0(5)=-(2.0D0*XX(1)+12.0D0*XX(3))
      A1(5)=-(1.0D0+7.0D0*XX(2)+24.0D0*XX(4))
      A0(6)=1.0D0+9.0D0*XX(2)+60.0D0*XX(4)
      A1(6)=3.0D0*XX(1)+33.0D0*XX(3)+120.0D0*XX(5)
      A0(7)=-(3.0D0*XX(1)+51.0D0*XX(3)+360.0D0*XX(5))
      A1(7)=-(1.0D0+15.0D0*XX(2)+192.0D0*XX(4)+720.0D0*XX(6))
      A0(8)=1.0D0+18.0D0*XX(2)+345.0D0*XX(4)+2520.0D0*XX(6)
      A1(8)=4.0D0*XX(1)+96.0D0*XX(3)+1320.0D0*XX(5)+5040.0D0*XX(7)
      A0(9)=-(4.0D0*XX(1)+132.0D0*XX(3)+2700.0D0*XX(5)+20160.0D0*XX(7))
      A1(9)=-(1.0D0+26.0D0*XX(2)+729.0D0*XX(4)+10440.0D0*XX(6)+
     140320.0D0*XX(8))
      A0(10)=1.0D0+30.0D0*XX(2)+1125.0D0*XX(4)+23940.0D0*XX(6)+
     1181440.0D0*XX(8)
      A1(10)=5.0D0*XX(1)+210.0D0*XX(3)+6345.0D0*XX(5)+93240.0D0*XX(7)+
     1362880.0D0*XX(9)
      A0(11)=-(5.0D0*XX(1)+270.0D0*XX(3)+10845.0D0*XX(5)+236880.0D0*
     1XX(7)+1814400.0D0*XX(9))
      A1(11)=-(1.0D0+40.0D0*XX(2)+1965.0D0*XX(4)+62010.0D0*XX(6)+
     1927360.0D0*XX(8)+3628800.0D0*XX(10))
      A0(12)=1.0D0+45.0D0*XX(2)+2775.0D0*XX(4)+116235.0D0*XX(6)+
     12585520.0D0*XX(8)+19958400.0D0*XX(10)
      A1(12)=6.0D0*XX(1)+390.0D0*XX(3)+20670.0D0*XX(5)+
     1670950.0D0*XX(7)+10160640.0D0*XX(9)+39916800.0D0*XX(11)
      A0(13)=-(6.0D0*XX(1)+480.0D0*XX(3)+31770.0D0*XX(5)+
     11368360.0D0*XX(7)+30844800.0D0*XX(9)+239500800.0D0*XX(11))
      A1(13)=-(1.0D0+57.0D0*XX(2)+4335.0D0*XX(4)+240255.0D0*XX(6)+
     17953120.0D0*XX(8)+121564800.0D0*XX(10)+479001600.0D0*XX(12))
      A0(14)=1.0D0+63.0D0*XX(2)+5775.0D0*XX(4)+399105.0D0*XX(6)+
     117531640.0D0*XX(8)+399168000.0D0*XX(10)+3113510400.0D0*XX(12)
      A1(14)=7.0D0*XX(1)+651.0D0*XX(3)+53445.0D0*XX(5)+3050145.0D0*
     1XX(7)+102422880.0D0*XX(9)+1576713600.0D0*XX(11)+6227020800.0D0
     2*XX(13)
      A0(15)=-(7.0D0*XX(1)+777.0D0*XX(3)+76545.0D0*XX(5)+5444775.0D0*
     1XX(7)+242676000.0D0*XX(9)+5568393600.0D0*XX(11)+43589145600.0D0
     2*XX(13))
      A1(15)=-(1.0D0+77.0D0*XX(2)+8379.0D0*XX(4)+719775.0D0*XX(6)+
     141932800.0D0*XX(8)+1423396800.0D0*XX(10)+22034073600.0D0*XX(12)+
     287178291200.0D0*XX(14))
      DK0TT=0.0D0
      DK1TT=0.0D0
      I1=0
      I2=0
      DO 9996 J=14,1,-1
      DK0TT=QFADD(DK0TT,(A0(J)*DK0T+A1(J)*DK1T)*DEL(J),I1,IEL(J))
      I1=IEXP
      DK1TT=QFADD(DK1TT,-(A0(J+1)*DK0T+A1(J+1)*DK1T)*DEL(J),I2,IEL(J))
      I2=IEXP
 9996 CONTINUE
      DK0T=DK0T+QFLOAT(DK0TT,I1)
      DK1T=DK1T+QFLOAT(DK1TT,I2)
   50 XTEST=XA
      IF(M-1)51,52,53
   51 DKN=DK0T
      IEXP=0
      GO TO 35
   52 DKN=DK1T
      IEXP=0
      GO TO 35
   53 DKR1=DK0T
      DKR2=DK1T
      NSS=1
      IER=0
   55 FCR=1.0D0
      DO 9988 NS=NSS,M-1
      DKR=FIEXP(((2.0D0*NS*DKR2)/XA)+DKR1*FCR)
      IER=IEXP+IER
      FCR=FEXPC
      DKR1=DKR2
      DKR2=DKR
 9988 CONTINUE
      DKN=DKR
      IEXP=IER
      GO TO 35
   23 DKN=DTN(3,1,M,XA)
      QX=QEXP(-XA)
      IES=IEXP
      DKN=FIEXP(DKN*QX)
      IEXP=IES+IEXP
      GO TO 35
   11 IF(XA-5.0D0)20,20,21
   12 IF(M-9)13,14,15
   13 IF(XA-6.0D0)20,21,21
   14 IF(XA-7.0D0)20,21,21
   15 IF(M-21)16,16,17
   16 IF(XA-8.0D0)20,20,21
   17 IF(XA-8.0D0)20,20,25
   25 IF(XA-16.0D0)22,26,26
   26 IF(XA-20.0D0)24,24,23
   24 QX=QEXP(-XA)
      IER=IEXP
      DKR1=DTN(3,1,20,XA)*QX
      DKR2=DTN(3,1,21,XA)*QX
      NSS=21
      GO TO 55
   35 DKNF=DKN
      IKNF=IEXP
      XTK=XA
      MTK=M
      GO TO 36
    6 DKN=DKNF
      IEXP=IKNF
   36 MEO=M-2*(M/2)
      IF((IXS.LT.0).AND.(MEO.GT.0)) DKN=-DKN
   39 RETURN
      END
      FUNCTION DTN(IB,IIF,N,X)
      COMMON/CIEXP/IEXP,FEXPC
      DOUBLE PRECISION DTN,X,X2,XT2,XT22,A,DTA,FNM2,DTNO,DTNE,THETA,
     1DTOA,DTEA,PI2,A1,PI,DSIN,DSQRT,DBLE,X22,FK1,FK2,QATAN2,FNRED2,C,
     2PSIF,PSIA,PSIB,TA,DTNT,DFLN,DTNTA,DLOG,GAMMA,GAMMA2,P0,P0T,
     3FEXPC,FC1,FIEXP,FC1A,QFADD,XTEST,DCOS,PIT2,AMP
      FN=FLOAT(N)
      FN2=FN**2
      XS=SNGL(X)
      XS2=XS**2
      X2=X/2.0D0
      X22=X2**2
      X22S=X22
      XT2=X*2.0D0
      XT2S=XT2
      XT22=XT2**2
      XT22S=XT22
      PI=3.141592653589793D0
      PI2=1.570796326794897D0
      GAMMA=5.7721566490153285D-1
      GAMMA2=1.154431329803065721D0
      IF(IIF-1)40,45,45
   40 LMAX=INT((-FN+SQRT(FN2+XS2))/2.0E0)
      FLMAX=FLOAT(LMAX)
      FLN=FLMAX+FN
      LT=1
      T=1.0E0
   50 LT=LT+1
      T=T*(X22S/((FLMAX+LT)*(FLN+LT)))
      IF(T-1.0E-18)55,55,50
   55 LCUT=LMAX+LT
      LS=LCUT+N
      GO TO (201,205,205,204),IB
  201 DTN=0.0D0
  501 DO 5999 L=LCUT,LMAX+1,-1
      FL=FLOAT(L)
      DTN=-(X22/DBLE(FL*(FL+FN)))*(1.0D0+DTN)
 5999 CONTINUE
      IF(LMAX-1)556,557,558
  556 DTN=1.0D0+DTN
      GO TO 105
  557 DTN=((DTN+1.0D0)*(-X22/(FN+1.0D0)))+1.0D0
      GO TO 105
  558 LCUT=LMAX
      LMAX=1
      GO TO 501
  204 DTN=0.0D0
      DO 6999 L=LCUT,LMAX+1,-1
      FL=FLOAT(L)
      DTN=(X22/DBLE(FL*(FL+FN)))*(1.0D0+DTN)
 6999 CONTINUE
      IF(LMAX-1)656,657,658
  656 DTN=FIEXP(1.0D0+DTN)
      GO TO 105
  657 DTN=FIEXP(((DTN+1.0D0)*(X22/(FN+1.0D0)))+1.0D0)
      GO TO 105
  658 A=1.0D0
      IEX1=0
      DO 6888 L=1,LMAX
      FL=FLOAT(L)
      A=FIEXP(A*(X22/DBLE(FL*(FL+FN))))
      IEX1=IEX1+IEXP
 6888 CONTINUE
      DTA=1.0D0
      DO 6998 L=1,LMAX
      FL=FLOAT(L)
      DTA=DTA*(DBLE(FL*(FL+FN))/X22)+1.0D0
 6998 CONTINUE
      DTN=FIEXP(A*(DTN+DTA))
      IEXP=IEXP+IEX1
      GO TO 105
  205 PSIA=PSIF(1,LCUT)
      PSIB=PSIF(1,LS)
      IF(IB-3)202,203,203
  202 C=2.0D0*DLOG(X2)+GAMMA2
      TA=C-(PSIA+PSIB)
      DTN=0.0D0
  401 DO 9999 L=LCUT,LMAX+1,-1
      FL=FLOAT(L)
      DFLN=DBLE(FL+FN)
      DTN=-(X22*(TA+DTN))/(FL*DFLN)
      TA=TA+(1.0D0/DBLE(FL)+1.0D0/DFLN)
 9999 CONTINUE
      C=C-PSIF(1,N)
      IF(LMAX-1)56,57,58
   56 DTN=C+DTN
      GO TO 38
   57 DFLN=1.0D0+FN
      DTN=((DTN+C-1.0D0-1.0D0/DFLN)*(-X22/DFLN))+C
      GO TO 38
   58 IF(LMAX-3)458,459,459
  459 LCUT=LMAX
      LMAX=2
      GO TO 401
  458 A=((X22)**2)/(2.0D0*(FN+1.0D0)*(FN+2.0D0))
      DFLN=1.0D0+FN
      TA=C-1.0D0-1.0D0/DFLN
      DTA=-(C*DFLN)/X22+TA
      DFLN=2.0D0+FN
      TA=TA-0.5D0-1.0D0/DFLN
      DTA=-((DTA*(2.0D0*DFLN))/X22)+TA
      DTN=A*(DTN+DTA)
   38 DTN=DTN/PI
      P0T=P0(N,X)
      IEX1=IEXP
      DTN=FIEXP(DTN*P0T)
      IEX1A=IEXP+IEX1
   37 IF(N-1)44,41,42
   44 DTA=0.0D0
      IEXP=IEX1A
      GO TO 105
   41 DTA=1.0D0
      GO TO 43
   42 DTA=0.0D0
      DO 9997 L=N-1,1,-1
      FL=FLOAT(L)
      DTA=(X22/DBLE(FL*(FN-FL)))*(1.0D0+DTA)
 9997 CONTINUE
      DTA=DTA+1.0D0
   43 DTA=FIEXP(-DTA/(P0T*FN*PI))
      IEX2=IEXP-IEX1
   39 DTN=QFADD(DTA,DTN,IEX2,IEX1A)
      GO TO 105
  203 C=DLOG(X2)+GAMMA
      TA=-C+0.5D0*(PSIA+PSIB)
      DTN=0.0D0
  701 DO 7999 L=LCUT,LMAX+1,-1
      FL=FLOAT(L)
      DFLN=DBLE(FL+FN)
      DTN=(X22*(TA+DTN))/(FL*DFLN)
      TA=TA-(0.5D0/DBLE(FL)+0.5D0/DFLN)
 7999 CONTINUE
      C=-C+0.5D0*PSIF(1,N)
      IF(LMAX-1)756,757,758
  756 DTN=C+DTN
      GO TO 738
  757 DFLN=1.0D0+FN
      DTN=((DTN+C+(0.5D0+0.5D0/DFLN))*(X22/DFLN))+C
      GO TO 738
  758 LCUT=LMAX
      LMAX=1
      GO TO 701
  738 DTN=DTN*((-1.0D0)**N)
      P0T=P0(N,X)
      IEX1=IEXP
      DTN=FIEXP(DTN*P0T)
      IEX1A=IEXP+IEX1
      IF(N-1)744,741,742
  744 DTA=0.0D0
      IEXP=IEX1A
      GO TO 105
  741 DTA=1.0D0
      GO TO 743
  742 DTA=0.0D0
      DO 7997 L=N-1,1,-1
      FL=FLOAT(L)
      DTA=(-X22/DBLE(FL*(FN-FL)))*(1.0D0+DTA)
 7997 CONTINUE
      DTA=DTA+1.0D0
  743 DTA=FIEXP(DTA/(P0T*FN*2.0D0))
      IEX2=IEXP-IEX1
  739 DTN=QFADD(DTA,DTN,IEX2,IEX1A)
      GO TO 105
   45 GO TO (144,144,145,145),IB
  144 IF(N-NTEST)145,216,145
  216 IF(X-XTEST)145,218,145
  145 IF(N-1)46,47,47
   46 KIN=0
      GO TO 60
   47 KIN=INT((FN2-0.25E0)/(XS-0.5E0+SQRT(XS*(XS-1.0E0)+FN2)))
   60 IMAX=KIN
      KOUT=INT(XS+0.5E0+SQRT(XS*(XS+1.0E0)+FN2))
      IF(KOUT-KIN-1)70,70,71
   70 KCUT=KOUT
      GO TO 81
   71 KT=1+KIN
      T=1.0E0
   76 KT=KT+1
      FK=FLOAT(KT)
      T=T*(((FN-0.5E0+FK)*(ABS(-FN-0.5E0+FK)))/(XT2S*FK))
      IF(T-1.0E-18)80,75,75
   75 IF(KOUT-KT)80,80,76
   80 KCUT=KT
   81 GO TO (215,215,213,214),IB
  215 IF(KCUT-(2*(KCUT/2)))82,83,82
   82 KCUTO=KCUT
      KCUTE=KCUT-1
      GO TO 85
   83 KCUTO=KCUT-1
      KCUTE=KCUT
   85 DTNO=0.0D0
      DO 9988 K=KCUTO,KIN+2,-2
      FK=FLOAT(K)
      DTNO=-(((FN2-((FK-0.5D0)**2))*(FN2-((FK-1.5D0)**2)))/(XT22*FK*
     1(FK-1.0D0)))*(1.0D0+DTNO)
 9988 CONTINUE
      DTNE=0.0D0
      DO 9987 K=KCUTE,KIN+2,-2
      FK=FLOAT(K)
      DTNE=-(((FN2-((FK-0.5D0)**2))*(FN2-((FK-1.5D0)**2)))/(XT22*FK*
     1(FK-1.0D0)))*(1.0D0+DTNE)
 9987 CONTINUE
      IF(KIN-1)86,87,88
   86 DTNE=DTNE+1.0D0
      GO TO 92
   87 DTNE=-((((FN2-0.25D0)*(FN2-2.25D0))/(2.0D0*XT22))*(1.0D0+DTNE))
     1+1.0D0
   92 DTNO=-((FN2-0.25D0)/XT2)*(1.0D0+DTNO)
      GO TO 97
   88 DTEA=1.0D0
      DTOA=1.0D0
      DO 9986 K=2,KIN,2
      FK=FLOAT(K)
      DTEA=-DTEA*((XT22*(FK-1.0D0)*FK)/((FN2-((FK-1.5D0)**2))*(FN2-((F
     1K-0.5D0)**2))))+1.0D0
      DTOA=-DTOA*((XT22*FK*(FK+1.0D0))/((FN2-((FK-0.5D0)**2))*(FN2-
     1((FK+0.5D0)**2))))+1.0D0
 9986 CONTINUE
      A=1.0D0
      C=1.0D0
      DO 9985 K=1,KIN
      FK=FLOAT(K)
      C=-C
      A=C*A*((FN2-((FK-0.5D0)**2))/(XT2*FK))
 9985 CONTINUE
      FK=FLOAT(KIN)
      C=-C
      A1=C*A*((FN2-((FK+0.5D0)**2))/(XT2*(FK+1.0D0)))
      IF(KIN-(2*(KIN/2)))96,95,96
   96 DTEA=-DTEA*((XT22*FK*(FK+1.0D0))/((FN2-((FK+0.5D0)**2))*(FN2-
     1((FK-0.5D0)**2))))+1.0D0
      DTNE=A1*(DTNE+DTEA)
      DTNO=A*(DTNO+DTOA)
      GO TO 97
   95 DTNE=A*(DTNE+DTEA)
      DTNO=A1*(DTNO+DTOA)
   97 KAOUT=INT(SQRT(FN2+XS2))
      KT=0
      T=1.0E0
  102 KT=KT+1
      FK=FLOAT(KT)
      FK2S=FK*2.0E0
      T=T*ABS(((FN2-((FK-0.5E0)**2))*FK2S*(FK2S-1.0E0))/(XT22S*(FK**2)))
      IF(T-1.0E-18)101,101,100
  100 IF(KAOUT-KT)101,101,102
  101 KACUT=KT
      A=0.0D0
      DO 9984 K=KACUT,1,-1
      FK=FLOAT(K)
      FK2=FK*2.0D0
      A=(((FN2-((FK-0.5D0)**2))*FK2*(FK2-1.0D0))/(XT22*(FK**2)))*(
     11.0D0+A)
 9984 CONTINUE
      AMP=DSQRT((1.0D0+A)/(PI2*X))
      FNRED2=FLOAT(N-(4*(N/4)))/2.0D0
      THETA=X-((FNRED2+0.25D0)*PI)+QATAN2(-DTNO,DTNE)
      NTEST=N
      XTEST=X
  218 IF(IB-2)211,212,212
  211 DTN=AMP*DCOS(THETA)
      IEXP=0
      GO TO 105
  212 DTN=AMP*DSIN(THETA)
      IEXP=0
      GO TO 105
  214 DTN=0.0D0
      PIT2=6.28318530717958648D0
      DO 6988 K=KCUT,KIN+1,-1
      FK=FLOAT(K)
      DTN=-(((FN2-((FK-0.5D0)**2)))/(XT2*FK))*(1.0D0+DTN)
 6988 CONTINUE
      IF(KIN-1)686,687,688
  686 DTN=DTN+1.0D0
      GO TO 697
  687 DTN=-(((FN2-0.25D0)/XT2)*(1.0D0+DTN))+1.0D0
      GO TO 697
  688 DTA=1.0D0
      DO 6986 K=1,KIN,1
      FK=FLOAT(K)
      DTA=-DTA*((XT2*FK)/((FN2-((FK-0.5D0)**2))))+1.0D0
 6986 CONTINUE
      A=1.0D0
      DO 6985 K=1,KIN
      FK=FLOAT(K)
      A=-A*((FN2-((FK-0.5D0)**2))/(XT2*FK))
 6985 CONTINUE
      DTN=A*(DTN+DTA)
  697 DTN=DTN/DSQRT(PIT2*X)
      GO TO 105
  213 DTN=0.0D0
      DO 7988 K=KCUT,KIN+1,-1
      FK=FLOAT(K)
      DTN=(((FN2-((FK-0.5D0)**2)))/(XT2*FK))*(1.0D0+DTN)
 7988 CONTINUE
      IF(KIN-1)786,787,788
  786 DTN=DTN+1.0D0
      GO TO 797
  787 DTN=(((FN2-0.25D0)/XT2)*(1.0D0+DTN))+1.0D0
      GO TO 797
  788 DTA=1.0D0
      DO 7986 K=1,KIN,1
      FK=FLOAT(K)
      DTA=DTA*((XT2*FK)/((FN2-((FK-0.5D0)**2))))+1.0D0
 7986 CONTINUE
      A=1.0D0
      DO 7985 K=1,KIN
      FK=FLOAT(K)
      A=A*((FN2-((FK-0.5D0)**2))/(XT2*FK))
 7985 CONTINUE
      DTN=A*(DTN+DTA)
  797 DTN=DTN*DSQRT(PI2/X)
      IEXP=0
  105 RETURN
      END
      DOUBLE PRECISION FUNCTION PSIF(I,J)
      DOUBLE PRECISION DBLE
      PSIF=0.0D0
      IF(I)15,15,5
    5 IF(J-I)15,10,10
   10 DO 9999 N=J,I,-1
      PSIF=(1.0D0/DBLE(FLOAT(N)))+PSIF
 9999 CONTINUE
   15 RETURN
      END
      FUNCTION P0(N,X)
      DOUBLE PRECISION P0,X,X2,DBLE,DFI,FIEXP,FEXPC
      COMMON/CIEXP/IEXP,FEXPC
      IF(N)10,15,20
   10 P0=0.0D0
      IEXP=0
      GO TO 30
   15 P0=1.0D0
      IEXP=0
      GO TO 30
   20 X2=X/2.0D0
      P0=1.0D0
      IEXP0=0
      DO 9999 I=1,N
      DFI=DBLE(FLOAT(I))
      P0=FIEXP((P0*X2)/DFI)
      IEXP0=IEXP0+IEXP
 9999 CONTINUE
      IEXP=IEXP0
   30 RETURN
      END
      BLOCK DATA
      COMMON/DKNET/DK0NET,DK1NET
      DOUBLE PRECISION DK0NET(12/48),DK1NET(12/48)
      DATA DK0NET/ 0.3473950438627915D-01,
     1 0.2605875525515488D-01, 0.1959889717036844D-01,
     1 0.1477425087712867D-01, 0.1115967608585300D-01,
     1 0.8444387724542933D-02, 0.6399857243233959D-02,
     1 0.4857204557887941D-02, 0.3691098334042585D-02,
     1 0.2808185109509765D-02, 0.2138708565950282D-02,
     1 0.1630399241396202D-02, 0.1243994328013120D-02,
     1 0.9499365977602106D-03, 0.7259317676293338D-03,
     1 0.5551343320208508D-03, 0.4247957418692310D-03,
     1 0.3252541730869297D-03, 0.2491776163561139D-03,
     1 0.1909954548670204D-03, 0.1464707052228151D-03,
     1 0.8625756634932497D-04, 0.5088131295645919D-04,
     1 0.3005788495793431D-04, 0.1778006231616764D-04,
     1 0.1052998814386532D-04, 0.6243020547653675D-05,
     1 0.3705038165956420D-05, 0.2200825397311491D-05,
     1 0.1308403696776977D-05, 0.7784543861420498D-06,
     1 0.4634841671408218D-06, 0.2761370823981621D-06,
     1 0.1646200520299790D-06, 0.9819536482396438D-07,
     1 0.5860481626637400D-07, 0.3499411663936501D-07/
      DATA DK1NET/ 0.4015643112819406D-01,
     1 0.2982552979604004D-01, 0.2223939292592377D-01,
     1 0.1663819175468886D-01, 0.1248349888726840D-01,
     1 0.9389680656043224D-02, 0.7078094908968072D-02,
     1 0.5345921817822826D-02, 0.4044613445452155D-02,
     1 0.3064795171512017D-02, 0.2325569008849000D-02,
     1 0.1766863714470692D-02, 0.1343919717735506D-02,
     1 0.1023285588634633D-02, 0.7798943982238019D-03,
     1 0.5949175971903856D-03, 0.4541824868848961D-03,
     1 0.3470006230989680D-03, 0.2652973901252890D-03,
     1 0.2029633036737783D-03, 0.1553692118050009D-03,
     1 0.9119724775006886D-04, 0.5363701637945189D-04,
     1 0.3160203411042672D-04, 0.1864877345382557D-04,
     1 0.1102047231135389D-04, 0.6520860674580883D-05,
     1 0.3862894146160996D-05, 0.2290757464767187D-05,
     1 0.1359767843821518D-05, 0.8078588412202348D-06,
     1 0.4803535332788457D-06, 0.2858343653440250D-06,
     1 0.1702048453059960D-06, 0.1014172936976210D-06,
     1 0.6046659442305888D-07, 0.3607157117528781D-07/
      END