I would like to modify the equation to calculate total reflection coefficient for N multilayer structure

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;

Answer 1

[Chris]

You could start with an approach like this (numbers are made up):

ep = [1 3.9 12 1].';  % The apostrophe transposes the array to vertical. 
% The .' specifically, ensures it's not a conjugate transpose, in case you're 
% working with complex numbers.
% Similarly, .* ./ and .^ are element-wise operations, as opposed to matrix
% math.
kx = [0 1e-6 .03 0].';
w = [1 .001 .1 1e10].'; % Close enough to infinite?
% What we have so far
table(ep, kx, w)

ans = 4×3 table
    ep      kx        w  
    ___    _____    _____

      1        0        1
    3.9    1e-06    0.001
     12     0.03      0.1
      1        0    1e+10

c0 = 3;
% Now the math. Note Small numbers look like zero when displayed with larger
% numbers
kz = sqrt(ep.*w.^2./c0^2 -kx.^2)

kz = 4×1

1.0e+09 * 0.0000 0.0000 0.0000 3.3333

rs = (kz(1:end-1)-kz(2:end))./(kz(1:end-1)+kz(2:end))

rs = 3×1

0.9961 -0.9883 -1.0000

rp = (ep(2:end).*kz(1:end-1)-kz(2:end))./(ep(2:end).*kz(1:end-1)-kz(2:end))

rp = 3×1

1 1 1

Does this give you an idea how to do it?

[Ambali Odebowale]

I still didn't get it

[Chris]

Could you post what you tried?