clc
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wl=550;     Z=0.001/100;              % Wave length
dt=1;                                 % Interval
D1=1000;    n(1)=1.0;                 % Air
D2=40;      n(2)=2.5495+0.016183i;    % MoO_3
D3=20;      n(3)=0.059582+3.5974i;    % sliver (Ag)
D4=100;     n(4)=n(2);                % MoO_3
D5=1000;    n(5)=1.5;                 % Glass
D6=D1;      n(6)=n(1);                % Air
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

d1=D1;  d2=D2;  d3=D3;  d4=0:dt:D4;  d5=D5;   d6=D6;

%Calculate Inter face Matrix j,j+1
for j=2:6
    t(j-1)=2*n(j-1.)/(n(j-1)+n(j));
    r(j-1)=(n(j-1)-n(j))./(n(j-1)+n(j));
end

I01=[1./t(1) r(1)./t(1);r(1)./t(1) 1./t(1)];
I12=[1./t(2) r(2)./t(2);r(2)./t(2) 1./t(2)];
I23=[1./t(3) r(3)./t(3);r(3)./t(3) 1./t(3)];
I34=[1./t(4) r(4)./t(4);r(4)./t(4) 1./t(4)];
I45=[1./t(5) r(5)./t(5);r(5)./t(5) 1./t(5)];

%Calculate Layer Matrix j

e1=-2*pi.*n(2).*d2*i./wl;
e2=-2*pi.*n(3).*d3*i./wl;

for j=1:D4/dt+1
    e3(j)=-2*pi.*n(4).*d4(j)*i./wl;
end

e4=-2*pi.*n(5).*d5*i./wl;
e5=-2*pi.*n(6).*d6*i./wl;
L2=[exp(e2) 0; 0 exp(-e2)];
L4=[exp(e5) 0; 0 exp(-e5)];

    for k=1:D4/dt+1
        L1=[exp(e1(j)) 0; 0 exp(-e1(j))];
        L3=[exp(e3(k)) 0; 0 exp(-e3(k))];
        S=I01*L1*I12*L2*I23*L3*I34*I45*L4;
        T(j,k)=(1./S(1,1))*(1./S(1,1))'*abs(n(6)./n(1));
        R(j,k)=(S(2,1)./S(1,1))*(S(2,1)./S(1,1))';
        A(j,k)=1-(T(j,k)+R(j,k));
    end

%Optimize
maxT=0;
    for k=1:D4/dt+1
        if maxT<(T(j,k)-Z)
            maxT=T(j,k);
            OR=R(j,k);
            OA=A(j,k);
            DT4=dt*(j-1);
        end
    end
    j=0:dt:D4;
maxT=maxT
OR=OR
OA=OA
DT4=DT4