I have a matrix H with dimensions N x M. I also have a few (1 x N) vectors denoted here as n1 and n2, and a few (1 x M) vectors denoted here as m1 and m2.
I want to define H such that H_(n,m) is a function of n1(n), n2(n), m1(m), and m2(m) without a for loop. Below is my attempt with arrayfun:
function [H] = solveH(N, M) % Initialize n1, n2, m1, and m2
n1 = <some code here>; n2 = <some code here>; m1 = <some code here>; m2 = <some code here>; % NMat(i) and MMat(i) contain all subscript pairs for H
[MMat, NMat] = meshgrid(m, n); H = arrayfun(@calcH, NMat, MMat); % inputs n and m are scalars
function outputH = calcH(n, m) outputH = <some function of n1(n), n2(n), m1(m), and m2(m)>; endend
This works fine, except that it is actually slower (!) than using a for loop because n1, n2, m1, and m2 have a scope in multiple functions and are slow like global variables.
I also considered turning n1, n2, m1, and m2 into matrices and doing element-by-element operations, but in reality I have something like five (1 x N) vectors and five (1 x M) vectors, and using repmat that many times wastes memory and makes it even slower.
What approach should I use so that my program executes the fastest?
EDIT: Here is most of the function for those who want to see exactly what operations I'm using.
% Solves the bold matrices H and E analytically
% ka wavenumber in air (1 x M vector)
% ks wavenumber in bar (1 x M vector)
% a air gap of HCG (scalar)
% s bar width of HCG (scalar)
% nBar refractive index (scalar)
% period period of HCG (scalar)
% lambda wavelength of light (scalar)
% M orders - region II (scalar)
% N orders - regions I, III (scalar)
function [H, E, beta] = solveHE(ka, ks, a, s, nBar, period, lambda, N, M, polarization) n = 1:1:N; m = 1:1:M; % Precalculate constants
gamma = conj( 2*pi * sqrt(1/lambda^2 - (n-1).^2/period^2) ); % (1 x N)
beta = conj( sqrt((2*pi*nBar/lambda)^2 - ks.^2) ); % (1 x M vector)
p = (2*pi/period) * (n-1); % (1 x N vector)
aa = a/2; ss = s/2; dif = period - aa; nInv = (1 / nBar^2); % NMat(i) and MMat(i) contain all subscript pairs for H and E
[MMat, NMat] = meshgrid(m, n); % Efficient vectorization?
[H, E] = arrayfun(@calcHE, NMat, MMat); % Calculate H_(n,m) and E_(n,m)
function [resH, resE] = calcHE(n, m) d = ~(n-1); % delta = 1 for n = 0 (index 1), 0 otherwise
d = 1 / period * (2-d); % precalculate
hAir = d*cos(ks(m)*ss) * 2 * eval(ka(m),aa,p(n),aa); hBar = d*cos(ka(m)*aa) * (eval(ks(m),ss,p(n),dif) - eval(ks(m),-ss,p(n),aa)); resH = hAir + hBar; if (strcmp(polarization, 'TM')) resE = beta(m)/gamma(n) * (hAir + nInv*hBar); else resE = gamma(n)/beta(m) * resH; end end % All parameters are scalars
%
% sin(KU - PA) sin(KU + PA)
% ------------ + ------------
% 2(K-P) 2(K+P)
function z = eval(K, U, P, A) z = sin(K*U-P*A)/(2*(K-P)) + sin(K*U+P*A)/(2*(K+P)); endend
Best Answer