Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
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