%% Authored by Joshua Soneson 2007
function[M,J,J_,dz,dr,z,r]=computational_grid(Z,R,G,a,d,gamma,N)
% generates node vectors for discretization in axial and radial directions
ppw_ax = 40; % est. points per wavelength in axial direction
M = round(ppw_ax*Z*G); % number of meshpoints in axial direction
dz = Z/(M-1); % axial stepsize
z = [0:dz:Z]; % axial node vector
ppw_rad = 50; % points per wavelength in radial direction
J = round(ppw_rad*R*G/pi); % number of meshpoints in radial direction
dr = (R+0.25)/(J-1); % [0,R] is physical; the extra 0.25 is for PML
r = [0:dr:R+0.25]'; % radial node vector
J_ = ceil(J*R/(R+0.25)); % number of nodes in [0,R]
% print parameters:
fprintf('\tdr = %2.3f mm\tJ = %d\n',10*a*dr,J)
fprintf('\tdz = %2.3f mm\tM = %d\n\n',10*d*dz,M)