依据导纳矩阵法用matlab编写多层介质膜系反射率和透射率的计算程序

function [R,T] = multilayer_reflectance_transmittance(n,d,k,theta_in,wavelength)

% n: vector of refractive indices of each layer
% d: vector of thicknesses of each layer (in nm)
% k: vector of extinction coefficients of each layer
% theta_in: incidence angle (in radians)
% wavelength: wavelength of incident light (in nm)

c = physconst('LightSpeed');
lambda = wavelength*1e-9;

lambda_p = 2*pi*c./sqrt(eps0.*n.^2 - k.^2 - j*2*pi*eps0.*n.*k./(lambda/(2*pi*c))));
lambda_s = 2*pi*c./sqrt(eps0.*n.^2 - k.^2 + j*2*pi*eps0.*n.*k./(lambda/(2*pi*c)));

theta_t = asin(n(1)*sin(theta_in)/n(end)); % transmission angle

Z0 = 1/376.73031;

for i = 1:length(n)
[r_p(i),t_p(i)] = fresnel_am_tf(lambda_p(i),theta_in,n(i),n(i+1),k(i),k(i+1));
[r_s(i),t_s(i)] = fresnel_am_tf(lambda_s(i),theta_in,n(i),n(i+1),k(i),k(i+1));

Z(i) = n(i)*sqrt(eps0)/cos(theta_in);
Gamma_p(i) = (Z(i) - n(i+1)*sqrt(eps0)*cos(theta_t))/(Z(i) + n(i+1)*sqrt(eps0)*cos(theta_t));
Gamma_s(i) = (n(i)*sqrt(eps0)*cos(theta_t) - n(i+1)*sqrt(eps0)*cos(theta_in))/(n(i)*sqrt(eps0)*cos(theta_t) + n(i+1)*sqrt(eps0)*cos(theta_in));
Delta_p(i) = exp(-2*j*lambda_p(i)*d(i)/cos(theta_in));
Delta_s(i) = exp(-2*j*lambda_s(i)*d(i)/cos(theta_in));
end

r_pp = r_p(1);
r_ss = r_s(1);
t_pp = t_p(1);
t_ss = t_s(1);

for i = 2:length(n)
r_pp = r_pp + (t_pp^2 * Gamma_p(i-1) * r_p(i) * Delta_p(i-1));
r_ss = r_ss + (t_ss^2 * Gamma_s(i-1) * r_s(i) * Delta_s(i-1));
t_pp = t_pp * t_p(i) * Delta_p(i-1) * (1 + Gamma_p(i-1) * r_p(i));
t_ss = t_ss * t_s(i) * Delta_s(i-1) * (1 + Gamma_s(i-1) * r_s(i));
end

R_p = abs(r_pp)^2;
R_s = abs(r_ss)^2;
T_p = real((Z0 * cos(theta_t) / n(end)) * (abs(t_pp)^2));
T_s = real((Z0 * cos(theta_t) / n(end)) * (abs(t_ss)^2));

R = (R_p + R_s)/2;
T = (T_p + T_s)/2;

end

function [r_p,t_p,r_s,t_s] = fresnel_am_tf(lambda,theta_in,n1,n2,k1,k2)

% lambda: wavelength (in m)
% theta_in: incidence angle (in radians)
% n1, n2: refractive indices of incident and transmitted media, respectively
% k1, k2: extinction coefficients of incident and transmitted media, respectively

theta_t = asin(n1*sin(theta_in)/n2); % transmission angle

Z0 = 1/376.73031;

Gamma_p = (n1*cos(theta_in) - n2*cos(theta_t))/(n1*cos(theta_in) + n2*cos(theta_t));
Gamma_s = (n1*cos(theta_t) - n2*cos(theta_in))/(n1*cos(theta_t) + n2*cos(theta_in));

Delta_p = exp(-2*j*pi*n2*k2*lambda*cos(theta_t)/(lambda^2*n2^2*cos(theta_t)^2 + k2^2));
Delta_s = exp(-2*j*pi*n2*k2*lambda*cos(theta_t)/(lambda^2*n2^2*cos(theta_in)^2 + k2^2));

r_p = (Gamma_p + (1-Gamma_p)*Delta_p);
r_s = (Gamma_s + (1-Gamma_s)*Delta_s);
t_p = 2*n1*cos(theta_in)/(n1*cos(theta_in) + n2*cos(theta_t)) * (1-Gamma_p)*sqrt(Delta_p);
t_s = 2*n1*cos(theta_in)/(n1*cos(theta_t) + n2*cos(theta_in)) * (1-Gamma_s)*sqrt(Delta_s);

end