%多層結構最終反射系數曲線計算

n0=1.65;%棱鏡折射率
e0=n0^2;%棱鏡介電常數

n1=1.56;%波片 
e1=n1^2;%波片介電常數

e2=-13.4+1.4i;%金屬膜介電常數

n3=1.35;%1.40;1.45;1.50;%被測物折射率
e3=n3^2;

n4=1.33;%緩沖液體
e4=n4^2;

lamta=632.8;%波長 

theta_degree=linspace(40,70,100);%入射光度數離散化
theta_deg=pi/180*theta_degree;  %化為弧度

for number=1:100;
theta=theta_deg(number);

kz0=2*pi/lamta*sqrt(e0-e0*sin(theta)^2);
kz1=2*pi/lamta*sqrt(e1-e0*sin(theta)^2);
kz2=2*pi/lamta*sqrt(e2-e0*sin(theta)^2);
kz3=2*pi/lamta*sqrt(e3-e0*sin(theta)^2);
kz4=2*pi/lamta*sqrt(e4-e0*sin(theta)^2);

zeta0=e0/kz0;
zeta1=e1/kz1;
zeta2=e2/kz2;
zeta3=e3/kz3;
zeta4=e4/kz4;


r34=(zeta4-zeta3)/(zeta4+zeta3);
r23=(zeta3-zeta2)/(zeta3+zeta2);
r12=(zeta2-zeta1)/(zeta2+zeta1);
r01=(zeta1-zeta0)/(zeta1+zeta0);


d3=0;  %被測物厚度
d2=60; %金屬膜厚度
d1=40; %波片厚度


r24=(r23+r34*exp(2*i*d3*kz3))/(1+r23*r34*exp(2*i*d3*kz3));
r14=(r12+r24*exp(2*i*d2*kz2))/(1+r12*r24*exp(2*i*d2*kz2));
r04=(r01+r14*exp(2*i*d1*kz1))/(1+r01*r14*exp(2*i*d1*kz1));

r=abs(r04)^2;

R(number)=r;%最終反射系數

end

plot(theta_degree,R,'r')