plot_modes.m

利用matlab编程实现含负折射材料的表面波和波导特性

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculates beta nd then plots the corresponding mode
% formulas used shown during lessons
% the program calls the subroutine disprel_modes in order
% to numerically evaluate the disp relation zero and thus
% find beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,Ey]=plot_modes(e1,u1,e2,u2,e3,u3,h,m)

 global h  e1 u1 e2 u2 e3 u3 lambda k0 kc ks kf phis phic m 
 c = 3e8*1e6*1e-15;  % in microns/fs
%  h = 5;
 
%  u1=1;
%  e1=1.2^2;
% %  n2=sqrt(e2)*sqrt(u2);
%  
%  u3=1;
%  e3=1.2^2;
% %  n3=sqrt(e3)*sqrt(u3);
%  
%  u2=-1;
%  e2=-1.8^2;
%  n1=sqrt(e1)*sqrt(u1);
 
% n3 = 1.45;              % substrate refractive index
% n2 = 1.45;               % cover refractive index
% n1= 1.48;                % film refractive index
% m =2;                  % mode order
 
lambda = 1.8;          % mode wavelength
k0 = 2*pi./lambda;

if e2*u2>e1*u1&e2*u2>e3*u3

        

beta0 = [1];           % Make a starting guess at the solution
[beta,fval] = fsolve(@disprel_modes1,beta0,optimset('fsolve')) ; % Call zero-finding subroutine

beta                   % print the vlaue of beta in the command window

neff = beta./k0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the waveguide mode spatial profile
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=(kf/u2-ks/u3)/2;
C=(kf/u2+kc/u1)/2;
D=B*exp(-i*kf*h)+C*exp(i*kf*h);
step = 1e-2;
ext = 2;
xc = [h+step:step:h+ext]; xf = [0:step:h]; xs =[-ext:step:-step];   
x = [xs xf xc];
Eyc = cos(kf*h-phis)*exp(-kc.*(xc-h)); Eyf = cos(kf.*xf-phis); Eys = cos(-phis).*exp(ks.*xs);
% Eyc =D*exp(-kc.*(xc-h)); Eyf = B*exp(-i*kf.*xf)+C*exp(i*kf.*xf); Eys = (B+C)*exp(ks.*xs);
% Eyc =D*exp(-kc.*(xc-h)); Eyf = B*cos(kf.*xf); Eys = (B+C)*exp(ks.*xs);
Ey = [Eys Eyf Eyc];
if m==1&e2<0&u2<0&Ey<1e-3;
figure, plot(x,Ey)
line([0 0],[-5 5]); line([h h],[-5 5]);
line([-2 h+2],[0 0]);
text(h/2,3.2,['TE_',num2str(m)])
text(-ext+0.5,1.8,['e_s=',num2str(e3)])
text(ext+h-1.5,1.8,['e_c=',num2str(e1)])
text(h/2-0.2,1.8,['e_f=',num2str(e2)])
text(-ext+0.5,1.2,['u_s=',num2str(u3)])
text(ext+h-1.5,1.2,['u_c=',num2str(u1)])
text(h/2-0.2,1.2,['u_f=',num2str(u2)])
text(h/2-1,-1.8,['h=',num2str(h)])
text(h/2+0.0,-1.8,['    lambda=',num2str(lambda)])
text(-1,-2.8,['TE_1的存在条件不满足，请重新设定参数'])
xlabel('x (\mum)')
ylabel('Ey ')    
else
% plot the mode profiles
if m==0&e2<0&u2<0
figure, plot(x,Ey)
line([0 0],[-5 5]); line([h h],[-5 5]);
line([-2 h+2],[0 0]);
text(h/2,3.2,['TE_',num2str(m)])
text(-ext+0.5,1.8,['e_s=',num2str(e3)])
text(ext+h-1.5,1.8,['e_c=',num2str(e1)])
text(h/2-0.2,1.8,['e_f=',num2str(e2)])
text(-ext+0.5,1.2,['u_s=',num2str(u3)])
text(ext+h-1.5,1.2,['u_c=',num2str(u1)])
text(h/2-0.2,1.2,['u_f=',num2str(u2)])
text(h/2-1,-1.8,['h=',num2str(h)])
text(h/2+0.0,-1.8,['    lambda=',num2str(lambda)])
text(0,-2.8,['       左手物质不存在基模'])
xlabel('x (\mum)')
ylabel('Ey ') 
else    
% figure, plot(x,abs(Ey).^2)
figure, plot(x,Ey)
line([0 0],[-5 5]); line([h h],[-5 5]);
line([-2 h+2],[0 0]);
text(h/2,3.2,['TE_',num2str(m)])
text(-ext+0.5,1.8,['e_s=',num2str(e3)])
text(ext+h-1.5,1.8,['e_c=',num2str(e1)])
text(h/2-0.2,1.8,['e_f=',num2str(e2)])
text(-ext+0.5,1.2,['u_s=',num2str(u3)])
text(ext+h-1.5,1.2,['u_c=',num2str(u1)])
text(h/2-0.2,1.2,['u_f=',num2str(u2)])
text(h/2-1,-1.8,['h=',num2str(h)])
text(h/2+0.0,-1.8,['    lambda=',num2str(lambda)])
text(-1,-2.8,['    条件是 e_f*u_f>e_c*u_c and e_f*u_f>e_s*u_s'])
xlabel('x (\mum)')
ylabel('Ey ')
end
end

    


 
else 
    x=0;
    y=0;
plot(x,y)   
 text(-1,0,'条件不能满足，请重新设定参数')

end