clc;clear all;
eps=1/36/pi*1e-9;
mu=4*pi*1e-7;
eps0=eps;
mu0=mu;
c=1/sqrt(eps*mu);
w=3*10^12;
lamda=c*2*pi/w;
d1=lamda/16;
d2=lamda/16;
eps1=-3;mu1=3;
eps2=3;mu2=-3;
for c5=0:899;
theta=(c5+0.1)*pi/1800;
kz1=(w/c)*sqrt(mu1*eps1)*sqrt(1-sin(theta)*sin(theta)/eps1/mu1);
q0=cos(theta);%第0层,真空
qs=q0;%基底
q1=sqrt(eps1*mu1-sin(theta)*sin(theta))/mu1;%第一层为负电 TE模式
o11=cos(kz1*d1);o12=i/q1*sin(kz1*d1);
o21=i*q1*sin(kz1*d1);o22=cos(kz1*d1);
M=[o11 o12;o21 o22];
q2=sqrt(eps2*mu2-sin(theta)*sin(theta))/mu2;%第2层为负磁
kz2=(w/c)*sqrt(mu2*eps2)*sqrt(1-sin(theta)*sin(theta)/eps2/mu2);
p11=cos(kz2*d2);p12=i*sin(kz2*d2)/q2;
p21=i*q2*sin(kz2*d2);p22=cos(kz2*d2);
N=[p11 p12;p21 p22];
ds=0.002*lamda;%夹层真空
kzs=w/c*cos(theta);
q0=cos(theta);%第0层,右手材料
qs=q0;%基底
s11=cos(kzs*ds);s12=i*sin(kzs*ds)/qs;
s21=i*qs*sin(kzs*ds);s22=cos(kzs*ds);
S=[s11 s12;s21 s22];%夹层真空
X=M*S*N;
x11=X(1,1);
x12=X(1,2);
x21=X(2,1);
x22=X(2,2);
r=((q0*x22-q0*x11)-q0*q0*x12+x21)/(q0*x22+q0*x11-q0*q0*x12-x21);
t=2*q0/((q0*x22+q0*x11)-(q0*q0*x12+x21));
phi=atan(imag(r)/real(r)) ;
s=c5+1;
xx11(1,s)=x11;xx12(1,s)=x12;xx21(1,s)=x21;xx22(1,s)=x22;
g1(1,s)=r;
k1(1,s)=phi;
rabs=abs(r);
R(1,s)=abs(r)*abs(r);
h1(1,s)=t;
tabs(1,s)=abs(t);
T(1,s)=abs(t)*abs(t);
ri(1,s)=imag(r);
rr(1,s)=real(r);
phit=atan(imag(t)/real(t));
t1(1,s)=phit;
end
for h=1:899;
k3(h)=-k1(h+1)+k1(h);%-d(phi)/d(theta)1
deltar(h)=k3(h)*c/w/(pi/1800);
t2(h)=-t1(h+1)+t1(h);
deltat(h)=t2(h)*c/w/(pi/1800);
end
k3(900)=k3(899);%-d(phi)/d(theta)2
deltar(900)=k3(900)*c/w/(pi/1800);
t2(900)=t2(899);
deltat(900)=t2(900)*c/w/(pi/1800);
c5=0:0.1:89.9;
subplot(2,1,1)
plot(c5,deltar);
subplot(2,1,2);
plot(c5,T,c5,R,'r');