x3=2.09;
for x2=8:12
x1=1./(1./x3-1./x2);
n1o=(4.2243+1.768*x1.^2./(x1.^2-0.227)+3.12*x1.^2./(x1.^2-3380)).^0.5;
n1e=(4.2009+1.8875*x1.^2./(x1.^2-0.2171)+3.6461*x1.^2./(x1.^2-3629)).^0.5;
n2o=(4.2243+1.768*x2.^2./(x2.^2-0.227)+3.12*x2.^2./(x2.^2-3380)).^0.5;
n2e=(4.2009+1.8875*x2.^2./(x2.^2-0.2171)+3.6461*x2.^2./(x2.^2-3629)).^0.5;
n3o=(4.2243+1.768*x3.^2./(x3.^2-0.227)+3.12*x3.^2./(x3.^2-3380)).^0.5;
w1=2*pi*3*10.^8./x1;
w2=2*pi*3*10.^8./x2;
w3=2*pi*3*10.^8./x3;
tm=asin((n2e.^2*((w3*n3o-w1*n1o).^2-w2.^2*n2o.^2)./((w3*n3o-w1*n1o).^2*(n2e.^2-n2o.^2))).^0.5);
t=atan(0.5*(n2e.^2-n1o.^2)*sin(2*tm)./((n1o*sin(tm)).^2+(n2e*cos(tm)).^2));
plot(x2,t);
hold on
end
for x2=8:12
x1=1./(1./x3-1./x2);
n1o=(4.2243+1.768*x1.^2./(x1.^2-0.227)+3.12*x1.^2./(x1.^2-3380)).^0.5;
n1e=(4.2009+1.8875*x1.^2./(x1.^2-0.2171)+3.6461*x1.^2./(x1.^2-3629)).^0.5;
n2o=(4.2243+1.768*x2.^2./(x2.^2-0.227)+3.12*x2.^2./(x2.^2-3380)).^0.5;
n2e=(4.2009+1.8875*x2.^2./(x2.^2-0.2171)+3.6461*x2.^2./(x2.^2-3629)).^0.5;
n3o=(4.2243+1.768*x3.^2./(x3.^2-0.227)+3.12*x3.^2./(x3.^2-3380)).^0.5;
w1=2*pi*3*10.^8./x1;
w2=2*pi*3*10.^8./x2;
w3=2*pi*3*10.^8./x3;
tm=asin((n2e.^2*((w3*n3o-w1*n1o).^2-w2.^2*n2o.^2)./((w3*n3o-w1*n1o).^2*(n2e.^2-n2o.^2))).^0.5);
t=atan(0.5*(n2e.^2-n1o.^2)*sin(2*tm)./((n1o*sin(tm)).^2+(n2e*cos(tm)).^2));
plot(x2,t);
hold on
end