s_slabinair.m

该matlab函数计算自由空间里有限厚度、无限大的各向同性媒质块（允许频率色散）在任意极化、任意入射角度的来波下的反射和透射系数

% Usage:
%   in exp(-i*w*t) system. i is reserved for representing imaginary sign
%   z is longitudinal direction.
%
% Format: [S11, S21]=S_SlabinAir(f,epr1,mur1,d1,deembed_inc,deembed_trans,theta_inc,cho, Fig, Keep)
%
% Input arguments:
%   f: frequency band concerned, in GHz;
%   epr1: relative permittivity of the slab;
%   mur1: relative permeability of the slab;
%       f,epr1,mur1: all column vectors (cf*1).
%   d1: slab thichness, in mm;
%   deembed_inc,deembed_trans: extruded background (air) length at the incident-wave/tansmitted-wave side, in mm (indicating reference plane); 
%   theta_inc: incident angle with respect to normal in deg (0 ~ 90); 
%   cho: string type, indicating polarization type of incident wave ('s': TEtoz or 'p': TMtoz).
%   Fig: boolean type, plot figure if its value is true.
%   Keep: Boolean type, indicating whether you want to keep all figures after the execution of this function when Fig=true. Set it as false when Fig=false.
%
% Output arguments:
%   S11: cf*1 colume vector, with respect to the incident interface between air and slab (when deembed_inc=0).
%   S21: cf*1 colume vector, with respect to the incident & transmit interface between air and slab (when deembed_inc=0,deembed_trans=0).
%       Note: If deembed_inc~=0, the reference plane at the incident side of the slabs move away from slabs by the amount indicated by deembed_inc.
%       Note: If deembed_trans~=0, the reference plane at the transmit side of the slabs move away from slabs by the amount indicated by deembed_trans.
function [S11, S21]=S_SlabinAir(f,epr1,mur1,d1,deembed_inc,deembed_trans,theta_inc,cho, Fig, Keep)
cf=length(f);
if length(epr1)==1
    epr1=epr1*ones(cf,1);
end
if length(mur1)==1
    mur1=mur1*ones(cf,1);
end

ep0=1/(36*pi*10^9);
mu0=4*pi*10^(-7);
w=2*pi*f*10^9; %rad/s
theta_inc=pi*theta_inc/180;

ep1=epr1*ep0;
mu1=mur1*mu0;

k0=w*sqrt(ep0*mu0); kx=k0*sin(theta_inc); k0z=k0*cos(theta_inc); 
k1=w.*sqrt(ep1.*mu1); k1z=sqrt(k1.^2-kx.^2);
for j=1:length(epr1)
    if real(sqrt(epr1(j)*mur1(j))/epr1(j))<0
        k1z(j)=-k1z(j);
    end
end

if cho=='s' %TEtoz
    p10=(mu1.*k0z)./(mu0.*k1z);
    p21=(mu0.*k1z)./(mu1.*k0z);
elseif cho=='p' %TMtoz
    p10=(ep1.*k0z)./(ep0.*k1z);
    p21=(ep0.*k1z)./(ep1.*k0z);
else
    error('incorrect value of "cho". "cho" should be "s" or "p" only!');
end
R10=(1-p10)./(1+p10);
R21=(1-p21)./(1+p21);

%R & T are with respect to the incident interface between air and slab.
%S11 is with respect to the incident interface between air and slab.
%S21 is with respect to the incident & transmit interface between air and slab.
for j=1:length(w)
    V10=(1+p10(j))*[exp(-i*k1z(j)*d1*10^-3), R10(j)*exp(-i*k1z(j)*d1*10^-3); R10(j)*exp(i*k1z(j)*d1*10^-3), exp(i*k1z(j)*d1*10^-3)]/2;
    V21=(1+p21(j))*[exp(i*k0z(j)*d1*10^-3), R21(j)*exp(i*k0z(j)*d1*10^-3); R21(j)*exp(-i*k0z(j)*d1*10^-3), exp(-i*k0z(j)*d1*10^-3)]/2;
    V20=V21*V10;
    R(j)=-V20(1,2)/V20(1,1);
    T(j)=V20(2,1)*R(j)+V20(2,2);
    S11(j)=R(j);
    S21(j)=T(j)*exp(i*k0z(j)*d1*10^-3);
end
S11=S11.';S21=S21.';

%Modify S11 and S21 to adjust them to the reference plane
for j=1:length(w)
    S11(j)=S11(j)*exp(2*i*k0z(j)*deembed_inc*10^-3);
    S21(j)=S21(j)*exp(i*k0z(j)*(deembed_inc+deembed_trans)*10^-3);
end

if (Fig)
    figure(201)
    plot(f,20*log10(abs(S11)),'-b',f,20*log10(abs(S21)),'-r','lineWidth',4)
    xlabel('f/GHz','fontsize',30)
    ylabel('dB','fontsize',30)
    set(gca,'linewidth',2,'fontsize',30);
    legend('dB(S11)','dB(S21)');
    grid on
    figure(202)
    plot(f,180*angle(S11)/pi,'-b',f,180*angle(S21)/pi,'-r','lineWidth',4)
    set(gca,'linewidth',2,'fontsize',30);
    xlabel('f/GHz','fontsize',30)
    ylabel('deg','fontsize',30)
    legend('deg(S11)','deg(S21)')
    grid on
    display('press "enter" to continue...');
    pause
    if ~Keep
        close all
    end
end