wgmodes.m

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

function [phix,phiy,neff] = wgmodes (lambda, guess, nmodes, dx, dy, varargin);

% This function computes the two transverse magnetic field
% components of a dielectric waveguide, using the finite
% difference method.  For details about the method, please
% consult:  
%
% A. B. Fallahkhair, K. S. Li, and T. E. Murphy, "Vector
% Finite Difference Modesolver for Anisotropic Dielectric
% Waveguides", submitted to J. Lightwave. Technol., 2007.
%
% USAGE:
% 
% [hx,hy,neff] = wgmodes(lambda, guess, nmodes, dx, dy, ...
%                        eps,boundary);
% [hx,hy,neff] = wgmodes(lambda, guess, nmodes, dx, dy, ...
%                        epsxx, epsyy, epszz, boundary);
% [hx,hy,neff] = wgmodes(lambda, guess, nmodes, dx, dy, ...
%                        epsxx, epsxy, epsyx, epsyy, epszz, boundary);
% 
% INPUT:
% 
% lambda - optical wavelength
% guess - scalar shift to apply when calculating the eigenvalues.
%     This routine will return the eigenpairs which have an
%     effective index closest to this guess
% nmodes - the number of modes to calculate
% 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:
% 
% hx - three-dimensional vector containing Hx for each
%      calculated mode 
% hy - three-dimensional vector containing Hy for each
%      calculated mode (e.g.: hy(:,k) = two dimensional Hy
%      matrix for the k-th mode 
% neff - vector of modal effective indices
%
% NOTES:
%
% 1) The units are arbitrary, but they must be self-consistent
% (e.g., if lambda is in um, then dx and dy should also be in
% um.
%
% 2) Unlike the E-field modesolvers, this method calculates
% the transverse MAGNETIC field components Hx and Hy.  Also,
% it calculates the components at the edges (vertices) of
% each cell, rather than in the center of each cell.  As a
% result, if size(eps) = [n,m], then the output eigenvectors
% will be have a size of [n+1,m+1].
%
% 3) This version of the modesolver can optionally support
% non-uniform grid sizes.  To use this feature, you may let dx
% and/or dy be vectors instead of scalars.
%
% 4) The modesolver can consider anisotropic materials, provided
% the permittivity of all constituent materials can be
% expressed in one of the following forms:   
%
%  [eps  0   0 ]  [epsxx   0     0  ]  [epsxx epsxy   0  ]
%  [ 0  eps  0 ]  [  0   epsyy   0  ]  [epsyx epsyy   0  ]
%  [ 0   0  eps]  [  0     0   epszz]  [  0     0   epszz]
%
% The program will decide which form is appropriate based upon
% the number of input arguments supplied.
%
% 5) Perfectly matched boundary layers can be accomodated by
% using the complex coordinate stretching technique at the
% edges of the computation window.  (stretchmesh.m can be used
% for complex or real-coordinate stretching.)
%
% AUTHORS:  Thomas E. Murphy ([e-mail])
%           Arman B. Fallahkhair ([e-mail])
%           Kai Sum Li ([e-mail])

if (nargin == 11)
  epsxx = varargin{1};
  epsxy = varargin{2};
  epsyx = varargin{3};
  epsyy = varargin{4};
  epszz = varargin{5};
  boundary = varargin{6};
elseif (nargin == 9)
  epsxx = varargin{1};
  epsxy = zeros(size(epsxx));
  epsyx = zeros(size(epsxx));
  epsyy = varargin{2};
  epszz = varargin{3};
  boundary = varargin{4};
elseif (nargin == 7)
  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

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);
ezz4 = ones(1,nx*ny);   ezz4(:) = epszz(2:nx+1,2:ny+1);

ns21 = n.*eyy2+s.*eyy1;
ns34 = n.*eyy3+s.*eyy4;
ew14 = e.*exx1+w.*exx4;
ew23 = e.*exx2+w.*exx3;

axxn = ((2*eyy4.*e-eyx4.*n).*(eyy3./ezz4)./ns34 + ...
        (2*eyy1.*w+eyx1.*n).*(eyy2./ezz1)./ns21)./(n.*(e+w));

axxs = ((2*eyy3.*e+eyx3.*s).*(eyy4./ezz3)./ns34 + ...
        (2*eyy2.*w-eyx2.*s).*(eyy1./ezz2)./ns21)./(s.*(e+w));

ayye = (2.*n.*exx4 - e.*exy4).*exx1./ezz4./e./ew14./(n+s) + ...
       (2.*s.*exx3 + e.*exy3).*exx2./ezz3./e./ew23./(n+s);

ayyw = (2.*exx1.*n + exy1.*w).*exx4./ezz1./w./ew14./(n+s) + ...
       (2.*exx2.*s - exy2.*w).*exx3./ezz2./w./ew23./(n+s);

axxe = 2./(e.*(e+w)) + ...
       (eyy4.*eyx3./ezz3 - eyy3.*eyx4./ezz4)./(e+w)./ns34;

axxw = 2./(w.*(e+w)) + ...
       (eyy2.*eyx1./ezz1 - eyy1.*eyx2./ezz2)./(e+w)./ns21;

ayyn = 2./(n.*(n+s)) + ...
       (exx4.*exy1./ezz1 - exx1.*exy4./ezz4)./(n+s)./ew14;

ayys = 2./(s.*(n+s)) + ...
       (exx2.*exy3./ezz3 - exx3.*exy2./ezz2)./(n+s)./ew23;

axxne = +eyx4.*eyy3./ezz4./(e+w)./ns34;
axxse = -eyx3.*eyy4./ezz3./(e+w)./ns34;
axxnw = -eyx1.*eyy2./ezz1./(e+w)./ns21;
axxsw = +eyx2.*eyy1./ezz2./(e+w)./ns21;

ayyne = +exy4.*exx1./ezz4./(n+s)./ew14;
ayyse = -exy3.*exx2./ezz3./(n+s)./ew23;
ayynw = -exy1.*exx4./ezz1./(n+s)./ew14;
ayysw = +exy2.*exx3./ezz2./(n+s)./ew23;

axxp = - axxn - axxs - axxe - axxw - axxne - axxse - axxnw - axxsw ...
       + k^2*(n+s).*(eyy4.*eyy3.*e./ns34 + eyy1.*eyy2.*w./ns21)./(e+w);

ayyp = - ayyn - ayys - ayye - ayyw - ayyne - ayyse - ayynw - ayysw ...
       + k^2*(e+w).*(exx1.*exx4.*n./ew14 + exx2.*exx3.*s./ew23)./(n+s);

axyn = (eyy3.*eyy4./ezz4./ns34 - ...
        eyy2.*eyy1./ezz1./ns21 + ...
        s.*(eyy2.*eyy4 - eyy1.*eyy3)./ns21./ns34)./(e+w);

axys = (eyy1.*eyy2./ezz2./ns21 - ...
        eyy4.*eyy3./ezz3./ns34 + ...
        n.*(eyy2.*eyy4 - eyy1.*eyy3)./ns21./ns34)./(e+w);

ayxe = (exx1.*exx4./ezz4./ew14 - ...