@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