lxsd / fvn

Download as
Browse Code ยป

Commit d55dcfb5ae31249d42c99c6af1da476b7e14d14e

Authored by wdaniau
1 parent e80b2ec787
Exists in master and in 3 other branches geevx, geevx_2020, multi

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

Showing 5 changed files with 438 additions and 6 deletions Side-by-side Diff

... ... @@ -32,6 +32,7 @@
32 32 dbie.o dbinom.o dbi.o dbsi0e.o \
33 33 dbsi1e.o dbsk0e.o dbsk1e.o dbskes.o \
34 34 db1slr.o dbesri.o dbesrj.o \
  35 +b1slr.o besri.o besrj.o \
35 36 dcbrt.o dchi.o dchu.o dcinh.o \
36 37 dcin.o dci.o dcosdg.o dcot.o \
37 38 dcsevl.o ddaws.o de1.o dei.o \
  1 + SUBROUTINE B1SLR(X,NB,IZE,B,NCALC)
  2 +C THIS ROUTINE CALCULATES BESSEL FUNCTIONS I AND J OF REAL
  3 +C ARGUMENT AND INTEGER ORDER.
  4 +C
  5 +C
  6 +C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE
  7 +C
  8 +C X REAL ARGUMENT FOR WHICH I*S OR J*S
  9 +C ARE TO BE CALCULATED. IF I*S ARE TO BE CALCULATED,
  10 +C ABS(X) MUST BE LESS THAN EXPARG (WHICH SEE BELOW).
  11 +C NB INTEGER TYPE. 1 + HIGHEST ORDER TO BE CALCULATED.
  12 +C IT MUST BE POSITIVE.
  13 +C IZE INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1
  14 +C IF I*S ARE TO BE CALCULATED.
  15 +C B REAL VECTOR OF LENGTH NB, NEED NOT BE
  16 +C INITIALIZED BY USER. IF THE ROUTINE TERMINATES
  17 +C NORMALLY (NCALC=NB), IT RETURNS J(OR I)-SUB-ZERO
  18 +C THROUGH J(OR I)-SUB-NB-MINUS-ONE OF X IN THIS
  19 +C VECTOR.
  20 +C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER.
  21 +C BEFORE USING THE RESULTS, THE USER SHOULD CHECK THAT
  22 +C NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO
  23 +C THE DESIRED ACCURACY. SEE ERROR RETURNS BELOW.
  24 +C
  25 +C
  26 +C EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
  27 +C
  28 +C NSIG DECIMAL SIGNIFICANCE DESIRED. SHOULD BE SET TO
  29 +C IFIX(ALOG10(2)*NBIT+1), WHERE NBIT IS THE NUMBER OF
  30 +C BITS IN THE MANTISSA OF A REAL VARIABLE.
  31 +C SETTING NSIG LOWER WILL RESULT IN DECREASED ACCURACY
  32 +C WHILE SETTING NSIG HIGHER WILL INCREASE CPU TIME
  33 +C WITHOUT INCREASING ACCURACY. THE TRUNCATION ERROR
  34 +C IS LIMITED TO T=.5*10**-NSIG FOR J*S OF ORDER LESS
  35 +C THAN ARGUMENT, AND TO A RELATIVE ERROR OF T FOR
  36 +C I*S AND THE OTHER J*S.
  37 +C NTEN LARGEST INTEGER K SUCH THAT 10**K IS MACHINE-
  38 +C REPRESENTABLE IN SINGLE PRECISION.
  39 +C LARGEX UPPER LIMIT ON THE MAGNITUDE OF X. BEAR IN MIND
  40 +C THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS OF THE
  41 +C BACKWARD RECURSION WILL BE EXECUTED.
  42 +C EXPARG LARGEST REAL ARGUMENT THAT THE LIBRARY
  43 +C EXP ROUTINE CAN HANDLE.
  44 +C
  45 +C PORT NOTE, SEPTEMBER 8,1976 -
  46 +C THE LARGEX AND EXPARG TESTS ARE MADE IN THE OUTER ROUTINES -
  47 +C BESRJ AND BESRI, WHICH CALL B1SLR.
  48 +C
  49 +C
  50 +C ERROR RETURNS
  51 +C
  52 +C PORT NOTE, SEPTEMBER 8, 1976 -
  53 +C THE NOTES BELOW ARE KEPT IN FOR THE RECORD, BUT, AS ABOVE,
  54 +C THE ACTUAL TESTS ARE NOW IN THE OUTER CALLING ROUTINES.
  55 +C
  56 +C LET G DENOTE EITHER I OR J.
  57 +C IN CASE OF AN ERROR, NCALC.NE.NB, AND NOT ALL G*S
  58 +C ARE CALCULATED TO THE DESIRED ACCURACY.
  59 +C IF NCALC.LT.0, AN ARGUMENT IS OUT OF RANGE. NB.LE.0
  60 +C OR IZE IS NEITHER 0 NOR 1 OR IZE=1 AND ABS(X).GE.EXPARG.
  61 +C IN THIS CASE, THE B-VECTOR IS NOT CALCULATED, AND NCALC
  62 +C IS SET TO MIN0(NB,0)-1 SO NCALC.NE.NB.
  63 +C NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(G-
  64 +C SUB-NB-OF-X/G-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB
  65 +C IS MUCH GREATER THAN MAGX. IN THIS CASE, B(N) IS CALCU-
  66 +C LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR
  67 +C NCALC.LT.N.LE.NB, PRECISION IS LOST. IF N.GT.NCALC AND
  68 +C ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT
  69 +C FIGURES OF B(N) ARE ERRONEOUS. IF THE USER WISHES TO
  70 +C CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN
  71 +C ASYMPTOTIC FORMULA FOR LARGE ORDER.
  72 +C
  73 + REAL
  74 + 1 X,B,P,TEST,TEMPA,TEMPB,TEMPC,SIGN,SUM,TOVER,
  75 + 2 PLAST,POLD,PSAVE,PSAVEL,R1MACH
  76 + DIMENSION B(NB)
  77 + DATA NSIG/0/, NTEN/0/
  78 + IF(NSIG .NE. 0) GO TO 1
  79 + NSIG = IFIX(-ALOG10(R1MACH(3))+1.)
  80 + NTEN = ALOG10(R1MACH(2))
  81 + 1 TEMPA=ABS(X)
  82 + MAGX=IFIX((TEMPA))
  83 +C
  84 + SIGN=FLOAT(1-2*IZE)
  85 + NCALC=NB
  86 +C USE 2-TERM ASCENDING SERIES FOR SMALL X
  87 + IF(TEMPA**4.LT..1E0**NSIG) GO TO 30
  88 +C INITIALIZE THE CALCULATION OF P*S
  89 + NBMX=NB-MAGX
  90 + N=MAGX+1
  91 + PLAST=1.E0
  92 + P=FLOAT(2*N)/TEMPA
  93 +C CALCULATE GENERAL SIGNIFICANCE TEST
  94 + TEST=2.E0*1.E1**NSIG
  95 + IF(IZE.EQ.1.AND.2*MAGX.GT.5*NSIG) TEST=SQRT(TEST*P)
  96 + IF(IZE.EQ.1.AND.2*MAGX.LE.5*NSIG) TEST=TEST/1.585**MAGX
  97 + M=0
  98 + IF(NBMX.LT.3) GO TO 4
  99 +C CALCULATE P*S UNTIL N=NB-1. CHECK FOR POSSIBLE OVERFLOW.
  100 + TOVER=1.E1**(NTEN-NSIG)
  101 + NSTART=MAGX+2
  102 + NEND=NB-1
  103 + DO 3 N=NSTART,NEND
  104 + POLD=PLAST
  105 + PLAST=P
  106 + P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD
  107 + IF(P-TOVER) 3,3,5
  108 + 3 CONTINUE
  109 +C CALCULATE SPECIAL SIGNIFICANCE TEST FOR NBMX.GT.2.
  110 + TEST=AMAX1(TEST,SQRT(PLAST*1.E1**NSIG)*SQRT(2.E0*P))
  111 +C CALCULATE P*S UNTIL SIGNIFICANCE TEST PASSES
  112 + 4 N=N+1
  113 + POLD=PLAST
  114 + PLAST=P
  115 + P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD
  116 + IF(P.LT.TEST) GO TO 4
  117 + IF(IZE.EQ.1.OR.M.EQ.1) GO TO 12
  118 +C FOR J*S, A STRONG VARIANT OF THE TEST IS NECESSARY.
  119 +C CALCULATE IT, AND CALCULATE P*S UNTIL THIS TEST IS PASSED.
  120 + M=1
  121 + TEMPB=P/PLAST
  122 + TEMPC=FLOAT(N+1)/TEMPA
  123 + IF(TEMPB+1.E0/TEMPB.GT.2.E0*TEMPC)TEMPB=TEMPC+SQRT(TEMPC**2-1.E0)
  124 + TEST=TEST/SQRT(TEMPB-1.E0/TEMPB)
  125 + IF(P-TEST) 4,12,12
  126 +C TO AVOID OVERFLOW, DIVIDE P*S BY TOVER. CALCULATE P*S
  127 +C UNTIL ABS(P).GT.1.
  128 + 5 TOVER=1.E1**NTEN
  129 + P=P/TOVER
  130 + PLAST=PLAST/TOVER
  131 + PSAVE=P
  132 + PSAVEL=PLAST
  133 + NSTART=N+1
  134 + 6 N=N+1
  135 + POLD=PLAST
  136 + PLAST=P
  137 + P=FLOAT(2*N)*PLAST/TEMPA-SIGN*POLD
  138 + IF(P.LE.1.E0) GO TO 6
  139 + TEMPB=FLOAT(2*N)/TEMPA
  140 + IF(IZE.EQ.1) GO TO 8
  141 + TEMPC=.5E0*TEMPB
  142 + TEMPB=PLAST/POLD
  143 + IF(TEMPB+1.E0/TEMPB.GT.2.E0*TEMPC)TEMPB=TEMPC+SQRT(TEMPC**2-1.E0)
  144 +C CALCULATE BACKWARD TEST, AND FIND NCALC, THE HIGHEST N
  145 +C SUCH THAT THE TEST IS PASSED.
  146 + 8 TEST=.5E0*POLD*PLAST*(1.E0-1.E0/TEMPB**2)/1.E1**NSIG
  147 + P=PLAST*TOVER
  148 + N=N-1
  149 + NEND=MIN0(NB,N)
  150 + DO 9 NCALC=NSTART,NEND
  151 + POLD=PSAVEL
  152 + PSAVEL=PSAVE
  153 + PSAVE=FLOAT(2*N)*PSAVEL/TEMPA-SIGN*POLD
  154 + IF(PSAVE*PSAVEL-TEST) 9,9,10
  155 + 9 CONTINUE
  156 + NCALC=NEND+1
  157 + 10 NCALC=NCALC-1
  158 +C THE SUM B(1)+2B(3)+2B(5)... IS USED TO NORMALIZE. M, THE
  159 +C COEFFICIENT OF B(N), IS INITIALIZED TO 2 OR 0.
  160 + 12 N=N+1
  161 + M=2*N-4*(N/2)
  162 +C INITIALIZE THE BACKWARD RECURSION AND THE NORMALIZATION
  163 +C SUM
  164 + TEMPB=0.E0
  165 + TEMPA=1.E0/P
  166 + SUM=FLOAT(M)*TEMPA
  167 + NEND=N-NB
  168 + IF(NEND) 17,15,13
  169 +C RECUR BACKWARD VIA DIFFERENCE EQUATION, CALCULATING (BUT
  170 +C NOT STORING) B(N), UNTIL N=NB.
  171 + 13 DO 14 L=1,NEND
  172 + N=N-1
  173 + TEMPC=TEMPB
  174 + TEMPB=TEMPA
  175 + TEMPA=FLOAT(2*N)*TEMPB/X-SIGN*TEMPC
  176 + M=2-M
  177 + 14 SUM=SUM+FLOAT(M)*TEMPA
  178 +C STORE B(NB)
  179 + 15 B(N)=TEMPA
  180 + IF(NB.GT.1) GO TO 16
  181 +C NB=1. SINCE 2*TEMPA WAS ADDED TO THE SUM, TEMPA MUST BE
  182 +C SUBTRACTED
  183 + SUM=SUM-TEMPA
  184 + GO TO 23
  185 +C CALCULATE AND STORE B(NB-1)
  186 + 16 N=N-1
  187 + B(N) =FLOAT(2*N)*TEMPA/X-SIGN*TEMPB
  188 + IF(N.EQ.1) GO TO 22
  189 + M=2-M
  190 + SUM=SUM+FLOAT(M)*B(N)
  191 + GO TO 19
  192 +C N.LT.NB, SO STORE B(N) AND SET HIGHER ORDERS TO ZERO
  193 + 17 B(N)=TEMPA
  194 + NEND=-NEND
  195 + DO 18 L=1,NEND
  196 + K=N+L
  197 + 18 B(K)=0.E0
  198 + 19 NEND=N-2
  199 + IF(NEND.EQ.0) GO TO 21
  200 +C CALCULATE VIA DIFFERENCE EQUATION AND STORE B(N),
  201 +C UNTIL N=2
  202 + DO 20 L=1,NEND
  203 + N=N-1
  204 + B(N)=(FLOAT(2*N)*B(N+1))/X-SIGN*B(N+2)
  205 + M=2-M
  206 + 20 SUM=SUM+FLOAT(M)*B(N)
  207 +C CALCULATE B(1)
  208 + 21 B(1)=2.E0*B(2)/X-SIGN*B(3)
  209 + 22 SUM=SUM+B(1)
  210 +C NORMALIZE--IF IZE=1, DIVIDE SUM BY COSH(X). DIVIDE ALL
  211 +C B(N) BY SUM.
  212 + 23 IF(IZE.EQ.0) GO TO 25
  213 + TEMPA=EXP(ABS(X))
  214 + SUM=2.E0*SUM/(TEMPA+1.E0/TEMPA)
  215 + 25 DO 26 N=1,NB
  216 + 26 B(N)=B(N)/SUM
  217 + RETURN
  218 +C
  219 +C TWO-TERM ASCENDING SERIES FOR SMALL X
  220 + 30 TEMPA=1.E0
  221 + TEMPB=-.25E0*X*X*SIGN
  222 + B(1)=1.E0+TEMPB
  223 + IF(NB.EQ.1) GO TO 32
  224 + DO 31 N=2,NB
  225 + TEMPA=TEMPA*X/FLOAT(2*N-2)
  226 + 31 B(N)=TEMPA*(1.E0+TEMPB/FLOAT(N))
  227 + 32 RETURN
  228 + END
  1 + SUBROUTINE BESRI(X, NB, B)
  2 +C
  3 +C THIS ROUTINE CALCULATES MODIFIED BESSEL FUNCTIONS I OF REAL
  4 +C ARGUMENT AND INTEGER ORDER.
  5 +C
  6 +C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE -
  7 +C
  8 +C X - REAL ARGUMENT FOR WHICH THE I*S
  9 +C ARE TO BE CALCULATED.
  10 +C
  11 +C NB - INTEGER TYPE. 1 + HIGHEST ORDER TO BE CALCULATED.
  12 +C IT MUST BE POSITIVE.
  13 +C
  14 +C B - REAL VECTOR OF LENGTH NB, NEED NOT BE
  15 +C INITIALIZED BY USER. IF THE ROUTINE TERMINATES
  16 +C NORMALLY IT RETURNS I-SUB-ZERO
  17 +C THROUGH I-SUB-NB-MINUS-ONE OF X IN THIS
  18 +C VECTOR.
  19 +C
  20 +C ACCURACY OF THE COMPUTED VALUES -
  21 +C
  22 +C IN CASE OF AN ERROR, NOT ALL I*S
  23 +C ARE CALCULATED TO THE DESIRED ACCURACY.
  24 +C
  25 +C THE SUBPROGRAM CALLED BY BESRI, B1SLR,
  26 +C RETURNS IN THE VARIABLE, NCALC, THE NUMBER CALCULATED CORRECTLY.
  27 +C
  28 +C LET NTEN BE THE LARGEST INTEGER K SUCH THAT 10**K IS MACHINE-
  29 +C REPRESENTABLE IN REAL.
  30 +C THEN NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(I-
  31 +C SUB-NB-OF-X/I-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB
  32 +C IS MUCH GREATER THAN MAGX. IN THIS CASE, B(N) IS CALCU-
  33 +C LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR
  34 +C NCALC.LT.N.LE.NB, PRECISION IS LOST. IF N.GT.NCALC AND
  35 +C ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT
  36 +C FIGURES OF B(N) ARE ERRONEOUS. IF THE USER WISHES TO
  37 +C CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN
  38 +C ASYMPTOTIC FORMULA FOR LARGE ORDER.
  39 +C
  40 + REAL X,B(NB),R1MACH
  41 +C
  42 +C CHECK INPUT VALUES
  43 +C
  44 +C AN UPPER LIMIT OF 10000 IS SET ON THE MAGNITUDE OF X.
  45 +C BEAR IN MIND THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS
  46 +C OF THE BACKWARD RECURSION WILL BE EXECUTED.
  47 +C
  48 + MAGX = ABS(X)
  49 + MEXP = ALOG(R1MACH(2))
  50 +C
  51 +C/6S
  52 +C IF (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR(
  53 +C 1 33H BESRI - X IS TOO BIG (MAGNITUDE),33,1,2)
  54 +C/7S
  55 + IF (MAGX .GT. 10000 .OR. MAGX .GT. MEXP) CALL SETERR(
  56 + 1 ' BESRI - X IS TOO BIG (MAGNITUDE)',33,1,2)
  57 +C/
  58 +C
  59 +C/6S
  60 +C IF (NB .LT. 1) CALL SETERR(
  61 +C 1 28H BESRI - NB SHOULD = ORDER+1,28,2,2)
  62 +C/7S
  63 + IF (NB .LT. 1) CALL SETERR(
  64 + 1 ' BESRI - NB SHOULD = ORDER+1',28,2,2)
  65 +C/
  66 +C
  67 +C BESRJ CALLS ON THE SUBPROGRAM,B1SLR,
  68 +C WHICH IS SOOKNES ORIGINAL BESLRI.
  69 +C
  70 +C THE ADDITIONAL INPUT ARGUMENTS REQUIRED FOR IT ARE -
  71 +C
  72 +C IZE INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1
  73 +C IF I*S ARE TO BE CALCULATED.(THIRD ARGUMENT BELOW)
  74 +C
  75 +C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER.
  76 +C BEFORE USING THE RESULTS, IT SHOULD BE CHECKED THAT
  77 +C NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO
  78 +C THE DESIRED ACCURACY.
  79 +C
  80 + CALL B1SLR (X, NB, 1, B, NCALC)
  81 +C
  82 +C TEST IF ALL GOT COMPUTED OK
  83 +C (SINCE SOME VALUES MAY BE OK, THIS IS A RECOVERABLE ERROR.)
  84 +C
  85 + IF (NB .EQ. NCALC) RETURN
  86 +C
  87 + NCALC = NCALC+10
  88 +C/6S
  89 +C CALL SETERR(
  90 +C 1 38H BESRI - ONLY THIS MANY ANSWERS ARE OK,38,NCALC,1)
  91 +C/7S
  92 + CALL SETERR(
  93 + 1 ' BESRI - ONLY THIS MANY ANSWERS ARE OK',38,NCALC,1)
  94 +C/
  95 +C
  96 + RETURN
  97 + END
  1 + SUBROUTINE BESRJ(X, NB, B)
  2 +C
  3 +C THIS ROUTINE CALCULATES BESSEL FUNCTIONS J OF REAL
  4 +C ARGUMENT AND INTEGER ORDER.
  5 +C
  6 +C EXPLANATION OF VARIABLES IN THE CALLING SEQUENCE -
  7 +C
  8 +C X - REAL ARGUMENT FOR WHICH THE J*S
  9 +C ARE TO BE CALCULATED.
  10 +C
  11 +C NB - INTEGER TYPE. 1 + HIGHEST ORDER TO BE CALCULATED.
  12 +C IT MUST BE POSITIVE.
  13 +C
  14 +C B - REAL VECTOR OF LENGTH NB, NEED NOT BE
  15 +C INITIALIZED BY USER. IF THE ROUTINE TERMINATES
  16 +C NORMALLY IT RETURNS J(OR I)-SUB-ZERO
  17 +C THROUGH J(OR I)-SUB-NB-MINUS-ONE OF X IN THIS
  18 +C VECTOR.
  19 +C
  20 +C ACCURACY OF THE COMPUTED VALUES -
  21 +C
  22 +C IN CASE OF AN ERROR, NOT ALL J*S
  23 +C ARE CALCULATED TO THE DESIRED ACCURACY.
  24 +C
  25 +C THE SUBPROGRAM CALLED BY BESRJ, B1SLR,
  26 +C RETURNS IN THE VARIABLE, NCALC, THE NUMBER CALCULATED CORRECTLY.
  27 +C
  28 +C LET NTEN BE THE LARGEST INTEGER K SUCH THAT 10**K IS MACHINE-
  29 +C REPRESENTABLE IN SINGLE PRECISION.
  30 +C THEN NB.GT.NCALC.GT.0 WILL OCCUR IF NB.GT.MAGX AND ABS(J-
  31 +C SUB-NB-OF-X/J-SUB-MAGX+NP-OF-X).LT.10.**(NTEN/2), I.E. NB
  32 +C IS MUCH GREATER THAN MAGX. IN THIS CASE, B(N) IS CALCU-
  33 +C LATED TO THE DESIRED ACCURACY FOR N.LE.NCALC, BUT FOR
  34 +C NCALC.LT.N.LE.NB, PRECISION IS LOST. IF N.GT.NCALC AND
  35 +C ABS(B(NCALC)/B(N)).EQ.10**-K, THEN THE LAST K SIGNIFICANT
  36 +C FIGURES OF B(N) ARE ERRONEOUS. IF THE USER WISHES TO
  37 +C CALCULATE B(N) TO HIGHER ACCURACY, HE SHOULD USE AN
  38 +C ASYMPTOTIC FORMULA FOR LARGE ORDER.
  39 +C
  40 + REAL X,B(NB)
  41 +C
  42 +C CHECK INPUT VALUES
  43 +C
  44 +C AN UPPER LIMIT OF 10000 IS SET ON THE MAGNITUDE OF X.
  45 +C BEAR IN MIND THAT IF ABS(X)=N, THEN AT LEAST N ITERATIONS
  46 +C OF THE BACKWARD RECURSION WILL BE EXECUTED.
  47 +C
  48 + MAGX = ABS(X)
  49 +C
  50 +C/6S
  51 +C IF (MAGX .GT. 10000) CALL SETERR(
  52 +C 1 33H BESRJ - X IS TOO BIG (MAGNITUDE),33,1,2)
  53 +C/7S
  54 + IF (MAGX .GT. 10000) CALL SETERR(
  55 + 1 ' BESRJ - X IS TOO BIG (MAGNITUDE)',33,1,2)
  56 +C/
  57 +C
  58 +C/6S
  59 +C IF (NB .LT. 1) CALL SETERR(
  60 +C 1 28H BESRJ - NB SHOULD = ORDER+1,28,2,2)
  61 +C/7S
  62 + IF (NB .LT. 1) CALL SETERR(
  63 + 1 ' BESRJ - NB SHOULD = ORDER+1',28,2,2)
  64 +C/
  65 +C
  66 +C BESRJ CALLS ON THE SUBPROGRAM,B1SLR,
  67 +C WHICH IS SOOKNES ORIGINAL BESLRI.
  68 +C
  69 +C THE ADDITIONAL INPUT ARGUMENTS REQUIRED FOR IT ARE -
  70 +C
  71 +C IZE INTEGER TYPE. ZERO IF J*S ARE TO BE CALCULATED, 1
  72 +C IF I*S ARE TO BE CALCULATED.(THIRD ARGUMENT BELOW)
  73 +C
  74 +C NCALC INTEGER TYPE, NEED NOT BE INITIALIZED BY USER.
  75 +C BEFORE USING THE RESULTS, IT SHOULD BE CHECKED THAT
  76 +C NCALC=NB, I.E. ALL ORDERS HAVE BEEN CALCULATED TO
  77 +C THE DESIRED ACCURACY.
  78 +C
  79 + CALL B1SLR (X, NB, 0, B, NCALC)
  80 +C
  81 +C TEST IF ALL GOT COMPUTED OK
  82 +C (SINCE SOME VALUES MAY BE OK, THIS IS A RECOVERABLE ERROR.)
  83 +C
  84 + IF (NB .EQ. NCALC) RETURN
  85 +C
  86 + NCALC = NCALC+10
  87 +C/6S
  88 +C CALL SETERR(
  89 +C 1 38H BESRJ - ONLY THIS MANY ANSWERS ARE OK,38,NCALC,1)
  90 +C/7S
  91 + CALL SETERR(
  92 + 1 ' BESRJ - ONLY THIS MANY ANSWERS ARE OK',38,NCALC,1)
  93 +C/
  94 +C
  95 + RETURN
  96 + END
fvn_fnlib/fvn_fnlib.f90
Diff comments   View file @ d55dcfb
... ... @@ -802,19 +802,29 @@
802 802 !!!!!!!!!!!!!!!!!!!!!
803 803 ! vector b of Bessel J values of x from order 0 to order (n-1)
804 804 interface besrj
  805 + subroutine besrj(x,n,b)
  806 + real(kind(1.e0)), intent(in) :: x
  807 + integer, intent(in) :: n
  808 + real(kind(1.e0)), intent(out) :: b(n)
  809 + end subroutine besrj
805 810 subroutine dbesrj(x,n,b)
806   - real(kind(1.d0)) :: x
807   - integer :: n
808   - real(kind(1.d0)) :: b(n)
  811 + real(kind(1.d0)), intent(in) :: x
  812 + integer, intent(in) :: n
  813 + real(kind(1.d0)), intent(out) :: b(n)
809 814 end subroutine dbesrj
810 815 end interface besrj
811 816  
812 817 ! vector b of Bessel I values of x from order 0 to order (n-1)
813 818 interface besri
  819 + subroutine besri(x,n,b)
  820 + real(kind(1.e0)), intent(in) :: x
  821 + integer, intent(in) :: n
  822 + real(kind(1.e0)), intent(out) :: b(n)
  823 + end subroutine besri
814 824 subroutine dbesri(x,n,b)
815   - real(kind(1.d0)) :: x
816   - integer :: n
817   - real(kind(1.d0)) :: b(n)
  825 + real(kind(1.d0)), intent(in) :: x
  826 + integer, intent(in) :: n
  827 + real(kind(1.d0)), intent(out) :: b(n)
818 828 end subroutine dbesri
819 829 end interface besri
820 830