mrtd2dcode.m

MRTD 2D code in matlab

% matrix indices preprocessing
ind = 1:(xmax*ymax);
ind = reshape(ind,xmax,ymax);
...
for t = 1:Nt, % time-marching cycle
% indices initialization
ind1 = [ind(1:xmax-1,:) ; -ind(xmax-1,:)];
ind2 = [ind(1,:) ; -ind(1:xmax-1,:)];
ind3 = [ind(:,1:ymax-1) -ind(:,ymax-1)];
ind4 = [ind(:,1) -ind(:,1:ymax-1)];
for i = 1:na, % alpha cycle
Ez = Ez + C1(i).*( Hy(abs(ind1)) - Hy(abs(ind2)) ...
- Hx(abs(ind3)) + Hx(abs(ind4)) );
ind1 = ind1 + (ind1>0);
mask = ismember(abs(ind1),PEC);
ind1(mask) = -ind1(mask);
ind1 = ind1 + (ind1<0);
ind2 = ind2 + (ind2>0);
mask = ismember(abs(ind2),PEC);
ind2(mask) = -ind2(mask);
ind2 = ind2 + (ind2<0);
ind3 = ind3 + xmax.*(ind3>0);
mask = ismember(abs(ind3),PEC);
ind3(mask) = -ind3(mask);
ind3 = ind3 + xmax.*(ind3<0);
ind4 = ind4 + xmax.*(ind4>0);
mask = ismember(abs(ind4),PEC);
ind4(mask) = -ind4(mask);
ind4 = ind4 + xmax.*(ind4<0);
end;
Ez(PEC) = 0; % PEC boundary condition
...
% indices initialization
ind1 = [ind(2:xmax-1,:) ; -ind([xmax xmax],:)];
ind2 = [ind(1,:) ; -ind(2:xmax,:)];
ind3 = [ind(:,2:ymax-1) -ind(:,[ymax ymax])];
ind4 = [ind(:,1) -ind(:,2:ymax)];
for i = 1:na; % alpha cycle
Hy = Hy + C2(i).*( sign(ind1).*Ez(abs(ind1)) ...
+ sign(ind2).*Ez(abs(ind2)) );
Hx = Hx - C2(i).*( sign(ind3).*Ez(abs(ind3)) ...
+ sign(ind4).*Ez(abs(ind4)) );
ind1 = ind1 + 1;
mask = ismember(abs(ind1),PEC);
ind1(mask) = -ind1(mask);
ind2 = ind2 + 1;
mask = ismember(abs(ind2),PEC);
ind2(mask) = -ind2(mask);
ind3 = ind3 + xmax;
mask = ismember(abs(ind3),PEC);
ind3(mask) = -ind3(mask);
ind4 = ind4 + xmax;
mask = ismember(abs(ind4),PEC);
ind4(mask) = -ind4(mask);
end;
end;