diff options
Diffstat (limited to 'Utilib/src/AKIN10.f')
| -rw-r--r-- | Utilib/src/AKIN10.f | 272 |
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 |