1      subroutine sjymec(x,j,y,mn,nout)
2c
3c
4c
5c     jsubn(x),ysubn(x)      spherical bessel functions
6c     do not use x=0.
7c     jsub0(0)=1.   jsub n (0)=0.   for all n.gt.0
8c     ysub n (0)= infinity  for all n.ge.0.
9c     y sub n (x)=-j sub (-(n+1))  (x)        for all n.ge.0
10c     mn  maximum number of subscript of bessel functions computed
11c     nout   maximum subscript of desired bessel functions
12c     recurrence techniques for bessel functions
13c     fr.  mechel   math. of comp.  vol. 22   no. 101  jan. 1968
14c
15c
16c
17      implicit none
18      integer mn,nout,nmin,nmax,ix,n,twon1,nm,nx,mnout
19      double precision j(0:1),y(0:1),jx,jnmp,jnm,jnm1,x,
20     &     ynm1,ynm,ynmp,x2,pmin1,qmin1,p0,q0,b0,p1,q1,gnmax1,c
21c
22c     the dimension of j and y may be varied to suit the users needs
23c     set dimension in calling routine  j((0,nout+5)),y((0,nout+5))
24c
25c
26c
27      double precision bound
28      integer nadd,nxadd
29      data  bound /1.e13/, nadd /10/,  nxadd /10/
30c
31c     bound  nadd   nxadd   may be varied
32c     nxadd controls the number of terms in the continued fraction
33c
34      mn=1
35      j(0)=sin(x)/x
36      y(0)=-cos(x)/x
37      y(1)=y(0)/x-j(0)
38      j(1)=j(0)/x+y(0)
39      nmin=1
40      nmax= nout
41c
42c     upward recursion valid for n<x for both bessel functions
43c
44      ix=x
45      nx=min0(nout,ix)
46      do 66 n=1,nx
47         mn=n+1
48         twon1=n+mn
49         j(mn)=twon1*j(n)/x-j(n-1)
50         y(mn)=twon1*y(n)/x-y(n-1)
51 66   continue
52      if(nout.le.mn) return
53      jx=j(mn)
54      mnout=mn
55      nmin=mn
56c
57c     upward recursion for ysubn(x) to obtain good guess for jsubn(x)
58c
59      do 1 n=nmin,nout
60         mn=n+1
61         twon1=n+mn
62         y(mn)=twon1*y(n)/x-y(n-1)
63 1    continue
64      ynm1=y(nout-1)
65      ynm=y(nout)
66      nmin=nout
67      nmax=nout+nadd
68 33   do 11 n=nmin,nmax
69         mn=n+1
70         twon1=n+mn
71         ynmp=twon1*ynm/x-ynm1
72         ynm1=ynm
73         ynm=ynmp
74 11   continue
75      if(abs(ynmp).gt.bound) go to 2
76 22   nmin=nmax+1
77      nmax=nmax+nadd
78      go to 33
79    2 x2=x*x
80      jnm=-1.0/(x2*ynmp)
81      pmin1=1.
82      qmin1=0.
83      q0=1.0
84c
85c     develop continued fraction to obtain good guess for jsubn+1(x)
86c
87      b0=mn+mn+1
88      p0=mn+mn+1
89      do 4 n=1,nxadd
90         b0=b0+2.
91         p1=b0*p0-x2*pmin1
92         q1=b0*q0-x2*qmin1
93         pmin1=p0
94         p0=p1
95         qmin1=q0
96         q0=q1
97    4 continue
98c
99      gnmax1=p1/q1
100      jnmp=x*jnm/gnmax1
101      do 55 n=mn-1,nout+1,-1
102         nm=n+1
103         twon1=n+nm
104         jnm1=twon1*jnm/x-jnmp
105         jnmp=jnm
106         jnm=jnm1
107 55   continue
108      j(nout+1)=jnmp
109      j(nout)=jnm1
110      do 5 n=nout,mnout+1,-1
111         nm=n+1
112         twon1=n+nm
113        j(n-1)=twon1*j(n)/x-j(n+1)
114 5    continue
115      c=jx/j(mnout)
116c
117c     correct bessel functions with weighting factor
118c
119      if (x.lt.1.e-3) return
120      do 6 n=mnout,nout
121         j(n)=c*j(n)
122 6    continue
123c
124      return
125      end
126c $Id$
127