From d55dcfb5ae31249d42c99c6af1da476b7e14d14e Mon Sep 17 00:00:00 2001 From: wdaniau Date: Wed, 2 Dec 2009 16:54:22 +0000 Subject: [PATCH] 1) Added single precision versions of besri and besrj 2) Modification of corresponding interface in fvn_fnlib.f90 git-svn-id: https://lxsd.femto-st.fr/svn/fvn@61 b657c933-2333-4658-acf2-d3c7c2708721 --- fvn_fnlib/Makefile | 1 + fvn_fnlib/b1slr.f | 228 ++++++++++++++++++++++++++++++++++++++++++++++++ fvn_fnlib/besri.f | 97 ++++++++++++++++++++ fvn_fnlib/besrj.f | 96 ++++++++++++++++++++ fvn_fnlib/fvn_fnlib.f90 | 22 +++-- 5 files changed, 438 insertions(+), 6 deletions(-) create mode 100644 fvn_fnlib/b1slr.f create mode 100644 fvn_fnlib/besri.f create mode 100644 fvn_fnlib/besrj.f diff --git a/fvn_fnlib/Makefile b/fvn_fnlib/Makefile index 805cbd8..08a9dfe 100644 --- a/fvn_fnlib/Makefile +++ b/fvn_fnlib/Makefile @@ -32,6 +32,7 @@ dbetai.o dbeta.o dbide.o dbid.o \ dbie.o dbinom.o dbi.o dbsi0e.o \ dbsi1e.o dbsk0e.o dbsk1e.o dbskes.o \ db1slr.o dbesri.o dbesrj.o \ +b1slr.o besri.o besrj.o \ dcbrt.o dchi.o dchu.o dcinh.o \ dcin.o dci.o dcosdg.o dcot.o \ dcsevl.o ddaws.o de1.o dei.o \ diff --git a/fvn_fnlib/b1slr.f b/fvn_fnlib/b1slr.f new file mode 100644 index 0000000..1ca7cc7 --- /dev/null +++ b/fvn_fnlib/b1slr.f @@ -0,0 +1,228 @@ + SUBROUTINE B1SLR(X,NB,IZE,B,NCALC) +C THIS ROUTINE CALCULATES BESSEL FUNCTIONS I AND J OF REAL +C ARGUMENT AND INTEGER ORDER. +C +C +C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE +C +C X REAL ARGUMENT FOR WHICH I*S OR J*S +C ARE TO BE CALCULATED. IF I*S ARE TO BE CALCULATED, +C ABS(X) MUST BE LESS THAN EXPARG (WHICH SEE BELOW). +C NB INTEGER TYPE. 1 + HIGHEST ORDER TO BE CALCULATED. +C IT MUST BE POSITIVE. +C IZE INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1 +C IF I*S ARE TO BE CALCULATED. +C B REAL VECTOR OF LENGTH NB, NEED NOT BE +C INITIALIZED BY USER. IF THE ROUTINE TERMINATES +C NORMALLY (NCALC=NB), IT RETURNS J(OR I)-SUB-ZERO +C THROUGH J(OR I)-SUB-NB-MINUS-ONE OF X IN THIS +C VECTOR. +C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER. +C BEFORE USING THE RESULTS, THE USER SHOULD CHECK THAT +C NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO +C THE DESIRED ACCURACY. SEE ERROR RETURNS BELOW. +C +C +C EXPLANATION OF MACHINE-DEPENDENT CONSTANTS +C +C NSIG DECIMAL SIGNIFICANCE DESIRED. SHOULD BE SET TO +C IFIX(ALOG10(2)*NBIT+1), WHERE NBIT IS THE NUMBER OF +C BITS IN THE MANTISSA OF A REAL VARIABLE. +C SETTING NSIG LOWER WILL RESULT IN DECREASED ACCURACY +C WHILE SETTING NSIG HIGHER WILL INCREASE CPU TIME +C WITHOUT INCREASING ACCURACY. THE TRUNCATION ERROR +C IS LIMITED TO T=.5*10**-NSIG FOR J*S OF ORDER LESS +C THAN ARGUMENT, AND TO A RELATIVE ERROR OF T FOR +C I*S AND THE OTHER J*S. +C NTEN LARGEST INTEGER K SUCH THAT 10**K IS MACHINE- +C REPRESENTABLE IN SINGLE PRECISION. +C LARGEX UPPER LIMIT ON THE MAGNITUDE OF X. BEAR IN MIND +C THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS OF THE +C BACKWARD RECURSION WILL BE EXECUTED. +C EXPARG LARGEST REAL ARGUMENT THAT THE LIBRARY +C EXP ROUTINE CAN HANDLE. +C +C PORT NOTE, SEPTEMBER 8,1976 - +C THE LARGEX AND EXPARG TESTS ARE MADE IN THE OUTER ROUTINES - +C BESRJ AND BESRI, WHICH CALL B1SLR. +C +C +C ERROR RETURNS +C +C PORT NOTE, SEPTEMBER 8, 1976 - +C THE NOTES BELOW ARE KEPT IN FOR THE RECORD, BUT, AS ABOVE, +C THE ACTUAL TESTS ARE NOW IN THE OUTER CALLING ROUTINES. +C +C LET G DENOTE EITHER I OR J. +C IN CASE OF AN ERROR, NCALC.NE.NB, AND NOT ALL G*S +C ARE CALCULATED TO THE DESIRED ACCURACY. +C IF NCALC.LT.0, AN ARGUMENT IS OUT OF RANGE. NB.LE.0 +C OR IZE IS NEITHER 0 NOR 1 OR IZE=1 AND ABS(X).GE.EXPARG. +C IN THIS CASE, THE B-VECTOR IS NOT CALCULATED, AND NCALC +C IS SET TO MIN0(NB,0)-1 SO NCALC.NE.NB. +C NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(G- +C SUB-NB-OF-X/G-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB +C IS MUCH GREATER THAN MAGX. IN THIS CASE, B(N) IS CALCU- +C LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR +C NCALC.LT.N.LE.NB, PRECISION IS LOST. IF N.GT.NCALC AND +C ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT +C FIGURES OF B(N) ARE ERRONEOUS. IF THE USER WISHES TO +C CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN +C ASYMPTOTIC FORMULA FOR LARGE ORDER. +C + REAL + 1 X,B,P,TEST,TEMPA,TEMPB,TEMPC,SIGN,SUM,TOVER, + 2 PLAST,POLD,PSAVE,PSAVEL,R1MACH + DIMENSION B(NB) + DATA NSIG/0/, NTEN/0/ + IF(NSIG .NE. 0) GO TO 1 + NSIG = IFIX(-ALOG10(R1MACH(3))+1.) + NTEN = ALOG10(R1MACH(2)) + 1 TEMPA=ABS(X) + MAGX=IFIX((TEMPA)) +C + SIGN=FLOAT(1-2*IZE) + NCALC=NB +C USE 2-TERM ASCENDING SERIES FOR SMALL X + IF(TEMPA**4.LT..1E0**NSIG) GO TO 30 +C INITIALIZE THE CALCULATION OF P*S + NBMX=NB-MAGX + N=MAGX+1 + PLAST=1.E0 + P=FLOAT(2*N)/TEMPA +C CALCULATE GENERAL SIGNIFICANCE TEST + TEST=2.E0*1.E1**NSIG + IF(IZE.EQ.1.AND.2*MAGX.GT.5*NSIG) TEST=SQRT(TEST*P) + IF(IZE.EQ.1.AND.2*MAGX.LE.5*NSIG) TEST=TEST/1.585**MAGX + M=0 + IF(NBMX.LT.3) GO TO 4 +C CALCULATE P*S UNTIL N=NB-1. CHECK FOR POSSIBLE OVERFLOW. + TOVER=1.E1**(NTEN-NSIG) + NSTART=MAGX+2 + NEND=NB-1 + DO 3 N=NSTART,NEND + POLD=PLAST + PLAST=P + P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD + IF(P-TOVER) 3,3,5 + 3 CONTINUE +C CALCULATE SPECIAL SIGNIFICANCE TEST FOR NBMX.GT.2. + TEST=AMAX1(TEST,SQRT(PLAST*1.E1**NSIG)*SQRT(2.E0*P)) +C CALCULATE P*S UNTIL SIGNIFICANCE TEST PASSES + 4 N=N+1 + POLD=PLAST + PLAST=P + P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD + IF(P.LT.TEST) GO TO 4 + IF(IZE.EQ.1.OR.M.EQ.1) GO TO 12 +C FOR J*S, A STRONG VARIANT OF THE TEST IS NECESSARY. +C CALCULATE IT, AND CALCULATE P*S UNTIL THIS TEST IS PASSED. + M=1 + TEMPB=P/PLAST + TEMPC=FLOAT(N+1)/TEMPA + IF(TEMPB+1.E0/TEMPB.GT.2.E0*TEMPC)TEMPB=TEMPC+SQRT(TEMPC**2-1.E0) + TEST=TEST/SQRT(TEMPB-1.E0/TEMPB) + IF(P-TEST) 4,12,12 +C TO AVOID OVERFLOW, DIVIDE P*S BY TOVER. CALCULATE P*S +C UNTIL ABS(P).GT.1. + 5 TOVER=1.E1**NTEN + P=P/TOVER + PLAST=PLAST/TOVER + PSAVE=P + PSAVEL=PLAST + NSTART=N+1 + 6 N=N+1 + POLD=PLAST + PLAST=P + P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD + IF(P.LE.1.E0) GO TO 6 + TEMPB=FLOAT(2*N)/TEMPA + IF(IZE.EQ.1) GO TO 8 + TEMPC=.5E0*TEMPB + TEMPB=PLAST/POLD + IF(TEMPB+1.E0/TEMPB.GT.2.E0*TEMPC)TEMPB=TEMPC+SQRT(TEMPC**2-1.E0) +C CALCULATE BACKWARD TEST, AND FIND NCALC, THE HIGHEST N +C SUCH THAT THE TEST IS PASSED. + 8 TEST=.5E0*POLD*PLAST*(1.E0-1.E0/TEMPB**2)/1.E1**NSIG + P=PLAST*TOVER + N=N-1 + NEND=MIN0(NB,N) + DO 9 NCALC=NSTART,NEND + POLD=PSAVEL + PSAVEL=PSAVE + PSAVE=FLOAT(2*N)*PSAVEL/TEMPA-SIGN*POLD + IF(PSAVE*PSAVEL-TEST) 9,9,10 + 9 CONTINUE + NCALC=NEND+1 + 10 NCALC=NCALC-1 +C THE SUM B(1)+2B(3)+2B(5)... IS USED TO NORMALIZE. M, THE +C COEFFICIENT OF B(N), IS INITIALIZED TO 2 OR 0. + 12 N=N+1 + M=2*N-4*(N/2) +C INITIALIZE THE BACKWARD RECURSION AND THE NORMALIZATION +C SUM + TEMPB=0.E0 + TEMPA=1.E0/P + SUM=FLOAT(M)*TEMPA + NEND=N-NB + IF(NEND) 17,15,13 +C RECUR BACKWARD VIA DIFFERENCE EQUATION, CALCULATING (BUT +C NOT STORING) B(N), UNTIL N=NB. + 13 DO 14 L=1,NEND + N=N-1 + TEMPC=TEMPB + TEMPB=TEMPA + TEMPA=FLOAT(2*N)*TEMPB/X-SIGN*TEMPC + M=2-M + 14 SUM=SUM+FLOAT(M)*TEMPA +C STORE B(NB) + 15 B(N)=TEMPA + IF(NB.GT.1) GO TO 16 +C NB=1. SINCE 2*TEMPA WAS ADDED TO THE SUM, TEMPA MUST BE +C SUBTRACTED + SUM=SUM-TEMPA + GO TO 23 +C CALCULATE AND STORE B(NB-1) + 16 N=N-1 + B(N) =FLOAT(2*N)*TEMPA/X-SIGN*TEMPB + IF(N.EQ.1) GO TO 22 + M=2-M + SUM=SUM+FLOAT(M)*B(N) + GO TO 19 +C N.LT.NB, SO STORE B(N) AND SET HIGHER ORDERS TO ZERO + 17 B(N)=TEMPA + NEND=-NEND + DO 18 L=1,NEND + K=N+L + 18 B(K)=0.E0 + 19 NEND=N-2 + IF(NEND.EQ.0) GO TO 21 +C CALCULATE VIA DIFFERENCE EQUATION AND STORE B(N), +C UNTIL N=2 + DO 20 L=1,NEND + N=N-1 + B(N)=(FLOAT(2*N)*B(N+1))/X-SIGN*B(N+2) + M=2-M + 20 SUM=SUM+FLOAT(M)*B(N) +C CALCULATE B(1) + 21 B(1)=2.E0*B(2)/X-SIGN*B(3) + 22 SUM=SUM+B(1) +C NORMALIZE--IF IZE=1, DIVIDE SUM BY COSH(X). DIVIDE ALL +C B(N) BY SUM. + 23 IF(IZE.EQ.0) GO TO 25 + TEMPA=EXP(ABS(X)) + SUM=2.E0*SUM/(TEMPA+1.E0/TEMPA) + 25 DO 26 N=1,NB + 26 B(N)=B(N)/SUM + RETURN +C +C TWO-TERM ASCENDING SERIES FOR SMALL X + 30 TEMPA=1.E0 + TEMPB=-.25E0*X*X*SIGN + B(1)=1.E0+TEMPB + IF(NB.EQ.1) GO TO 32 + DO 31 N=2,NB + TEMPA=TEMPA*X/FLOAT(2*N-2) + 31 B(N)=TEMPA*(1.E0+TEMPB/FLOAT(N)) + 32 RETURN + END diff --git a/fvn_fnlib/besri.f b/fvn_fnlib/besri.f new file mode 100644 index 0000000..225b5ab --- /dev/null +++ b/fvn_fnlib/besri.f @@ -0,0 +1,97 @@ + SUBROUTINE BESRI(X, NB, B) +C +C THIS ROUTINE CALCULATES MODIFIED BESSEL FUNCTIONS I OF REAL +C ARGUMENT AND INTEGER ORDER. +C +C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE - +C +C X - REAL ARGUMENT FOR WHICH THE I*S +C ARE TO BE CALCULATED. +C +C NB - INTEGER TYPE. 1 + HIGHEST ORDER TO BE CALCULATED. +C IT MUST BE POSITIVE. +C +C B - REAL VECTOR OF LENGTH NB, NEED NOT BE +C INITIALIZED BY USER. IF THE ROUTINE TERMINATES +C NORMALLY IT RETURNS I-SUB-ZERO +C THROUGH I-SUB-NB-MINUS-ONE OF X IN THIS +C VECTOR. +C +C ACCURACY OF THE COMPUTED VALUES - +C +C IN CASE OF AN ERROR, NOT ALL I*S +C ARE CALCULATED TO THE DESIRED ACCURACY. +C +C THE SUBPROGRAM CALLED BY BESRI, B1SLR, +C RETURNS IN THE VARIABLE, NCALC, THE NUMBER CALCULATED CORRECTLY. +C +C LET NTEN BE THE LARGEST INTEGER K SUCH THAT 10**K IS MACHINE- +C REPRESENTABLE IN REAL. +C THEN NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(I- +C SUB-NB-OF-X/I-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB +C IS MUCH GREATER THAN MAGX. IN THIS CASE, B(N) IS CALCU- +C LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR +C NCALC.LT.N.LE.NB, PRECISION IS LOST. IF N.GT.NCALC AND +C ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT +C FIGURES OF B(N) ARE ERRONEOUS. IF THE USER WISHES TO +C CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN +C ASYMPTOTIC FORMULA FOR LARGE ORDER. +C + REAL X,B(NB),R1MACH +C +C CHECK INPUT VALUES +C +C AN UPPER LIMIT OF 10000 IS SET ON THE MAGNITUDE OF X. +C BEAR IN MIND THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS +C OF THE BACKWARD RECURSION WILL BE EXECUTED. +C + MAGX = ABS(X) + MEXP = ALOG(R1MACH(2)) +C +C/6S +C IF (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR( +C 1 33H BESRI - X IS TOO BIG (MAGNITUDE),33,1,2) +C/7S + IF (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR( + 1 ' BESRI - X IS TOO BIG (MAGNITUDE)',33,1,2) +C/ +C +C/6S +C IF (NB .LT. 1) CALL SETERR( +C 1 28H BESRI - NB SHOULD = ORDER+1,28,2,2) +C/7S + IF (NB .LT. 1) CALL SETERR( + 1 ' BESRI - NB SHOULD = ORDER+1',28,2,2) +C/ +C +C BESRJ CALLS ON THE SUBPROGRAM,B1SLR, +C WHICH IS SOOKNES ORIGINAL BESLRI. +C +C THE ADDITIONAL INPUT ARGUMENTS REQUIRED FOR IT ARE - +C +C IZE INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1 +C IF I*S ARE TO BE CALCULATED.(THIRD ARGUMENT BELOW) +C +C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER. +C BEFORE USING THE RESULTS, IT SHOULD BE CHECKED THAT +C NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO +C THE DESIRED ACCURACY. +C + CALL B1SLR (X, NB, 1, B, NCALC) +C +C TEST IF ALL GOT COMPUTED OK +C (SINCE SOME VALUES MAY BE OK, THIS IS A RECOVERABLE ERROR.) +C + IF (NB .EQ. NCALC) RETURN +C + NCALC = NCALC+10 +C/6S +C CALL SETERR( +C 1 38H BESRI - ONLY THIS MANY ANSWERS ARE OK,38,NCALC,1) +C/7S + CALL SETERR( + 1 ' BESRI - ONLY THIS MANY ANSWERS ARE OK',38,NCALC,1) +C/ +C + RETURN + END diff --git a/fvn_fnlib/besrj.f b/fvn_fnlib/besrj.f new file mode 100644 index 0000000..1cd0333 --- /dev/null +++ b/fvn_fnlib/besrj.f @@ -0,0 +1,96 @@ + SUBROUTINE BESRJ(X, NB, B) +C +C THIS ROUTINE CALCULATES BESSEL FUNCTIONS J OF REAL +C ARGUMENT AND INTEGER ORDER. +C +C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE - +C +C X - REAL ARGUMENT FOR WHICH THE J*S +C ARE TO BE CALCULATED. +C +C NB - INTEGER TYPE. 1 + HIGHEST ORDER TO BE CALCULATED. +C IT MUST BE POSITIVE. +C +C B - REAL VECTOR OF LENGTH NB, NEED NOT BE +C INITIALIZED BY USER. IF THE ROUTINE TERMINATES +C NORMALLY IT RETURNS J(OR I)-SUB-ZERO +C THROUGH J(OR I)-SUB-NB-MINUS-ONE OF X IN THIS +C VECTOR. +C +C ACCURACY OF THE COMPUTED VALUES - +C +C IN CASE OF AN ERROR, NOT ALL J*S +C ARE CALCULATED TO THE DESIRED ACCURACY. +C +C THE SUBPROGRAM CALLED BY BESRJ, B1SLR, +C RETURNS IN THE VARIABLE, NCALC, THE NUMBER CALCULATED CORRECTLY. +C +C LET NTEN BE THE LARGEST INTEGER K SUCH THAT 10**K IS MACHINE- +C REPRESENTABLE IN SINGLE PRECISION. +C THEN NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(J- +C SUB-NB-OF-X/J-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB +C IS MUCH GREATER THAN MAGX. IN THIS CASE, B(N) IS CALCU- +C LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR +C NCALC.LT.N.LE.NB, PRECISION IS LOST. IF N.GT.NCALC AND +C ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT +C FIGURES OF B(N) ARE ERRONEOUS. IF THE USER WISHES TO +C CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN +C ASYMPTOTIC FORMULA FOR LARGE ORDER. +C + REAL X,B(NB) +C +C CHECK INPUT VALUES +C +C AN UPPER LIMIT OF 10000 IS SET ON THE MAGNITUDE OF X. +C BEAR IN MIND THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS +C OF THE BACKWARD RECURSION WILL BE EXECUTED. +C + MAGX = ABS(X) +C +C/6S +C IF (MAGX .GT. 10000) CALL SETERR( +C 1 33H BESRJ - X IS TOO BIG (MAGNITUDE),33,1,2) +C/7S + IF (MAGX .GT. 10000) CALL SETERR( + 1 ' BESRJ - X IS TOO BIG (MAGNITUDE)',33,1,2) +C/ +C +C/6S +C IF (NB .LT. 1) CALL SETERR( +C 1 28H BESRJ - NB SHOULD = ORDER+1,28,2,2) +C/7S + IF (NB .LT. 1) CALL SETERR( + 1 ' BESRJ - NB SHOULD = ORDER+1',28,2,2) +C/ +C +C BESRJ CALLS ON THE SUBPROGRAM,B1SLR, +C WHICH IS SOOKNES ORIGINAL BESLRI. +C +C THE ADDITIONAL INPUT ARGUMENTS REQUIRED FOR IT ARE - +C +C IZE INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1 +C IF I*S ARE TO BE CALCULATED.(THIRD ARGUMENT BELOW) +C +C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER. +C BEFORE USING THE RESULTS, IT SHOULD BE CHECKED THAT +C NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO +C THE DESIRED ACCURACY. +C + CALL B1SLR (X, NB, 0, B, NCALC) +C +C TEST IF ALL GOT COMPUTED OK +C (SINCE SOME VALUES MAY BE OK, THIS IS A RECOVERABLE ERROR.) +C + IF (NB .EQ. NCALC) RETURN +C + NCALC = NCALC+10 +C/6S +C CALL SETERR( +C 1 38H BESRJ - ONLY THIS MANY ANSWERS ARE OK,38,NCALC,1) +C/7S + CALL SETERR( + 1 ' BESRJ - ONLY THIS MANY ANSWERS ARE OK',38,NCALC,1) +C/ +C + RETURN + END diff --git a/fvn_fnlib/fvn_fnlib.f90 b/fvn_fnlib/fvn_fnlib.f90 index fb06c80..8a8236d 100644 --- a/fvn_fnlib/fvn_fnlib.f90 +++ b/fvn_fnlib/fvn_fnlib.f90 @@ -802,19 +802,29 @@ end interface bskn !!!!!!!!!!!!!!!!!!!!! ! vector b of Bessel J values of x from order 0 to order (n-1) interface besrj + subroutine besrj(x,n,b) + real(kind(1.e0)), intent(in) :: x + integer, intent(in) :: n + real(kind(1.e0)), intent(out) :: b(n) + end subroutine besrj subroutine dbesrj(x,n,b) - real(kind(1.d0)) :: x - integer :: n - real(kind(1.d0)) :: b(n) + real(kind(1.d0)), intent(in) :: x + integer, intent(in) :: n + real(kind(1.d0)), intent(out) :: b(n) end subroutine dbesrj end interface besrj ! vector b of Bessel I values of x from order 0 to order (n-1) interface besri + subroutine besri(x,n,b) + real(kind(1.e0)), intent(in) :: x + integer, intent(in) :: n + real(kind(1.e0)), intent(out) :: b(n) + end subroutine besri subroutine dbesri(x,n,b) - real(kind(1.d0)) :: x - integer :: n - real(kind(1.d0)) :: b(n) + real(kind(1.d0)), intent(in) :: x + integer, intent(in) :: n + real(kind(1.d0)), intent(out) :: b(n) end subroutine dbesri end interface besri -- 2.16.4