slab.m

Calculate the reflection and transmission coefficients using the S-Matrix approach in Matlab.

clc
clear all
close all
%%RCWA INPUT VARIABLES

input_var_ISKANDER;

kc=k0*nc;
kx_inc=kc*sin(theta0)*cos(phi0);
ky_inc=kc*sin(theta0)*sin(phi0);
kzc_inc=kc*cos(theta0);

%%Transmitted Medium
ks=k0*ns;
kzs_tra=sqrt(ks^2-kx_inc^2-ky_inc^2);
kzs_tra = real(kzs_tra) - i*abs(imag(kzs_tra));
kxy = (kx_inc*ky_inc/kzs_tra)*(-i/k0);
kx2 = (((kx_inc/k0)^2 - ns^2)/kzs_tra)*(-i*k0);
ky2 = (((ky_inc/k0)^2 - ns^2)./kzs_tra)*(i*k0);
A = [kxy ky2; kx2 -kxy];
mat2=[diag([1 1]) diag([1 1]); -A A];

%%Slab Medium
kr=k0*nr;
kzr_tra=sqrt(kr^2-kx_inc^2-ky_inc^2);
kzr_tra = real(kzr_tra) - i*abs(imag(kzr_tra));
kxy = (kx_inc*ky_inc/kzr_tra)*(-i/k0);
kx2 = (((kx_inc/k0)^2-nr^2)/kzr_tra)*(-i*k0);
ky2 = (((ky_inc/k0)^2-nr^2)/kzr_tra)*(i*k0);
A = [kxy ky2; kx2 -kxy];
mat1=[diag([1 1]) diag([1 1]); -A A];

AC1 = inv(mat1)*mat2;

A=AC1(1:2,1:2);
C=AC1(3:4,1:2);
kappa_L=([j*kzr_tra/k0 j*kzr_tra/k0]);
ft=inv(A)*diag(exp(-k0*kappa_L*depth));
b1=diag(exp(-k0*kappa_L*depth))*C*inv(A)*diag(exp(-k0*kappa_L*depth));

%%Incident Medium
kzc_inc = real(kzc_inc) - i*abs(imag(kzc_inc));
kxy = (kx_inc*ky_inc/kzc_inc)*(-i/k0);
kx2 = (((kx_inc/k0)^2-nc^2)/kzc_inc)*(-i*k0);
ky2 = (((ky_inc/k0)^2-nc^2)/kzc_inc)*(i*k0);
A = [kxy ky2; kx2 -kxy];
mat0=[diag([1 1]) diag([1 1]); -A A];

AC0 = inv(mat0)*mat1;
A1=AC0(1:2,1:2);
B1=AC0(1:2,3:4);
C1=AC0(3:4,1:2);
D1=AC0(3:4,3:4);

f1=inv((A1+B1*b1));
ft=ft*f1;
bi=(C1+D1*b1)*f1;

ux0=cos(phi0)*cos(theta0)*cos(psi)-sin(phi0)*sin(psi);
uy0=sin(phi0)*cos(theta0)*cos(psi)+cos(phi0)*sin(psi);
uz0=-sin(theta0)*cos(psi);
fi=[uy0;ux0];

R1 = bi*fi;
T1 = ft*fi;

Ry = R1(1)
Rx = R1(2)
Rz = (kx_inc*Rx + ky_inc*Ry)./kzc_inc;
    
Ty = T1(1);  
Tx = T1(2);
Tz = -(kx_inc*Tx + ky_inc*Ty)./kzs_tra;

n1=1;
for ide = 1:length(kzc_inc)       
    Rfm = abs(Rx(ide))^2 + abs(Ry(ide))^2 + abs(Rz(ide))^2;
    pr(ide) = (real(kzc_inc(ide))/kzc_inc(n1))*Rfm;
end

for ide = 1:length(kzs_tra)       
    Tfm = abs(Tx(ide))^2 + abs(Ty(ide))^2 + abs(Tz(ide))^2;
    pt(ide) = (real(kzs_tra(ide))/kzc_inc(n1))*Tfm;
end

SUM= pr+pt;