Initial commit from Polytechnique Montreal

1 files changed, 68 insertions, 0 deletions

diff --git a/Utilib/src/AK0BES.f b/Utilib/src/AK0BES.f
new file mode 100644
index 0000000..ab78a21
--- /dev/null
+++ b/Utilib/src/AK0BES.f
@@ -0,0 +1,68 @@
+*DECK AK0BES
+      FUNCTION AK0BES(X)
+C
+C FUNCTION TO COMPUTE BESSEL FUNCTION K0(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 , K0(X)=-ALOG I0(X)+P3(X**2)/Q3(X**2)
+C          FOR X.GE.1, K0(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 I0(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 AK0BES,X,P(4),Q(3),R(4),S(3),H(7),V(7)
+      DOUBLE PRECISION A,B,C,D,E,F,G,U,Z,XSQ,SQRTX
+      SAVE
+      DATA P/-3.95157677473120D4, -9.39640131154667D4,
+     >       -5.90621631658113D3, -8.03837079456277D1 /
+      DATA Q/-3.40854404627503D5,  9.73445375796915D3,
+     >       -1.35884289391378D2 /
+      DATA R/-4.61386437648537D5, -1.03155855077145D5,
+     >       -4.31585225221548D3, -4.73684318951378D1/
+      DATA S/-4.61386437648537D5,  1.21907543350014D4,
+     >       -1.54377747904204D2 /
+      DATA H/ 1.75426986999764D2,  1.81039368258996D3,
+     >        6.39989885473113D3,  9.57620781628533D3,
+     >        6.07187112523892D3,  1.42043536504859D3,
+     >        8.10958280123701D1 /
+      DATA V/ 1.39970484475278D2,  1.46198147887226D3,
+     >        5.27928649717172D3,  8.20807534468755D3,
+     >        5.59084277299757D3,  1.50955451467318D3,
+     >        1.21349566357318D2 /
+C
+C TEST FOR ARGUMENT.LE.0
+C
+      AK0BES=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(4)*XSQ+R(3))*XSQ+R(2))*XSQ+R(1)
+      E=((XSQ+S(3))*XSQ+S(2))*XSQ+S(1)
+      F=D/E
+      AK0BES=-LOG(X)*F+C
+      RETURN
+C
+C  THIS PART OF THE PROGRAM CALCULATES K0(X) FOR X.GT.1
+C
+    5 SQRTX=SQRT(X)
+      Z=1.D0/X
+      G=(((((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)
+      AK0BES=(1.D0/SQRTX)*(EXP(-X))*G/U
+      RETURN
+C
+C ERROR SECTION FOR X.LE.0
+   10 CALL XABORT('AK0BES: ILLEGAL ARGUMENT FOR BESSEL K0')
+      RETURN
+      END