12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576 |
- function [FE] = calc_emntd_bin(xno, yno, cx, cy, dx, dy, eps)
- nx = length(cx);
- ny = length(cy);
- if (length(cx)~=length(dx)) || (length(cy)~=length(dy)) || (length(eps)~=(nx*ny))
- error("incorrect binary grating definition");
- end
- FE = cellmat(1,2,2*yno-1,2*xno-1);
- [CX,CY] = meshgrid(cx,cy);
- [DX,DY] = meshgrid(dx,dy);
- ix = linspace(1,xno-1,xno-1);
- iy = linspace(1,yno-1,yno-1);
- [IX,IY] = meshgrid(ix,iy);
-
- for ip = 1:nx*ny
- fx = (sin(ix*pi*DX(ip))./(pi*ix)).*exp((-2*pi*1i*CX(ip))*ix);
- fy = (sin(iy*pi*DY(ip))./(pi*iy)).*exp((-2*pi*1i*CY(ip))*iy);
- FX = (sin(IX*pi*DX(ip))./(pi*IX)).*exp((-2*pi*1i*CX(ip))*IX);
- FY = (sin(IY*pi*DY(ip))./(pi*IY)).*exp((-2*pi*1i*CY(ip))*IY);
- M = zeros(2*yno-1,2*xno-1);
- M(yno+1:2*yno-1,xno) = DX(ip)*fy;
- M(yno-1:-1:1,xno) = conj(M(yno+1:2*yno-1,xno));
- M(yno,xno+1:2*xno-1) = DY(ip)*fx;
- M(yno,xno-1:-1:1) = conj(M(yno,xno+1:2*xno-1));
- M(yno+1:2*yno-1,xno+1:2*xno-1) = FX.*FY;
- M(yno+1:2*yno-1,xno-1:-1:1) = conj(FX).*FY;
- M(yno-1:-1:1,xno+1:2*xno-1) = FX.*conj(FY);
- M(yno-1:-1:1,xno-1:-1:1) = conj(FX.*FY);
- FE{1,1} = FE{1,1} + eps(ip)*M;
- FE{1,2} = FE{1,2} + (1/eps(ip))*M;
- FE{1,1}(yno,xno) = FE{1,1}(yno,xno) + DX(ip)*DY(ip)*eps(ip);
- FE{1,2}(yno,xno) = FE{1,2}(yno,xno) + DX(ip)*DY(ip)/eps(ip);
- end
- end
|