postprocess.m

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

function [Hz,Ex,Ey,Ez] = postprocess (lambda,neff,Hx,Hy,dx,dy,varargin);

% This function takes the two computed transverse magnetic
% fields (Hx and Hy) of an optical waveguide structure and
% solves for the remaining 4 vield components:  Hz, Ex, Ey,
% and Ez.
%
% USAGE:
% 
% [Hz,Ex,Ey,Ez] = postprocess(lambda, neff, Hx, Hy, dx, dy, ...
%                     eps, boundary);
% [Hz,Ex,Ey,Ez] = postprocess(lambda, neff, Hx, Hy, dx, dy, ...
%                     epsxx, epsyy, epszz, boundary);
% [Hz,Ex,Ey,Ez] = postprocess(lambda, neff, Hx, Hy, dx, dy, ...
%                     epsxx, epsxy, epsyx, epsyy, epszz, boundary);
% 
% INPUT:
% 
% lambda - optical wavelength at which mode was calculated
% neff - the calculated effective index of the optial mode
% Hx, Hy - the calculated transverse magnetic fields of the mode
% dx - horizontal grid spacing (vector or scalar)
% dy - vertical grid spacing (vector or scalar)
% eps - index mesh (isotropic materials)  OR:
% epsxx, epsxy, epsyx, epsyy, epszz - index mesh (anisotropic)
% boundary - 4 letter string specifying boundary conditions to be
% applied at the edges of the computation window.  
%   boundary(1) = North boundary condition
%   boundary(2) = South boundary condition
%   boundary(3) = East boundary condition
%   boundary(4) = West boundary condition
% The following boundary conditions are supported: 
%   'A' - Hx is antisymmetric, Hy is symmetric.
%   'S' - Hx is symmetric and, Hy is antisymmetric.
%   '0' - Hx and Hy are zero immediately outside of the
%         boundary. 
% 
% OUTPUT:
% 
% Hz - calculated longitudinal magnetic field.  This output will 
%   have the same dimensions as Hx and Hy.
% Ex, Ey, Ez - calculated electric field.  These field components 
%   are computed at the center of each element instead of on the
%   edges or vertices.
%
% NOTES:
%
% 1) This routine is meant to be used in conjunction with
% wgmodes.m, the vector eigenmode solver.  Please consult the
% help file for wgmodes.m for more information.
%
% 2) The boundary conditions and waveguide specifications
% (given in dx, dy, eps, and boundary) should be the same as
% what was used in wgmodes.m to compute the mode.
%
% 3) The magnetic field components (Hx, Hy, and Hz) are
% calculated at the edges of each cell, whereas the electric
% field components are computed at the center of each cell.
% Therefore if size(eps) = [n,m], then the magnetic fields
% will have a size of [n+1,m+1] while the computed electric
% fields will have a size of [n,m].
%
% 4) Even though wgmodes.m will optionally calculate more than
% one mode at a time, this postprocessing routine must be
% invoked separately for each computed mode.
%
% AUTHORS:  Thomas E. Murphy ([e-mail])

if (nargin == 12)
  epsxx = varargin{1};
  epsxy = varargin{2};
  epsyx = varargin{3};
  epsyy = varargin{4};
  epszz = varargin{5};
  boundary = varargin{6};
elseif (nargin == 10)
  epsxx = varargin{1};
  epsxy = zeros(size(epsxx));
  epsyx = zeros(size(epsxx));
  epsyy = varargin{2};
  epszz = varargin{3};
  boundary = varargin{4};
elseif (nargin == 8)
  epsxx = varargin{1};
  epsxy = zeros(size(epsxx));
  epsyx = zeros(size(epsxx));
  epsyy = epsxx;
  epszz = epsxx;
  boundary = varargin{2};
else
  error('Incorrect number of input arguments.\n');
end

[nx,ny] = size(epsxx);
nx = nx + 1;
ny = ny + 1;

% now we pad eps on all sides by one grid point
epsxx = [epsxx(:,1),epsxx,epsxx(:,ny-1)];
epsxx = [epsxx(1,1:ny+1);epsxx;epsxx(nx-1,1:ny+1)];

