Vectorization - Sum and Bessel function

Question

Can anyone help vectorize this Matlab code? The specific problem is the sum and bessel function with vector inputs. Thank you!

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];
for ii = 1:length(rho_g)
    for jj = 1:length(phi_g)
        % Coordinates
        rho_o = rho_g(ii);
        phi_o = phi_g(jj);
        % factors
        fc = cos(n.*(phi_o-phi_s));
        fs = sin(n.*(phi_o-phi_s));

        Ez_t(ii,jj) = sum(tau.*besselj(n,k(3)*rho_s).*besselh(n,2,k(3)*rho_o).*fc);
    end
end

Answer 1

Initialization -

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

Nested loops form (Copy from your code and shown here for comparison only) -

for ii = 1:length(rho_g)
    for jj = 1:length(phi_g)
        % Coordinates
        rho_o = rho_g(ii);
        phi_o = phi_g(jj);
        % factors
        fc = cos(n.*(phi_o-phi_s));
        fs = sin(n.*(phi_o-phi_s));

        Ez_t(ii,jj) = sum(tau.*besselj(n,k(3)*rho_s).*besselh(n,2,k(3)*rho_o).*fc);
    end
end

Vectorized solution -

%%// Term - 1
term1 = repmat(tau.*besselj(n,k(3)*rho_s),[N*N 1]);

%%// Term - 2
[n1,rho_g1] = meshgrid(n,rho_g);
term2_intm = besselh(n1,2,k(3)*rho_g1);
term2 = transpose(reshape(repmat(transpose(term2_intm),[N 1]),N,N*N));

%%// Term -3
angle1 = repmat(bsxfun(@times,bsxfun(@minus,phi_g,phi_s')',n),[N 1]);
fc = cos(angle1);

%%// Output
Ez_t = sum(term1.*term2.*fc,2);
Ez_t = transpose(reshape(Ez_t,N,N));

Points to note about this vectorization or code simplification –

1. ‘fs’ doesn’t change the output of the script, Ez_t, so it could be removed for now.
2. The output seems to be ‘Ez_t’,which requires three basic terms in the code as – tau.*besselj(n,k(3)*rho_s), besselh(n,2,k(3)*rho_o) and fc. These are calculated separately for vectorization as terms1,2 and 3 respectively.
3. All these three terms appear to be of 1xN sizes. Our aim thus becomes to calculate these three terms without loops. Now, the two loops run for N times each, thus giving us a total loop count of NxN. Thus, we must have NxN times the data in each such term as compared to when these terms were inside the nested loops.
4. This is basically the essence of the vectorization done here, as the three terms are represented by ‘term1’,’term2’ and ‘fc’ itself.