I have attached an equation for calculating total reflection coefficient for multilayer structure. I have used the equation to calculate total reflection coefficient for five layer structure (made of metal and dielectrics). However, I would like to modify the code to be used for N layer structure. Below is the matlab I have written for 5 layers; the fift layer is semi-infinite. The thickness of the metal and dielectrics are denoted by dsm and ds respectively.
for ind=1:n+1
kx(ind)= min+(ind-1)*dkx;
% vertical component of wavevector k
kz1 = sqrt(ep_1*w*w/(c0*c0)-kx(ind)^2);
kz2 = sqrt(ep_2*w*w/(c0*c0)-kx(ind)^2);
kz3 = sqrt(ep_3*w*w/(c0*c0)-kx(ind)^2);
kz4 = sqrt(ep_4*w*w/(c0*c0)-kx(ind)^2);
kz5 = sqrt(ep_5*w*w/(c0*c0)-kx(ind)^2);
kz0 = sqrt(w*w/(c0*c0)-kx(ind)^2); %vacuum gap
% frensnel coefficients
rs01 = ((kz0-kz1)/(kz0+kz1)); %s-polarization interface between 1 and 0
rs12 = ((kz1-kz2)/(kz1+kz2)); %s-polarization interface between 2 and 1
rs23 = ((kz2-kz3)/(kz2+kz3));
rs34 = ((kz3-kz4)/(kz3+kz4));
rs45 = ((kz4-kz5)/(kz4+kz5));
rp01 = ((ep_1*kz0-kz1)/(ep_1*kz0+kz1)); %p-polarization interface between 1 and 0
rp12 = ((ep_2*kz1-kz2)/(ep_2*kz1+kz2)); %p-polarization interface between 2 and 0
rp23 = ((ep_3*kz2-kz3)/(ep_3*kz2+kz3));
rp34 = ((ep_4*kz3-kz4)/(ep_4*kz3+kz4));
rp45 = ((ep_5*kz4-kz5)/(ep_5*kz4+kz5));
%total reflection
%s polarization
rs56=rs45;
rs45=(rs45+rs56*exp(1i*2*kz5*ds))/(1+rs45*rs56*exp(1i*2*kz5*ds));
rs34=(rs34+rs34*exp(1i*2*kz4*dsm))/(1+rs34*rs45*exp(1i*2*kz4*dsm));
rs23=(rs23+rs34*exp(1i*2*kz3*ds))/(1+rs23*rs34*exp(1i*2*kz3*ds));
rs12=(rs12+rs23*exp(1i*2*kz2*dsm))/(1+rs12*rs23*exp(1i*2*kz2*dsm));
rs_total=rs56+rs45+rs34+rs23+rs12+rs01;
rp56=rp45;
rp45=(rp45+rp56*exp(1i*2*kz5*ds))/(1+rp45*rs56*exp(1i*2*kz5*ds));
rp34=(rp34+rp34*exp(1i*2*kz4*dsm))/(1+rp34*rp45*exp(1i*2*kz4*dsm));
rp23=(rp23+rp34*exp(1i*2*kz3*ds))/(1+rp23*rp34*exp(1i*2*kz3*ds));
rp12=(rp12+rp23*exp(1i*2*kz2*dsm))/(1+rp12*rp23*exp(1i*2*kz2*dsm));
rp_total=rp56+rp45+rp34+rp23+rp12+rp01;