epsyy = [epsyy(:,1),epsyy,epsyy(:,ny-1)];
epsyy = [epsyy(1,1:ny+1);epsyy;epsyy(nx-1,1:ny+1)];

epsxy = [epsxy(:,1),epsxy,epsxy(:,ny-1)];
epsxy = [epsxy(1,1:ny+1);epsxy;epsxy(nx-1,1:ny+1)];

epsyx = [epsyx(:,1),epsyx,epsyx(:,ny-1)];
epsyx = [epsyx(1,1:ny+1);epsyx;epsyx(nx-1,1:ny+1)];

epszz = [epszz(:,1),epszz,epszz(:,ny-1)];
epszz = [epszz(1,1:ny+1);epszz;epszz(nx-1,1:ny+1)];

k = 2*pi/lambda;  % free-space wavevector
b = neff*k;       % propagation constant (eigenvalue)

if isscalar(dx)
  dx = dx*ones(nx+1,1);             % uniform grid
else
  dx = dx(:);                       % convert to column vector
  dx = [dx(1);dx;dx(length(dx))];   % pad dx on top and bottom
end

if isscalar(dy)
  dy = dy*ones(1,ny+1);             % uniform grid
else
  dy = dy(:);                       % convert to column vector
  dy = [dy(1);dy;dy(length(dy))]';  % pad dy on top and bottom
end

% distance to neighboring points to north south east and west,
% relative to point under consideration (P), as shown below.

n = ones(1,nx*ny);      n(:) = ones(nx,1)*dy(2:ny+1);
s = ones(1,nx*ny);      s(:) = ones(nx,1)*dy(1:ny);
e = ones(1,nx*ny);      e(:) = dx(2:nx+1)*ones(1,ny);
w = ones(1,nx*ny);      w(:) = dx(1:nx)*ones(1,ny);

% epsilon tensor elements in regions 1,2,3,4, relative to the
% mesh point under consideration (P), as shown below.
%
%                 NW------N------NE
%                 |       |       |
%                 |   1   n   4   |
%                 |       |       |
%                 W---w---P---e---E
%                 |       |       |
%                 |   2   s   3   |
%                 |       |       |
%                 SW------S------SE

exx1 = ones(1,nx*ny);   exx1(:) = epsxx(1:nx,2:ny+1);
exx2 = ones(1,nx*ny);   exx2(:) = epsxx(1:nx,1:ny);
exx3 = ones(1,nx*ny);   exx3(:) = epsxx(2:nx+1,1:ny);
exx4 = ones(1,nx*ny);   exx4(:) = epsxx(2:nx+1,2:ny+1);

eyy1 = ones(1,nx*ny);   eyy1(:) = epsyy(1:nx,2:ny+1);
eyy2 = ones(1,nx*ny);   eyy2(:) = epsyy(1:nx,1:ny);
eyy3 = ones(1,nx*ny);   eyy3(:) = epsyy(2:nx+1,1:ny);
eyy4 = ones(1,nx*ny);   eyy4(:) = epsyy(2:nx+1,2:ny+1);

exy1 = ones(1,nx*ny);   exy1(:) = epsxy(1:nx,2:ny+1);
exy2 = ones(1,nx*ny);   exy2(:) = epsxy(1:nx,1:ny);
exy3 = ones(1,nx*ny);   exy3(:) = epsxy(2:nx+1,1:ny);
exy4 = ones(1,nx*ny);   exy4(:) = epsxy(2:nx+1,2:ny+1);

eyx1 = ones(1,nx*ny);   eyx1(:) = epsyx(1:nx,2:ny+1);
eyx2 = ones(1,nx*ny);   eyx2(:) = epsyx(1:nx,1:ny);
eyx3 = ones(1,nx*ny);   eyx3(:) = epsyx(2:nx+1,1:ny);
eyx4 = ones(1,nx*ny);   eyx4(:) = epsyx(2:nx+1,2:ny+1);

ezz1 = ones(1,nx*ny);   ezz1(:) = epszz(1:nx,2:ny+1);
ezz2 = ones(1,nx*ny);   ezz2(:) = epszz(1:nx,1:ny);
ezz3 = ones(1,nx*ny);   ezz3(:) = epszz(2:nx+1,1:ny);