new.zip

clear % read data %cd('F:\study\meta\matlab') s11=dlmread('s11.txt','',0,0); s21=dlmread('s21.txt','',0,0); s22=dlmread('s22.txt','',0,0); arg_s11=dlmread('arg11.txt','',0,0); arg_s21=dlmread('arg21.txt','',0,0); arg_s22=dlmread('arg22.txt','',0,0); %============================ % input d d =0.0000003; %if isempty(d)==1 % disp('d = default(??????)') % d=10e-3; %end %============================== k=2*pi*s11(:,1)*1e12/3e8; %?????? arg1= arg_s11(:,2)./360*(2*pi); arg2= arg_s21(:,2)./360*(2*pi); arg3= arg_s22(:,2)./360*(2*pi); S11=s11(:,2).*exp(-1i*arg1); S21=s21(:,2).*exp(-1i*arg2); S22=s22(:,2).*exp(-1i*arg3); S11=sqrt(S11.*S22); [row,column]=size(S11); z=sqrt(((1+S11).*(1+S11)-S21.*S21)./((1-S11).*(1-S11)-S21.*S21)); % determine z delta_z=0.1; for ii=1:row if abs(real(z(ii)))>=delta_z if real(z(ii))<0 z(ii)=-z(ii); end else if abs(S21(ii)/(1-S11(ii)*(z(ii)-1)/(z(ii)+1)))>1 z(ii)=-z(ii); end end end % determine n imag_n=-real(log(S21./(1-S11.*(z-1)./(z+1))))./(k*d); real_n=imag(log(S21./(1-S11.*(z-1)./(z+1))))./(k*d); n=real_n+1i*imag_n; fig=1; count=0; n_temp=n(1); for nn=-100:200 n(1)=n_temp+2*nn*pi/(k(1)*d); if abs(real(n(1))*imag(z(1)))<=imag(n(1))*real(z(1))+100 % determine the first real(n) %???????????????eps??mu妊??0??????????妊????? % determine the branch of real(n) nnn=nn;count=count+1; for ii=1:row-1 % two roots delta1=-1+sqrt(1-2*(1-(S21(ii+1)/(1-S11(ii+1)*(z(ii+1)-1)/(z(ii+1)+1)))/(S21(ii)/(1-S11(ii)*(z(ii)-1)/(z(ii)+1))))); delta2=-1-sqrt(1-2*(1-(S21(ii+1)/(1-S11(ii+1)*(z(ii+1)-1)/(z(ii+1)+1)))/(S21(ii)/(1-S11(ii)*(z(ii)-1)/(z(ii)+1))))); n1=(delta1+1i*n(ii)*k(ii)*d)/(1i*k(ii+1)*d); n2=(delta2+1i*n(ii)*k(ii)*d)/(1i*k(ii+1)*d); % determine which root is the correct one, denoted as n0 if abs(imag(n1)-imag(n(ii+1)))<abs(imag(n2)-imag(n(ii+1))) n0=n1; else n0=n2; end % determine m by comparing with n0 m=0; delta_n=real_n(ii+1)-real(n0); if delta_n<0 mm=1; while abs(real(n0)-(real_n(ii+1)+mm*2*pi/(k(ii+1)*d)))<abs(delta_n) delta_n=(real_n(ii+1)+mm*2*pi/(k(ii+1)*d))-real(n0); m=mm; mm=mm+1; end else mm=-1; while abs(real(n0)-(real_n(ii+1)+mm*2*pi/(k(ii+1)*d)))<abs(delta_n) delta_n=(real_n(ii+1)+mm*2*pi/(k(ii+1)*d))-real(n0); m=mm; mm=mm-1; end end n(ii+1)=real_n(ii+1)+m*2*pi/(k(ii+1)*d)+1i*imag_n(ii+1); end % eps & mue eps=n./z; mu=n.*z; % plot figure(fig) subplot(2,2,1); plot(s11(:,1),real(n),'k.',s11(:,1),imag(n),'y.') %axis([9.3 9.4 -5 5]) xlabel('Frequency/THz') title('n') subplot(2,2,2); plot(s11(:,1),real(z),'k.',s11(:,1),imag(z),'y.') title('z') subplot(2,2,3); plot(s11(:,1),real(eps),'k',s11(:,1),imag(eps),'r--') legend('real(Eps)','imag(Eps)') %axis([2 4 -10 100]) title('Eps') subplot(2,2,4); plot(s11(:,1),real(mu),'k',s11(:,1),imag(mu),'r--') legend('Real(Mu)','imag(Mu)') %axis([9.2 9.6 -10 10]) title('Mu') fig=fig+1; end end