clcclear allclose all%%RCWA INPUT VARIABLESinput_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 Mediumks=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 Mediumkr=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 Mediumkzc_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;endfor 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;endSUM= pr+pt;