fdmtm.m

finite difference method有限差分法求波导的TM模式

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

% define constants
nsize=size(nref) ;
kmag=2*pi/lamda ;
permitivity=nref.^2 ;
jone=complex(0,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(1)*kmag)) ;
end

max=zeros(nsize(1)) ;
max(1,1)=-2/delta^2*permitivity(1)*(1/(permitivity(1)*s(1)^2+permitivity(2)*s(1)*s(2))+1/(2*permitivity(1)*s(1)^2)) ;
max(1,2)=2/delta^2*permitivity(1)*(1/(permitivity(1)*s(1)^2+permitivity(2)*s(1)*s(2))) ;

for i=2:nsize(1)-1
    max(i,i-1)=2/delta^2*permitivity(i)*(1/(permitivity(i)*s(i)^2+permitivity(i-1)*s(i)*s(i-1))) ;
    max(i,i)=-2/delta^2*permitivity(i)*(1/(permitivity(i)*s(i)^2+permitivity(i+1)*s(i)*s(i+1)))-2/delta^2*permitivity(i)*(1/(permitivity(i)*s(i)^2+permitivity(i-1)*s(i)*s(i-1))) ;
    max(i,i+1)=2/delta^2*permitivity(i)*(1/(permitivity(i)*s(i)^2+permitivity(i+1)*s(i)*s(i+1))) ;
end

max(nsize(1),nsize(1)-1)=2/delta^2*permitivity(nsize(1))*(1/(permitivity(nsize(1))*s(nsize(1))^2+permitivity(nsize(1)-1)*s(nsize(1))*s(nsize(1)-1))) ; 
max(nsize(1),nsize(1))=-2/delta^2*permitivity(nsize(1))*(1/(2*permitivity(nsize(1))*s(nsize(1))^2)+1/(permitivity(nsize(1))*s(nsize(1))^2+permitivity(nsize(1)-1)*s(nsize(1))*s(nsize(1)-1))) ;

for i=1:nsize(1)
   k(i,i)=kmag^2*permitivity(i) ;
end


u=k+max ;

%define maximum bata 

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

[V,d]=eigs(u,10,lamsqmax,opts) ;

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

%delta=0.05
%lamda=1.5
%npml=20
%nref=[ones(80,1);1.7*ones(40,1);1.5*ones(80,1)] ;