Initial commit from Polytechnique Montreal

1 files changed, 272 insertions, 0 deletions

diff --git a/Utilib/src/AKIN10.f b/Utilib/src/AKIN10.f
new file mode 100644
index 0000000..b7fe2eb
--- /dev/null
+++ b/Utilib/src/AKIN10.f
@@ -0,0 +1,272 @@
+*DECK AKIN10
+      SUBROUTINE AKIN10(X,AKIN)
+C
+C THE SUBROUTINE AKIN10 EVALUATES THE BICKLEY FUNCTIONS KIN(X),
+C N = 1,2,...,10, X \ 0, TO A MINIMUM PRECISION OF 12.75 SIGNIFICANT
+C FIGURES AND RETURNS THE VALUES KIN(X),N = 1,2,...,10 IN THE ARRAY AKIN
+C
+C
+C ARGUMENT LIST:
+C
+C     X - REAL, X \ 0.
+C     AKIN - REAL, ONE DIMENSIONAL ARRAY TO STORE KIN(X),N = 1,2,...,10
+C
+C
+C METHOD:
+C
+C THE BICKLEY FUNCTIONS SATISFY A RECURSIVE RELATION ENABLING A FUNCTION
+C OF ANY ORDER TO BE FOUND IF THREE BICKLEY FUNCTIONS OF CONSECUTIVE
+C ORDER ARE KNOWN.  THE FORMULA FOR FORWARD RECURSION IS:
+C     (KI(N))(X) = (X * ((KI(N-3))(X) - (KI(N-1))(X)) / (N-1))
+C                  + ((N-2) * (KI(N-2))(X)) / (N-1)
+C THE FORMULA FOR BACKWARD RECURSION IS:
+C     (KI(N))(X) = ((N+2) * (KI(N+3))(X)) / X + (KI(N+2))(X)
+C                  - ((N+1) * (KI(N+1))(X)) / X
+C FORWARD AND BACKWARD RECURSION WERE TESTED \ SEE "APPROXIMATING THE
+C BICKLEY FUNCTIONS KIN(X),N=1,2,...,10",C.A.EDWARDS!.  LEAST LOSS OF
+C PRECISION WAS OBTAINED BY USING FORWARD RECURSION FOR 0@X@6, AND
+C BACKWARD RECUSION FOR X\6.
+C
+C X \ 6.     KI,(X),KI9(X),KI10(X) ARE CALCULATED AND USED IN BACKWARD
+C RECURSION TO DERIVE KIN(X),N=7,6,...,1
+C
+C 0 < X < 6.  KI1(X),KI2(X),KI3(X) ARE CALCULATED AND USED IN FORWARD
+C RECURSION TO DERIVE KIN(X),N=4,5,...,10
+C
+C RATIONAL POLYNOMIAL APPROXIMATIONS TO THE BICKLEY FUNCTIONS WITH A
+C MINIMUM ACCURACY OF 15D WERE OBTAINED FROM "RATIONAL CHEBYSHEV
+C APPROXIMATIONS TO THE BICKLEY FUNCTIONS KIN(X)",AECL-     AND
+C "RATIONAL CHEBYSHEV  APPROXIMATIONS TO THE BICKLEY FUNCTION KI3(X),
+C X \ 6",AECL-5820.
+C
+C     X = 0.  KIN(X) ,N=1,2,...,10 IS APPROXIMATED BY CONSTANTS.
+C
+C     0 < X < 1.  KIN(X),N = 1,2,3 IS APPROXIMATED BY
+C          R(X) + X**N * ALOG(X) * S(X**2)
+C     WHERE R AND S ARE RATIONAL POLYNOMIALS IN X AND X**2 RESPECTIVELY.
+C
+C     1 @ X < 6.  KIN(X),N=1,2,3 IS APPROXIMATED BY
+C          EXP(-X) * R(1. / X) ) / SQRT(X)
+C     WHERE R IS A RATIONAL POLYNOMIAL IN (1. / X)
+C
+C     6 @ X < 673.5  KIN(X),N = 8,9,10 IS APPROXIMATED BY
+C          EXP(-X) * R(1. / X) ) / SQRT(X)
+C     WHERE R IS A RATIONAL POLYNOMIAL IN (1. / X)
+C
+C     X \ 673.5  KIN(X),1=1,10 IS APPROXIMATED BY ZERO
+C
+C
+C WRITTEN BY P.CHRISTIE, MARCH,1977.
+C
+C TAKEN FROM *AELIB* LIBRARY (SUBROUTINE NAME: KIN)
+C MODIFY BY R.ROY (JULY 1994) IN DOUBLE PRECISION
+C
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      DOUBLE PRECISION AKIN(10),LNX
+      SAVE
+C
+      DATA R1P01,R1P11/-.3916967515498982D+02 , .1501664584981108D+03 /
+      DATA R1P21,R1P31/ .1531257133183402D+02 , .1398682937631850D+02 /
+      DATA R1P41,R1P51/ .5454840912170553D+00 , .5720051721276178D+00 /
+      DATA R1P61,R1P71/ .8766262885587739D-02 , .1001220446004980D-01 /
+      DATA R1P81,R1P91/ .5167460909383400D-04 , .7945473334662959D-04 /
+      DATA R1Q01,R1Q11/-.1193155300614385D+03 , .1219596245756669D+01 /
+      DATA S1P01,S1P11/ .1752373600092810D+05 , .1234654843555450D+04 /
+      DATA S1P21,S1P31/ .3695696751512241D+02 , .4646979722852471D+00 /
+      DATA S1P41      / .2337898258663651D-02 /
+      DATA S1Q01,S1Q11/ .1752373600092810D+05 ,-.2256564898552151D+03 /
+      DATA R2P01,R2P11/ .1833164538368226D+02 , .2427580524508585D+03 /
+      DATA R2P21,R2P31/ .1059016416894119D+04 , .1868734192859498D+04 /
+      DATA R2P41,R2P51/ .1326511766009986D+04 , .3272001530672078D+03 /
+      DATA R2P61      / .1822929159877549D+02 /
+      DATA R2Q01,R2Q11/ .1462653804563246D+02 , .2028344160151355D+03 /
+      DATA R2Q21,R2Q31/ .9570056876282360D+03 , .1922624187690926D+04 /
+      DATA R2Q41,R2Q51/ .1687760486772990D+04 , .5952592332227032D+03 /
+      DATA R2Q61      / .6590511376539962D+02 /
+      DATA R1P02,R1P12/ .4811961706232723D+03 , .1403058999999999D-10 /
+      DATA R1P22,R1P32/ .1693807200769639D+04 , .7521727532834893D+02 /
+      DATA R1P42,R1P52/ .4905564077627560D+02 , .9084569646859357D+00 /
+      DATA R1P62,R1P72/ .1190096804348251D+01 ,-.9999999999999957D-15 /
+      DATA R1P82      / .1086764750096697D-01 /
+      DATA R1Q02,R1Q12/ .1810502008060146D+04 ,-.7855291318496802D+02 /
+      DATA R1Q22      /-.1893578319929816D+02 /
+      DATA S1P02      /-.1520774316867189D+09 /
+      DATA S1Q02,S1Q12/ .3041548633734379D+09 ,-.1267311930720872D+08 /
+      DATA S1Q22,S1Q32/ .2112186548895240D+06 ,-.3143123637802091D+03 /
+      DATA S1Q42      /-.5631701819761997D+02 /
+      DATA R2P02,R2P12/ .5766817227841408D+00 , .9534267279889207D+01 /
+      DATA R2P22,R2P32/ .5320941946830476D+02 , .1244041140268300D+03 /
+      DATA R2P42,R2P52/ .1239740741934670D+03 , .4801733781249936D+02 /
+      DATA R2P62,R2P72/ .5596498537189973D+01 , .8407469297501269D-01 /
+      DATA R2Q02,R2Q12/ .4601255143693006D+00 , .8124881079392082D+01 /
+      DATA R2Q22,R2Q32/ .5035544525458363D+02 , .1383006201574071D+03 /
+      DATA R2Q42,R2Q52/ .1754089481769652D+03 , .9747773947009136D+02 /
+      DATA R2Q62      / .2014290370371339D+02 /
+      DATA R1P03,R1P13/-.7043581454636306D+06 , .1269499275481224D+07 /
+      DATA R1P23,R1P33/-.3790509999999980D-08 , .1095667013274141D+07 /
+      DATA R1P43      / .1173330390873767D+05 /
+      DATA R1Q03,R1Q13/-.2961415636470914D+07 , .2053835980116203D+06 /
+      DATA R1Q23,R1Q33/ .6077621585261822D+05 ,-.3942586515380026D+04 /
+      DATA R1Q43,R1Q53/-.1335161914247710D+03 ,-.1755175839841900D+02 /
+      DATA R1Q63      /-.5812262590904993D+01 /
+      DATA S1P03,S1P13/ .3027965765686327D+04 , .3721363059831219D+02 /
+      DATA S1P23,S1P33/ .5562992588150486D+00 , .2631126488553487D-02 /
+      DATA S1Q03,S1Q13/ .1816779459411791D+05 ,-.2309130812632629D+03 /
+      DATA R2P03,R2P13/ .3057730071278506D+00 , .5484114159179143D+01 /
+      DATA R2P23,R2P33/ .3295675793097689D+02 , .8221915520880930D+02 /
+      DATA R2P43,R2P53/ .8663638447767516D+02 , .3523089103388543D+02 /
+      DATA R2P63,R2P73/ .4296855364707056D+01 , .6753131438170179D-01 /
+      DATA R2Q03,R2Q13/ .2439716682658748D+00 , .4772136519722949D+01 /
+      DATA R2Q23,R2Q33/ .3279841354751966D+02 , .9980505157913694D+02 /
+      DATA R2Q43,R2Q53/ .1397165642636965D+03 , .8502251629078864D+02 /
+      DATA R2Q63      / .1898628899818902D+02 /
+      DATA R1P08,R1P18/ .2613746296365794D-04 , .1418798411613728D-02 /
+      DATA R1P28,R1P38/ .2815600768166693D-01 , .2548459851246907D+00 /
+      DATA R1P48,R1P58/ .1072628589602992D+01 , .1911501312353472D+01 /
+      DATA R1P68,R1P78/ .1112575911492456D+01 , .8986422681366424D-01 /
+      DATA R1Q08,R1Q18/ .2085467815725935D-04 , .1218062894916834D-02 /
+      DATA R1Q28,R1Q38/ .2689392859262307D-01 , .2845497743872052D+00 /
+      DATA R1Q48,R1Q58/ .1507878566796656D+01 , .3799408842694803D+01 /
+      DATA R1Q68      / .3846344279573060D+01 /
+      DATA R1P09,R1P19/ .4238334538408917D-04 , .2673933491122202D-02 /
+      DATA R1P29,R1P39/ .6307055224810387D-01 , .7008898871037985D+00 /
+      DATA R1P49,R1P59/ .3808588292608732D+01 , .9546586415461775D+01 /
+      DATA R1P69,R1P79/ .9312166436472808D+01 , .2271732886275686D+01 /
+      DATA R1Q09,R1Q19/ .3381701691714022D-04 , .2289893952421851D-02 /
+      DATA R1Q29,R1Q39/ .5971301111405558D-01 , .7653242147414652D+00 /
+      DATA R1Q49,R1Q59/ .5098663790807535D+01 , .1715435708724829D+02 /
+      DATA R1Q69,R1Q79/ .2599209729530514D+02 , .1346934010366220D+02 /
+      DATA R1P010,R1P110/.2984857391503590D-04, .1966652996666221D-02 /
+      DATA R1P210,R1P310/.4842493199976428D-01, .5612690964007783D+00 /
+      DATA R1P410,R1P510/.3176720369338028D+01, .8279846050989022D+01 /
+      DATA R1P610,R1P710/.8380471295924849D+01, .2115457456609434D+01 /
+      DATA R1Q010,R1Q110/.2381571628879027D-04, .1691217608476688D-02 /
+      DATA R1Q210,R1Q310/.4627626853651856D-01, .6224227393849020D+00 /
+      DATA R1Q410,R1Q510/.4348767092235221D+01, .1531699714433142D+02 /
+      DATA R1Q610,R1Q710/.2421477559533582D+02, .1303324507927971D+02 /
+C
+C
+      IF (X.LT.0.D0) GO TO 400
+C
+      IF (X.EQ.0.D0) GO TO 300
+      IF (X.LT.1.D0) GO TO 200
+      IF (X.GE.673.5D0) GO TO 110
+      SQRTX = SQRT(X)
+      EXPX = EXP(-X)
+      XREC = 1.D0/ X
+      IF (X.LT.6.D0) GO TO 210
+C
+C
+C     6 @ X < 673.5  CALCULATE KI(X),KI9(X),KI10(X)
+C
+      AKIN(10) =(((((((((R1P710 * XREC + R1P610) * XREC + R1P510) * XREC
+     ,+ R1P410) * XREC + R1P310) * XREC + R1P210) * XREC + R1P110) *
+     ,XREC + R1P010)         /  ((((((((XREC + R1Q710) * XREC + R1Q610)
+     ,* XREC + R1Q510) * XREC + R1Q410) * XREC + R1Q310) * XREC + R1Q210
+     ,) * XREC + R1Q110) * XREC + R1Q010) ) / SQRTX) * EXPX
+C
+      AKIN(9) =(((((((((R1P79 * XREC + R1P69) * XREC + R1P59) * XREC
+     ,+ R1P49) * XREC + R1P39) * XREC + R1P29) * XREC + R1P19) * XREC
+     ,+ R1P09)         /  ((((((((XREC + R1Q79) * XREC + R1Q69) * XREC
+     ,+ R1Q59) * XREC + R1Q49) * XREC + R1Q39) * XREC + R1Q29) * XREC
+     ,+ R1Q19) * XREC + R1Q09) ) / SQRTX) * EXPX
+C
+      AKIN(8) =(((((((((R1P78 * XREC + R1P68) * XREC + R1P58) * XREC
+     ,+ R1P48) * XREC + R1P38) * XREC + R1P28) * XREC + R1P18) * XREC
+     ,+ R1P08)         /  (((((((XREC + R1Q68) * XREC + R1Q58) * XREC
+     ,+ R1Q48) * XREC + R1Q38) * XREC + R1Q28) * XREC + R1Q18) * XREC
+     ,+ R1Q08) ) / SQRTX) * EXPX
+C
+C     BACKWARD RECURSION
+C
+      DO 100 I = 1,7
+           J = 8 - I
+           AKIN(J) = ((DBLE(J + 2) * AKIN(J + 3)) - (DBLE(J + 1)
+     ,               * AKIN(J + 1))) / X + AKIN(J + 2)
+  100 CONTINUE
+      RETURN
+C
+C
+C     X \ 673.5  AUTOMATICALLY SET AKIN(I),I=1,10 EQUAL TO ZERO
+C
+  110 DO 120 I = 1,10
+           AKIN(I)=0.D0
+  120 CONTINUE
+      RETURN
+C
+C
+C     0 < X < 1.  CALCULATE KI1(X),KI2(X),KI3(X)
+C
+  200 XSQ = X*X
+      LNX = LOG(X)
+      XI1 = X - .91169076342161D0
+      XI2 = X - .2451000192866D0
+      XI3 = X - 1.D0
+      XI4 = X - .68448452306295D0
+C
+      AKIN(1) = (((((((((R1P91 * X + R1P81) * X + R1P71) * X + R1P61)
+     ,* X + R1P51) * X + R1P41) * X + R1P31) * X + R1P21) * X + R1P11)
+     ,*XI3 + R1P01) / ((X + R1Q11) * XI3 + R1Q01) + (X * LNX * ((((S1P41
+     ,* XSQ + S1P31) * XSQ + S1P21) * XSQ + S1P11) * XSQ + S1P01) /
+     ,((XSQ + S1Q11) * XSQ + S1Q01))
+C
+      AKIN(2) = ((((R1P12 * XI1 + R1P02) + ((((R1P52 * X + R1P42) * X
+     ,+ R1P32) * X + R1P22) + (((R1P82 * XI2 + R1P72) * XI2 + R1P62)
+     ,* (X**4))) * (XI1 * XI1)) / (R1Q02 + (((        XI3 + R1Q22) * X
+     ,+ R1Q12) * XI3))) + (XSQ * LNX * ((S1P02) / (((((        XSQ
+     ,+ S1Q42) * XSQ + S1Q32) * XSQ + S1Q22) * XSQ + S1Q12) * XSQ
+     ,+ S1Q02))))
+C
+      AKIN(3) = (((((R1P43 * X + R1P33) * XI4 + R1P23) * XI4 + R1P13)
+     ,* XI3 + R1P03) /  (((((((XI3 + R1Q63) * X + R1Q53) * X + R1Q43)
+     ,* X + R1Q33) * XI3 + R1Q23) * X + R1Q13) * XI3 + R1Q03)) + (X**3
+     ,* LNX * ((((S1P33 * XSQ + S1P23) * XSQ + S1P13) * XSQ + S1P03) /
+     ,((XSQ + S1Q13) * XSQ + S1Q03)))
+C
+      GO TO 220
+C
+C
+C     1 @ X < 6.  CALCULATE KI1(X),KI2(X),KI3(X)
+C
+  210 AKIN(1) = (((((((R2P61 * XREC + R2P51) * XREC + R2P41) * XREC
+     ,+ R2P31) * XREC + R2P21) * XREC + R2P11) * XREC + R2P01) * EXPX) /
+     ,((((((((XREC + R2Q61) *XREC + R2Q51) * XREC + R2Q41) * XREC
+     ,+ R2Q31) * XREC + R2Q21) * XREC + R2Q11) * XREC + R2Q01) * SQRTX)
+C
+      AKIN(2) = ((((((((R2P72 * XREC + R2P62) * XREC + R2P52) * XREC
+     ,+ R2P42) * XREC + R2P32) * XREC + R2P22) * XREC + R2P12) * XREC
+     ,+ R2P02) * EXPX) / ((((((((XREC + R2Q62) * XREC + R2Q52) * XREC
+     ,+ R2Q42) * XREC + R2Q32) * XREC + R2Q22) * XREC + R2Q12) * XREC
+     ,+ R2Q02) * SQRTX)
+C
+      AKIN(3) = ((((((((R2P73 * XREC + R2P63) * XREC + R2P53) * XREC
+     ,+ R2P43) * XREC + R2P33) * XREC + R2P23) * XREC + R2P13) * XREC
+     ,+ R2P03) * EXPX) / ((((((((XREC + R2Q63) * XREC + R2Q53) * XREC
+     ,+ R2Q43) * XREC + R2Q33) * XREC + R2Q23) * XREC + R2Q13) * XREC
+     ,+ R2Q03) * SQRTX)
+C
+C     FORWARD RECURSION
+C
+  220 DO 230 I = 1,7
+           J = I + 3
+           AKIN(J) = (AKIN(J - 3) - AKIN(J - 1)) * X / DBLE(J - 1)
+     ,               + DBLE(J - 2) * AKIN(J - 2) / DBLE(J - 1)
+  230 CONTINUE
+      RETURN
+C
+C
+C     X = 0.
+C
+  300 AKIN(1) = 1.57079632679489D0
+      AKIN(2) = 1.D0
+      AKIN(3) = .78539816339745D0
+      AKIN(4) = .66666666666667D0
+      AKIN(5) = .58904862254808D0
+      AKIN(6) = .53333333333333D0
+      AKIN(7) = .49087385212340D0
+      AKIN(8) = .45714285714286D0
+      AKIN(9) = .42951462060798D0
+      AKIN(10)= .40634920634920D0
+      RETURN
+  400 CALL XABORT('AKIN10: ILLEGAL ARGUMENT FOR BICKLEY FUNCTIONS')
+      END