function z = BesselJPrimeRoots(M,P)
%z = BesselJPrimeRoots(m,p) returns the root where the m-th
%  Bessel function of the First Kind is either maximum of minimum.
%  That is, d[ J_m(z) ] / dz =0 with z_{mp} being the p-th root.
%  This program is useful in various electromagnetics problems 
%  where high order modes, such as TE_{22,6} are employed.
%  
%INPUTS:
%  'm': the order of the Bessel function.
%  'p': the p-th root where the Bessel function slope is zero.
%    Note: either 'm' or 'p' or both can be vectors, but the computation
%          time will be that much longer.
%
%OUTPUT:
%  'z': the root z_{mp}.  z=[] if p==0. 'z' is of size Nm x Np, where
%       Nm = length(m) and Np = length(p).
%
%NOTES: Accurate to 5 significant figures if the internal variable
% is dx=8e-5;  Accuracy can be extended to six sig figs with 
% dx=1e-6, but computation time is longer.
%   Requires the function "find_intersection.m", which can be 
% downloaded from 
%   http://web.mit.edu/cjoye/www/research/matlab/find_intersection.m
%
% Last edit by Colin Joye, 2/18/04 ; 8/6/04 (added vector capability).

% ------- step size (and ultimate accuracy) ------------
dx = 8e-6; 
% ------- --------------------------------- ------------

for h = 1:length(M),
  m = M(h);
for k = 1:length(P),
  p = P(k);

% Define range of interest to get certain root:
%x  = linspace( m+5*(p-1) , m+7*p ,1/dx); 
x  = linspace( 0 , m+7*p ,1/dx); 

% generate besselj function 
y  = besselj(m,x); 

% find derivative
dy_dx = diff(y)/dx;

% find intersection of discrete slope with zero line:
[xi,yi] = find_intersection(x(1:end-1)+dx/2, dy_dx, zeros(1,length(dy_dx)));

% assign output
if p~=0,
  z(h,k) = xi(p);
else
  z(h,k)=[];
end

end % P loop
end % M loop
% --------- End BesselJPrimeRoots.m --------------