Initial commit from Polytechnique Montreal

1 files changed, 69 insertions, 0 deletions

diff --git a/Utilib/src/AK1BES.f b/Utilib/src/AK1BES.f
new file mode 100644
index 0000000..55ff28d
--- /dev/null
+++ b/Utilib/src/AK1BES.f
@@ -0,0 +1,69 @@
+*DECK AK1BES
+      FUNCTION AK1BES(X)
+C
+C FUNCTION TO COMPUTE BESSEL FUNCTION K1(X)
+C X   ARGUMENT, REAL AND POSITIVE
+C IF X INDEFINITE OR INFINITE, ROUTINE SETS RESULT TO INDEFINITE,
+C ISSUES MESSAGE, AND RETURNS
+C METHOD-FOR X.LT.1 , K1(X)=ALOG I1(X)+1/X+X*P3(X**2)/Q3(X**2)
+C          FOR X.GE.1, K1(X)=(1/SQRT(X))*EXP(-X)*P6(Z)/Q7(Z) WHERE Z=1/X
+C COEFFICIENTS OF P AND Q GIVEN IN AECL-3461
+C COEFFICIENTS OF I1(X) GIVEN IN AECL-3461
+C PROGRAMMED BY B. PURCELL     JULY 1975
+C MODIFICATION OF AIKBES
+C
+C TAKEN FROM *AELIB* LIBRARY
+C MODIFY BY R.ROY (JULY 1994) IN DOUBLE PRECISION
+C
+      DOUBLE PRECISION AK1BES,X,P(4),Q(3),R(5),H(8),V(7)
+      DOUBLE PRECISION A,B,C,D,E,F,G,U,Z,XSQ,SQRTX
+      SAVE
+      DATA P/ 2.05653783869908D5, 5.21490277398084D4,
+     >        1.80714545607523D3, 1.66644511883896D1 /
+      DATA Q/-6.67781331663377D5, 1.57807971896305D4,
+     >       -1.76644575202506D2 /
+      DATA R/-6.04290103419901D1, -7.05362629269085D0,
+     >       -2.52234429536922D-1,-3.9527977417420D-3,
+     >       -2.77141112612D-5   /
+      DATA H/ 5.05657074674356D1, 5.86274965274502D2,
+     >        2.43240313288218D3, 4.58937398896551D3,
+     >        4.14893728249324D3, 1.73179927703691D3,
+     >        2.87833733475985D2, 1.24115199887274D1 /
+      DATA V/ 4.03455972943410D1, 4.52650144192445D2,
+     >        1.77576110097515D3, 3.04478817648148D3,
+     >        2.33607774590470D3, 7.34544788193384D2,
+     >        7.34518391966755D1 /
+C
+C TEST FOR ARGUMENT.LE.0
+C
+      AK1BES=0.0D0
+      IF(X.LE.0.D0) GOTO 10
+      IF(X.GT.1.D0) GOTO 5
+C
+C THIS BRANCH FOR X.LT.1
+C
+      XSQ=X*X
+      A=((P(4)*XSQ+P(3))*XSQ+P(2))*XSQ+P(1)
+      B=((XSQ+Q(3))*XSQ+Q(2))*XSQ+Q(1)
+      C=A/B
+      D=(((R(5)*XSQ+R(4))*XSQ+R(3))*XSQ+R(2))*XSQ+R(1)
+      E=-1.20858020683979D2+XSQ
+      F=D/E
+      AK1BES= LOG(X)*X*F+1.D0/X+X*C
+      RETURN
+C
+C  THIS PART OF THE PROGRAM CALCULATES K1(X) FOR X.GT.1
+C
+    5 SQRTX=SQRT(X)
+      Z=1.D0/X
+      G=((((((H(8)*Z+H(7))*Z+H(6))*Z+H(5))*Z+H(4))*Z+H(3))*Z+H(2))*Z
+     >   + H(1)
+      U=((((((Z+V(7))*Z+V(6))*Z+V(5))*Z+V(4))*Z+V(3))*Z+V(2))*Z+V(1)
+      AK1BES=(1.D0/SQRTX)*(EXP(-X))*G/U
+      RETURN
+C
+C
+C ERROR SECTION FOR X.LE.0
+   10 CALL XABORT('AK1BES: ILLEGAL ARGUMENT FOR BESSEL K1')
+      RETURN
+      END