subroutine zeros(f, z, nnm, nnz, id) * * Find the Zeros of a J and J' * Inputs: * ------- * f: Bessel function J or its derivative J' * nnm, nnz: Dimensions of the array holding the zeros * id: 0 if the function f is J * 1 if the function f is J' * * Outputs: * -------- * z() : An array of dimensions 0:nnm, nnz holding the zeros * * Subroutines: * ------------ * zbrent: A root finder from page 253 of Numerical Recipes * by Press, Flannery, Teukolsky, and Vetterling that * finds the root known to lie between two numbers. * * bessj: (Passed as f) is the bessel function J from the * above reference. J' is computed using the recurrence * relations (9.1.27) on page 361 of Handbook of * Mathematical Functions by M. Abramowitz and I. A. Stegun. * real z(0:nnm, nnz) external f * * Stepping increment * del = 0.1 * * Tolerance for the zero finder * tol = 1e-9 xfirst = del do no = 0, nnm n = no - 1 nn = 1 if(n.eq.1.and.id.eq.1)then x = del else x = xfirst endif * * Bracket the zero (check for zero transition) * 1 x1 = x x2 = x + del f1 = f(n,x1) f2 = f(n,x2) x = x + del if(f1*f2.gt.0)goto 1 * * Call the root finder from Numerical Recipes * an extra argument "n" has been added for the Bessel function. * root = zbrent (f,x1,x2,tol,n) if(nn.eq.1)xfirst = x1 if(nn.le.nnz) then z(no, nn) = root else goto 2 endif * * Set the next starting value so we don't bump * into the same zero again. * x = x2 * * Increment the zero index * nn = nn + 1 goto 1 2 continue enddo end