solve_slab_belta.m

matlab例程

function belta_eff=solve_slab_belta(h,perm)

%global k0;
k0=2*pi/1.55;
belta=linspace(min(perm)*k0,max(perm)*k0,101);
layer_num=length(h);
gamma=zeros(length(belta),layer_num);
for t1=1:layer_num
   for t2=1:length(belta)
      gamma(t2,t1)=sqrt(belta(t2)^2-perm(t1)^2*k0^2);
   end
end

for t1=1:length(belta)
   M=[1 0;0 1];
   for t2=1:layer_num-1
      if gamma(t1,t2+1)==0
         gamma(t1,t2+1)=1;
      end
      D=[gamma(t1,t2+1)+gamma(t1,t2),gamma(t1,t2+1)-gamma(t1,t2);
      	gamma(t1,t2+1)-gamma(t1,t2),gamma(t1,t2+1)+gamma(t1,t2)]/(2*(gamma(t1,t2+1)));
		E=[exp(gamma(t1,t2+1)*h(t2+1)),0;0,exp(-gamma(t1,t2+1)*h(t2+1))];
     	M=E*D*M;
   end
   M11(t1)=M(1,1);
   M12(t1)=M(1,2);
   if abs(M12(t1))==0
      M12(t1)=1;
   end
   rho(t1)=abs(M11(t1))/abs(M12(t1));
end
%figure
%plot(belta,log(rho))
%[temp,t1]=min(rho);
t1=length(belta)-1;
while ~(rho(t1+1)>rho(t1)&rho(t1-1)>rho(t1))
   t1=t1-1;
end
belta_eff=belta(t1);

clear belta gamma rho M12 M12
belta=linspace(belta_eff-k0*(max(perm)-min(perm))/100,...
			belta_eff+k0*(max(perm)-min(perm))/100,101);
gamma=zeros(length(belta),layer_num);
for t1=1:layer_num
   for t2=1:length(belta)
      gamma(t2,t1)=sqrt(belta(t2)^2-perm(t1)^2*k0^2);
   end
end

for t1=1:length(belta)
   M=[1 0;0 1];
   for t2=1:layer_num-1
      if gamma(t1,t2+1)==0
         gamma(t1,t2+1)=1;
      end
      D=[gamma(t1,t2+1)+gamma(t1,t2),gamma(t1,t2+1)-gamma(t1,t2);
      	gamma(t1,t2+1)-gamma(t1,t2),gamma(t1,t2+1)+gamma(t1,t2)]/(2*abs(gamma(t1,t2+1)));
		E=[exp(gamma(t1,t2+1)*h(t2+1)),0;0,exp(-gamma(t1,t2+1)*h(t2+1))];
     	M=E*D*M;
   end
   M11(t1)=M(1,1);
   M12(t1)=M(1,2);
   if abs(M12(t1))==0
      rho(t1)=1e11;
   else
      rho(t1)=abs(M11(t1))/abs(M12(t1));
   end
end
%figure
%plot(belta,log(rho))
%[temp,t1]=min(rho);
t1=length(belta)-1;
while ~(rho(t1+1)>=rho(t1)&rho(t1-1)>=rho(t1))
   t1=t1-1;
end
belta_eff=belta(t1);