Disk Harmonic Functions

Computes Disk Harmonic Functions; uses Bessel functions and derivatives
Requires BessDerivZerosBisect2.m
%function dnm = DiskHarmonic(m,n,r,theta)
% if(m == 0)
% dnm = 1*Rnm
% elseif(m < 0)
% dnm = sqrt(2)*Rnm.*sin(mu*theta)
% elseif(m > 0)
% dnm = sqrt(2)*Rnm.*cos(mu*theta)
% end
% m<0: dnm(r,theta) = sqrt(2)*Rnm(r)*sin(m*theta) (mu = |m|)
% m=0: dnm(r,theta) = 1*R_n0(r) (mu = |m|)
% m>0: dnm(r,theta) = sqrt(2)*Rnm(r)*cos(m*theta) (m = |mu|)
%
% n = 0,1,2,...
% m = ...-2,-1,0,1,2... (Does not appear to be a requirement the \m\<=n
% r & theta are vectors of polar coordinates of individual points (not a grid)
% Rnm(r) = anm*J_m(knm*r)
% where
% J_m(r) is a Bessel function of the first kind
% knm = a spatial freq. = the nth positive root of d/dr(J_m(r)) = 0
% anm = 1/sqrt((1-(m/knm)^2)*(J_m(knm))) for n>0 or m>0 (a_00=1)
% For n=0, m=0: d_00=1 (a00=1, J0(0)=1, R00=1) (piston term)