wgmodes.m

来自「强大的计算电磁场本征函数与本征模的程序」· M 代码 · 共 499 行 · 第 1/2 页

        exx2.*exx3./ezz3./ew23 + ...
        w.*(exx2.*exx4 - exx1.*exx3)./ew23./ew14)./(n+s);

ayxw = (exx3.*exx2./ezz2./ew23 - ...
        exx4.*exx1./ezz1./ew14 + ...
        e.*(exx4.*exx2 - exx1.*exx3)./ew23./ew14)./(n+s);

axye = (eyy4.*(1-eyy3./ezz3) - eyy3.*(1-eyy4./ezz4))./ns34./(e+w) - ...
       2*(eyx1.*eyy2./ezz1.*n.*w./ns21 + ...
          eyx2.*eyy1./ezz2.*s.*w./ns21 + ...
          eyx4.*eyy3./ezz4.*n.*e./ns34 + ...
          eyx3.*eyy4./ezz3.*s.*e./ns34 + ...
          eyy1.*eyy2.*(1./ezz1-1./ezz2).*w.^2./ns21 + ...
          eyy3.*eyy4.*(1./ezz4-1./ezz3).*e.*w./ns34)./e./(e+w).^2;

axyw = (eyy2.*(1-eyy1./ezz1) - eyy1.*(1-eyy2./ezz2))./ns21./(e+w) - ...
       2*(eyx4.*eyy3./ezz4.*n.*e./ns34 + ...
          eyx3.*eyy4./ezz3.*s.*e./ns34 + ...
          eyx1.*eyy2./ezz1.*n.*w./ns21 + ...
          eyx2.*eyy1./ezz2.*s.*w./ns21 + ...
          eyy4.*eyy3.*(1./ezz3-1./ezz4).*e.^2./ns34 + ...
          eyy2.*eyy1.*(1./ezz2-1./ezz1).*w.*e./ns21)./w./(e+w).^2;

ayxn = (exx4.*(1-exx1./ezz1) - exx1.*(1-exx4./ezz4))./ew14./(n+s) - ...
       2*(exy3.*exx2./ezz3.*e.*s./ew23 + ...
          exy2.*exx3./ezz2.*w.*s./ew23 + ...
          exy4.*exx1./ezz4.*e.*n./ew14 + ...
          exy1.*exx4./ezz1.*w.*n./ew14 + ...
          exx3.*exx2.*(1./ezz3-1./ezz2).*s.^2./ew23 + ...
          exx1.*exx4.*(1./ezz4-1./ezz1).*n.*s./ew14)./n./(n+s).^2;

ayxs = (exx2.*(1-exx3./ezz3) - exx3.*(1-exx2./ezz2))./ew23./(n+s) - ...
       2*(exy4.*exx1./ezz4.*e.*n./ew14 + ...
          exy1.*exx4./ezz1.*w.*n./ew14 + ...
          exy3.*exx2./ezz3.*e.*s./ew23 + ...
          exy2.*exx3./ezz2.*w.*s./ew23 + ...
          exx4.*exx1.*(1./ezz1-1./ezz4).*n.^2./ew14 + ...
          exx2.*exx3.*(1./ezz2-1./ezz3).*s.*n./ew23)./s./(n+s).^2;

axyne = +eyy3.*(1-eyy4./ezz4)./(e+w)./ns34;
axyse = -eyy4.*(1-eyy3./ezz3)./(e+w)./ns34;
axynw = -eyy2.*(1-eyy1./ezz1)./(e+w)./ns21;
axysw = +eyy1.*(1-eyy2./ezz2)./(e+w)./ns21;

ayxne = +exx1.*(1-exx4./ezz4)./(n+s)./ew14;
ayxse = -exx2.*(1-exx3./ezz3)./(n+s)./ew23;
ayxnw = -exx4.*(1-exx1./ezz1)./(n+s)./ew14;
ayxsw = +exx3.*(1-exx2./ezz2)./(n+s)./ew23;

