function E = irrad(W,r)

n=length(r);
u = zeros(size(r));
for k=1:n
u(k) = quad(@mhankel,0,1,[],[],W,r(k));
end
E = 4*u.*conj(u);

Published with MATLAB® 7.13