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