axyp = -(axyn + axys + axye + axyw + axyne + axyse + axynw + axysw) ...
       - k^2.*(w.*(n.*eyx1.*eyy2 + s.*eyx2.*eyy1)./ns21 + ...
               e.*(s.*eyx3.*eyy4 + n.*eyx4.*eyy3)./ns34)./(e+w);

ayxp = -(ayxn + ayxs + ayxe + ayxw + ayxne + ayxse + ayxnw + ayxsw) ...
       - k^2.*(n.*(w.*exy1.*exx4 + e.*exy4.*exx1)./ew14 + ...
               s.*(w.*exy2.*exx3 + e.*exy3.*exx2)./ew23)./(n+s);  

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

axxs(ib)  = axxs(ib)  + sign*axxn(ib);
axxse(ib) = axxse(ib) + sign*axxne(ib);
axxsw(ib) = axxsw(ib) + sign*axxnw(ib);
ayxs(ib)  = ayxs(ib)  + sign*ayxn(ib);
ayxse(ib) = ayxse(ib) + sign*ayxne(ib);
ayxsw(ib) = ayxsw(ib) + sign*ayxnw(ib);
ayys(ib)  = ayys(ib)  - sign*ayyn(ib);
ayyse(ib) = ayyse(ib) - sign*ayyne(ib);
ayysw(ib) = ayysw(ib) - sign*ayynw(ib);
axys(ib)  = axys(ib)  - sign*axyn(ib);
axyse(ib) = axyse(ib) - sign*axyne(ib);
axysw(ib) = axysw(ib) - sign*axynw(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

axxn(ib)  = axxn(ib)  + sign*axxs(ib);
axxne(ib) = axxne(ib) + sign*axxse(ib);
axxnw(ib) = axxnw(ib) + sign*axxsw(ib);
ayxn(ib)  = ayxn(ib)  + sign*ayxs(ib);
ayxne(ib) = ayxne(ib) + sign*ayxse(ib);
ayxnw(ib) = ayxnw(ib) + sign*ayxsw(ib);
ayyn(ib)  = ayyn(ib)  - sign*ayys(ib);
ayyne(ib) = ayyne(ib) - sign*ayyse(ib);
ayynw(ib) = ayynw(ib) - sign*ayysw(ib);
axyn(ib)  = axyn(ib)  - sign*axys(ib);
axyne(ib) = axyne(ib) - sign*axyse(ib);
axynw(ib) = axynw(ib) - sign*axysw(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

axxw(ib)  = axxw(ib)  + sign*axxe(ib);
axxnw(ib) = axxnw(ib) + sign*axxne(ib);
axxsw(ib) = axxsw(ib) + sign*axxse(ib);
ayxw(ib)  = ayxw(ib)  + sign*ayxe(ib);
ayxnw(ib) = ayxnw(ib) + sign*ayxne(ib);
ayxsw(ib) = ayxsw(ib) + sign*ayxse(ib);
ayyw(ib)  = ayyw(ib)  - sign*ayye(ib);
ayynw(ib) = ayynw(ib) - sign*ayyne(ib);
ayysw(ib) = ayysw(ib) - sign*ayyse(ib);
axyw(ib)  = axyw(ib)  - sign*axye(ib);
axynw(ib) = axynw(ib) - sign*axyne(ib);
axysw(ib) = axysw(ib) - sign*axyse(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

axxe(ib)  = axxe(ib)  + sign*axxw(ib);
axxne(ib) = axxne(ib) + sign*axxnw(ib);
axxse(ib) = axxse(ib) + sign*axxsw(ib);
ayxe(ib)  = ayxe(ib)  + sign*ayxw(ib);
ayxne(ib) = ayxne(ib) + sign*ayxnw(ib);
ayxse(ib) = ayxse(ib) + sign*ayxsw(ib);
ayye(ib)  = ayye(ib)  - sign*ayyw(ib);
ayyne(ib) = ayyne(ib) - sign*ayynw(ib);
ayyse(ib) = ayyse(ib) - sign*ayysw(ib);
axye(ib)  = axye(ib)  - sign*axyw(ib);
axyne(ib) = axyne(ib) - sign*axynw(ib);
axyse(ib) = axyse(ib) - sign*axysw(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);

Axx = sparse ([iall,iw,ie,is,in,ine,ise,isw,inw], ...
	[iall,ie,iw,in,is,isw,inw,ine,ise], ...
	[axxp(iall),axxe(iw),axxw(ie),axxn(is),axxs(in), ...
     axxsw(ine),axxnw(ise),axxne(isw),axxse(inw)]);

Axy = sparse ([iall,iw,ie,is,in,ine,ise,isw,inw], ...
	[iall,ie,iw,in,is,isw,inw,ine,ise], ...
	[axyp(iall),axye(iw),axyw(ie),axyn(is),axys(in), ...
     axysw(ine),axynw(ise),axyne(isw),axyse(inw)]);

Ayx = sparse ([iall,iw,ie,is,in,ine,ise,isw,inw], ...
	[iall,ie,iw,in,is,isw,inw,ine,ise], ...
	[ayxp(iall),ayxe(iw),ayxw(ie),ayxn(is),ayxs(in), ...
     ayxsw(ine),ayxnw(ise),ayxne(isw),ayxse(inw)]);

Ayy = sparse ([iall,iw,ie,is,in,ine,ise,isw,inw], ...
	[iall,ie,iw,in,is,isw,inw,ine,ise], ...
	[ayyp(iall),ayye(iw),ayyw(ie),ayyn(is),ayys(in), ...
     ayysw(ine),ayynw(ise),ayyne(isw),ayyse(inw)]);

A = [[Axx Axy];[Ayx Ayy]];

fprintf(1,'nnz(A) = %d\n',nnz(A));

shift = (guess*k)^2;
options.tol = 1e-8;
options.disp = 0;						% suppress output

clear Axx Axy Ayx Ayy ...
    axxnw axxne axxne ...
    axxw  axxp  axxe ...
    axxsw axxse axxse ...
    axynw axyne axyne ...
    axyw  axyp  axye ...
    axysw axyse axyse ...
    ayynw ayyne ayyne ...
    ayyw  ayyp  ayye ...
    ayysw ayyse ayyse ...
    ayxnw ayxne ayxne ...
    ayxw  ayxp  ayxe ...
    ayxsw ayxse ayxse ...
    iall ie iw in iw ...
    isw inw ine ise ...
    exx1 exx2 exx3 exx4 ...
    exy1 exy2 exy3 exy4 ...
    eyx1 eyx2 eyx3 eyx4 ...
    eyy1 eyy2 eyy3 eyy4 ...
    ezz1 ezz2 ezz3 ezz4 ...
    ns21 ns34 ew14 ew23;

[v,d] = eigs(A,speye(size(A)),nmodes,shift,options);
neff = lambda*sqrt(diag(d))/(2*pi);

phix = zeros(nx,ny,nmodes);
phiy = zeros(nx,ny,nmodes);
temp = zeros(nx,2*ny);

% Normalize modes

temp = zeros(nx*ny,2);
for kk = 1:nmodes;
  temp(:) = v(:,kk);
  [mag,ii] = max(sqrt(sum(abs(temp).^2,2)));
  if abs(temp(ii,1)) > abs(temp(ii,2)),
    jj = 1;
  else 
    jj = 2;
  end
  mag = mag*temp(ii,jj)/abs(temp(ii,jj));
  temp = temp/mag;
  phix(:,:,kk) = reshape(temp(:,1),nx,ny);
  phiy(:,:,kk) = reshape(temp(:,2),nx,ny);
end;

return;