femte.asv

有限元法(finite element method)求波导中的TE模

function[V,d,neff] =femte(nref,delta,lamda,npml) 

% define constants
nsize=size(nref) ;
kmag=2*pi/lamda ;
permitivity=nref.^2 ;
jone=complex(0,1) ;
le=delta(1); 

% define PML
s=ones(nsize(1),1) ;
for i=1:npml 
    sigmax=((npml+1-i)/npml).^4 ;
    s(i)=(1+jone*sigmax.*376.7343/(permitivity(1)*kmag)) ;
    s(nsize(1)-i+1)=(1+jone*sigmax.*376.7343/(permitivity(nsize(1))*kmag)) ;
end

% define element matrix
m1=[1,-1;-1,1] ;
m2=[2,1;1,2] ;

p=zeros(nsize(1),nsize(1)) ;
for i=1:nsize(1)-1
    p(i:i+1,i:i+1)=p(i:i+1,i:i+1)+m1 ;
end
r=zeros(nsize(1),nsize(1)) ;
for i=1:nsize(1)-1
    r(i:i+1,i:i+1)=r(i:i+1,i:i+1)+m2 ;
end
q=zeros(nsize(1),nsize(1)) ;
for i=1:nsize(1)-1
    q(i:i+1,i:i+1)=q(i:i+1,i:i+1)+m2*nref(i)^2; 
end

% define global matrix
s2=1./s ;
pglobal=zeros(nsize(1),nsize(1));
for i=1:nsize(1)
  pglobal(i,:)=p(i,:).*s2(i)*(-1/le) ;
end

rglobal=le/6*r ;

qglobal=kmag^2*le/6*q ;

% define maximum bata 
lamsqmax=(1.7*kmag)^2;
opts.tol=1.e-19 ;

%[v,d]=eigs((pglobal+qglobal),rglobal,10,lamsqmax,opts) ;
[V,d]=eig((pglobal+qglobal),rglobal) ;

for i=1:260
    neff(i,1)=sqrt(d(i,i))/kmag ;
     if (imag(neff(i,1))>1.e-4)
       neff(i,1)=0.0 ;
     end
end