www.pudn.com > Slab.rar > solve_slab_belta.m, change:2002-03-20,size:1970b


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);