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