postprocess.m

强大的计算电磁场本征函数与本征模的程序

ezz4 = ones(1,nx*ny);   ezz4(:) = epszz(2:nx+1,2:ny+1);

bzxne = (1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*eyx4./ezz4./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy3.*eyy1.*w.*eyy2+1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*(1-exx4./ezz4)./ezz3./ezz2./(w.*exx3+e.*exx2)./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*exx1.*s)./b;

bzxse = (-1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*eyx3./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy1.*w.*eyy2+1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*(1-exx3./ezz3)./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*n.*exx1.*exx4)./b;

bzxnw = (-1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*eyx1./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy2.*e-1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*(1-exx1./ezz1)./ezz3./ezz2./(w.*exx3+e.*exx2)./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*exx4.*s)./b;

bzxsw = (1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*eyx2./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*e-1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*(1-exx2./ezz2)./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx3.*n.*exx1.*exx4)./b;

bzxn = ((1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*n.*ezz1.*ezz2./eyy1.*(2.*eyy1./ezz1./n.^2+eyx1./ezz1./n./w)+1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*n.*ezz4.*ezz3./eyy4.*(2.*eyy4./ezz4./n.^2-eyx4./ezz4./n./e))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+((ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*(1./2.*ezz4.*((1-exx1./ezz1)./n./w-exy1./ezz1.*(2./n.^2-2./n.^2.*s./(n+s)))./exx1.*ezz1.*w+(ezz4-ezz1).*s./n./(n+s)+1./2.*ezz1.*(-(1-exx4./ezz4)./n./e-exy4./ezz4.*(2./n.^2-2./n.^2.*s./(n+s)))./exx4.*ezz4.*e)-(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*(-ezz3.*exy2./n./(n+s)./exx2.*w+(ezz3-ezz2).*s./n./(n+s)-ezz2.*exy3./n./(n+s)./exx3.*e))./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzxs =((1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*s.*ezz2.*ezz1./eyy2.*(2.*eyy2./ezz2./s.^2-eyx2./ezz2./s./w)+1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*s.*ezz3.*ezz4./eyy3.*(2.*eyy3./ezz3./s.^2+eyx3./ezz3./s./e))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+((ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*(-ezz4.*exy1./s./(n+s)./exx1.*w-(ezz4-ezz1).*n./s./(n+s)-ezz1.*exy4./s./(n+s)./exx4.*e)-(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*(1./2.*ezz3.*(-(1-exx2./ezz2)./s./w-exy2./ezz2.*(2./s.^2-2./s.^2.*n./(n+s)))./exx2.*ezz2.*w-(ezz3-ezz2).*n./s./(n+s)+1./2.*ezz2.*((1-exx3./ezz3)./s./e-exy3./ezz3.*(2./s.^2-2./s.^2.*n./(n+s)))./exx3.*ezz3.*e))./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzxe = ((n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*(1./2.*n.*ezz4.*ezz3./eyy4.*(2./e.^2-eyx4./ezz4./n./e)+1./2.*s.*ezz3.*ezz4./eyy3.*(2./e.^2+eyx3./ezz3./s./e))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+(-1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*ezz1.*(1-exx4./ezz4)./n./exx4.*ezz4-1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*ezz2.*(1-exx3./ezz3)./s./exx3.*ezz3)./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzxw = ((-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*(1./2.*n.*ezz1.*ezz2./eyy1.*(2./w.^2+eyx1./ezz1./n./w)+1./2.*s.*ezz2.*ezz1./eyy2.*(2./w.^2-eyx2./ezz2./s./w))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+(1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*ezz4.*(1-exx1./ezz1)./n./exx1.*ezz1+1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*ezz3.*(1-exx2./ezz2)./s./exx2.*ezz2)./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzxp = (((-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*(1./2.*n.*ezz1.*ezz2./eyy1.*(-2./w.^2-2.*eyy1./ezz1./n.^2+k.^2.*eyy1-eyx1./ezz1./n./w)+1./2.*s.*ezz2.*ezz1./eyy2.*(-2./w.^2-2.*eyy2./ezz2./s.^2+k.^2.*eyy2+eyx2./ezz2./s./w))+(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*(1./2.*n.*ezz4.*ezz3./eyy4.*(-2./e.^2-2.*eyy4./ezz4./n.^2+k.^2.*eyy4+eyx4./ezz4./n./e)+1./2.*s.*ezz3.*ezz4./eyy3.*(-2./e.^2-2.*eyy3./ezz3./s.^2+k.^2.*eyy3-eyx3./ezz3./s./e)))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+((ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*(1./2.*ezz4.*(-k.^2.*exy1-(1-exx1./ezz1)./n./w-exy1./ezz1.*(-2./n.^2-2./n.^2.*(n-s)./s))./exx1.*ezz1.*w+(ezz4-ezz1).*(n-s)./n./s+1./2.*ezz1.*(-k.^2.*exy4+(1-exx4./ezz4)./n./e-exy4./ezz4.*(-2./n.^2-2./n.^2.*(n-s)./s))./exx4.*ezz4.*e)-(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*(1./2.*ezz3.*(-k.^2.*exy2+(1-exx2./ezz2)./s./w-exy2./ezz2.*(-2./s.^2+2./s.^2.*(n-s)./n))./exx2.*ezz2.*w+(ezz3-ezz2).*(n-s)./n./s+1./2.*ezz2.*(-k.^2.*exy3-(1-exx3./ezz3)./s./e-exy3./ezz3.*(-2./s.^2+2./s.^2.*(n-s)./n))./exx3.*ezz3.*e))./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzyne = (1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*(1-eyy4./ezz4)./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy3.*eyy1.*w.*eyy2+1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*exy4./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*exx1.*s)./b;

bzyse = (-1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*(1-eyy3./ezz3)./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy1.*w.*eyy2+1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*exy3./ezz3./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*n.*exx1.*exx4)./b;

bzynw = (-1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*(1-eyy1./ezz1)./ezz4./ezz3./(n.*eyy3+s.*eyy4)./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy2.*e-1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*exy1./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*exx4.*s)./b;

bzysw = (1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*(1-eyy2./ezz2)./ezz4./ezz3./(n.*eyy3+s.*eyy4)./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*e-1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*exy2./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx3.*n.*exx1.*exx4)./b;

bzyn = ((1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*ezz1.*ezz2./eyy1.*(1-eyy1./ezz1)./w-1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*ezz4.*ezz3./eyy4.*(1-eyy4./ezz4)./e)./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*(1./2.*ezz4.*(2./n.^2+exy1./ezz1./n./w)./exx1.*ezz1.*w+1./2.*ezz1.*(2./n.^2-exy4./ezz4./n./e)./exx4.*ezz4.*e)./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzys = ((-1./2.*(-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*ezz2.*ezz1./eyy2.*(1-eyy2./ezz2)./w+1./2.*(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*ezz3.*ezz4./eyy3.*(1-eyy3./ezz3)./e)./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e-(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*(1./2.*ezz3.*(2./s.^2-exy2./ezz2./s./w)./exx2.*ezz2.*w+1./2.*ezz2.*(2./s.^2+exy3./ezz3./s./e)./exx3.*ezz3.*e)./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzye = (((-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*(-n.*ezz2./eyy1.*eyx1./e./(e+w)+(ezz1-ezz2).*w./e./(e+w)-s.*ezz1./eyy2.*eyx2./e./(e+w))+(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*(1./2.*n.*ezz4.*ezz3./eyy4.*(-(1-eyy4./ezz4)./n./e-eyx4./ezz4.*(2./e.^2-2./e.^2.*w./(e+w)))+1./2.*s.*ezz3.*ezz4./eyy3.*((1-eyy3./ezz3)./s./e-eyx3./ezz3.*(2./e.^2-2./e.^2.*w./(e+w)))+(ezz4-ezz3).*w./e./(e+w)))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+(1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*ezz1.*(2.*exx4./ezz4./e.^2-exy4./ezz4./n./e)./exx4.*ezz4.*e-1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*ezz2.*(2.*exx3./ezz3./e.^2+exy3./ezz3./s./e)./exx3.*ezz3.*e)./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzyw = (((-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*(1./2.*n.*ezz1.*ezz2./eyy1.*((1-eyy1./ezz1)./n./w-eyx1./ezz1.*(2./w.^2-2./w.^2.*e./(e+w)))-(ezz1-ezz2).*e./w./(e+w)+1./2.*s.*ezz2.*ezz1./eyy2.*(-(1-eyy2./ezz2)./s./w-eyx2./ezz2.*(2./w.^2-2./w.^2.*e./(e+w))))+(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*(-n.*ezz3./eyy4.*eyx4./w./(e+w)-s.*ezz4./eyy3.*eyx3./w./(e+w)-(ezz4-ezz3).*e./w./(e+w)))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+(1./2.*(ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*ezz4.*(2.*exx1./ezz1./w.^2+exy1./ezz1./n./w)./exx1.*ezz1.*w-1./2.*(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*ezz3.*(2.*exx2./ezz2./w.^2-exy2./ezz2./s./w)./exx2.*ezz2.*w)./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

bzyp = (((-n.*ezz4.*ezz3./eyy4-s.*ezz3.*ezz4./eyy3).*(1./2.*n.*ezz1.*ezz2./eyy1.*(-k.^2.*eyx1-(1-eyy1./ezz1)./n./w-eyx1./ezz1.*(-2./w.^2+2./w.^2.*(e-w)./e))+(ezz1-ezz2).*(e-w)./e./w+1./2.*s.*ezz2.*ezz1./eyy2.*(-k.^2.*eyx2+(1-eyy2./ezz2)./s./w-eyx2./ezz2.*(-2./w.^2+2./w.^2.*(e-w)./e)))+(n.*ezz1.*ezz2./eyy1+s.*ezz2.*ezz1./eyy2).*(1./2.*n.*ezz4.*ezz3./eyy4.*(-k.^2.*eyx4+(1-eyy4./ezz4)./n./e-eyx4./ezz4.*(-2./e.^2-2./e.^2.*(e-w)./w))+1./2.*s.*ezz3.*ezz4./eyy3.*(-k.^2.*eyx3-(1-eyy3./ezz3)./s./e-eyx3./ezz3.*(-2./e.^2-2./e.^2.*(e-w)./w))+(ezz4-ezz3).*(e-w)./e./w))./ezz4./ezz3./(n.*eyy3+s.*eyy4)./ezz2./ezz1./(n.*eyy2+s.*eyy1)./(e+w).*eyy4.*eyy3.*eyy1.*w.*eyy2.*e+((ezz3./exx2.*ezz2.*w+ezz2./exx3.*ezz3.*e).*(1./2.*ezz4.*(-2./n.^2-2.*exx1./ezz1./w.^2+k.^2.*exx1-exy1./ezz1./n./w)./exx1.*ezz1.*w+1./2.*ezz1.*(-2./n.^2-2.*exx4./ezz4./e.^2+k.^2.*exx4+exy4./ezz4./n./e)./exx4.*ezz4.*e)-(ezz4./exx1.*ezz1.*w+ezz1./exx4.*ezz4.*e).*(1./2.*ezz3.*(-2./s.^2-2.*exx2./ezz2./w.^2+k.^2.*exx2+exy2./ezz2./s./w)./exx2.*ezz2.*w+1./2.*ezz2.*(-2./s.^2-2.*exx3./ezz3./e.^2+k.^2.*exx3-exy3./ezz3./s./e)./exx3.*ezz3.*e))./ezz3./ezz2./(w.*exx3+e.*exx2)./ezz4./ezz1./(w.*exx4+e.*exx1)./(n+s).*exx2.*exx3.*n.*exx1.*exx4.*s)./b;

ii = zeros(nx,ny);
ii(:) = (1:nx*ny); 

% NORTH boundary

ib = zeros(nx,1);  ib(:) = ii(1:nx,ny);

switch (boundary(1))
  case 'S',   sign = +1;
  case 'A',   sign = -1;
  case '0',   sign = 0;
  otherwise,  
    error('Unrecognized north boundary condition: %s.\n', boundary(1));
end

bzxs(ib)  = bzxs(ib)  + sign*bzxn(ib);
bzxse(ib) = bzxse(ib) + sign*bzxne(ib);
bzxsw(ib) = bzxsw(ib) + sign*bzxnw(ib);
bzys(ib)  = bzys(ib)  - sign*bzyn(ib);
bzyse(ib) = bzyse(ib) - sign*bzyne(ib);
bzysw(ib) = bzysw(ib) - sign*bzynw(ib);

% SOUTH boundary

ib = zeros(nx,1);  ib(:) = ii(1:nx,1);

switch (boundary(2))
  case 'S',   sign = +1;
  case 'A',   sign = -1;
  case '0',   sign = 0;
  otherwise,  
    error('Unrecognized south boundary condition: %s.\n', boundary(2));
end

bzxn(ib)  = bzxn(ib)  + sign*bzxs(ib);
bzxne(ib) = bzxne(ib) + sign*bzxse(ib);
bzxnw(ib) = bzxnw(ib) + sign*bzxsw(ib);
bzyn(ib)  = bzyn(ib)  - sign*bzys(ib);
bzyne(ib) = bzyne(ib) - sign*bzyse(ib);
bzynw(ib) = bzynw(ib) - sign*bzysw(ib);

% EAST boundary

ib = zeros(1,ny);  ib(:) = ii(nx,1:ny);

switch (boundary(3))
  case 'S',   sign = +1;
  case 'A',   sign = -1;
  case '0',   sign = 0;
  otherwise,  
    error('Unrecognized east boundary condition: %s.\n', boundary(3));
end

bzxw(ib)  = bzxw(ib)  + sign*bzxe(ib);
bzxnw(ib) = bzxnw(ib) + sign*bzxne(ib);
bzxsw(ib) = bzxsw(ib) + sign*bzxse(ib);
bzyw(ib)  = bzyw(ib)  - sign*bzye(ib);
bzynw(ib) = bzynw(ib) - sign*bzyne(ib);
bzysw(ib) = bzysw(ib) - sign*bzyse(ib);

% WEST boundary

ib = zeros(1,ny);  ib(:) = ii(1,1:ny);

switch (boundary(4))
  case 'S',   sign = +1;
  case 'A',   sign = -1;
  case '0',   sign = 0;
  otherwise,  
    error('Unrecognized west boundary condition: %s.\n', boundary(4));
end

bzxe(ib)  = bzxe(ib)  + sign*bzxw(ib);
bzxne(ib) = bzxne(ib) + sign*bzxnw(ib);
bzxse(ib) = bzxse(ib) + sign*bzxsw(ib);
bzye(ib)  = bzye(ib)  - sign*bzyw(ib);
bzyne(ib) = bzyne(ib) - sign*bzynw(ib);
bzyse(ib) = bzyse(ib) - sign*bzysw(ib);

% Assemble sparse matrix

iall = zeros(1,nx*ny);          iall(:) = ii;
is = zeros(1,nx*(ny-1));        is(:) = ii(1:nx,1:(ny-1));
in = zeros(1,nx*(ny-1));        in(:) = ii(1:nx,2:ny);
ie = zeros(1,(nx-1)*ny);        ie(:) = ii(2:nx,1:ny);
iw = zeros(1,(nx-1)*ny);        iw(:) = ii(1:(nx-1),1:ny);
ine = zeros(1,(nx-1)*(ny-1));   ine(:) = ii(2:nx, 2:ny);
ise = zeros(1,(nx-1)*(ny-1));   ise(:) = ii(2:nx, 1:(ny-1));
isw = zeros(1,(nx-1)*(ny-1));   isw(:) = ii(1:(nx-1), 1:(ny-1));
inw = zeros(1,(nx-1)*(ny-1));   inw(:) = ii(1:(nx-1), 2:ny);

Bzx = sparse ([iall,iw,ie,is,in,ine,ise,isw,inw], ...
	[iall,ie,iw,in,is,isw,inw,ine,ise], ...
	[bzxp(iall),bzxe(iw),bzxw(ie),bzxn(is),bzxs(in), ...
     bzxsw(ine),bzxnw(ise),bzxne(isw),bzxse(inw)]);

Bzy = sparse ([iall,iw,ie,is,in,ine,ise,isw,inw], ...
	[iall,ie,iw,in,is,isw,inw,ine,ise], ...
	[bzyp(iall),bzye(iw),bzyw(ie),bzyn(is),bzys(in), ...
     bzysw(ine),bzynw(ise),bzyne(isw),bzyse(inw)]);

B = [Bzx Bzy];

Hz = zeros(size(Hx));
Hz(:) = B*reshape([Hx,Hy],2*nx*ny,1)/j;

nx = nx-1;
ny = ny-1;

exx = epsxx(2:nx+1,2:ny+1);
exy = epsxy(2:nx+1,2:ny+1);
eyx = epsyx(2:nx+1,2:ny+1);
eyy = epsyy(2:nx+1,2:ny+1);
ezz = epszz(2:nx+1,2:ny+1);
edet = (exx.*eyy - exy.*eyx);

h = dx(2:nx+1)*ones(1,ny);
v = ones(nx,1)*dy(2:ny+1);

i1 = ii(1:nx,2:ny+1);
i2 = ii(1:nx,1:ny);
i3 = ii(2:nx+1,1:ny);
i4 = ii(2:nx+1,2:ny+1);

Dx = +neff*(Hy(i1) + Hy(i2) + Hy(i3) + Hy(i4))/4 + ...
     (Hz(i1) + Hz(i4) - Hz(i2) - Hz(i3))./(j*2*k*v);
Dy = -neff*(Hx(i1) + Hx(i2) + Hx(i3) + Hx(i4))/4 - ...
     (Hz(i3) + Hz(i4) - Hz(i1) - Hz(i2))./(j*2*k*h);
Dz = ((Hy(i3) + Hy(i4) - Hy(i1) - Hy(i2))./(2*h) - ...
      (Hx(i1) + Hx(i4) - Hx(i2) - Hx(i3))./(2*v))/(j*k);

Ex = (eyy.*Dx - exy.*Dy)./edet;
Ey = (exx.*Dy - eyx.*Dx)./edet;
Ez = Dz./ezz;