MATLAB仿真传输矩阵法(TMM)计算光学薄膜/一维光子晶体的透射谱/反射谱 可定制s波，p波，增透膜，高反膜，随入射角变化等。

传输矩阵法这玩意儿在光学薄膜设计里就像炒菜时的锅铲，虽然不起眼但没它真不行。今天咱们用MATLAB整点活，手把手教你怎么玩转多层膜系的透反射谱计算。先别急着掏公式，咱直接上代码更实在。

先来点基础设定：

lambda = 400:800; % 波长范围400-800nm
theta = 30; % 入射角
polarization = 's'; % 偏振类型
n = [1.0, 2.35, 1.38, 1.52]; % 介质折射率序列
thickness = [0, 100, 200, 0]; % 各层厚度(nm)

这个n数组里首尾是入射和出射介质，中间才是真正的膜层。厚度为0表示半无限介质，这设计挺巧妙吧？

核心在传输矩阵计算，咱们拆开看：

function M = TMM_layer(n, d, lambda, theta, pol)
    k0 = 2*pi./lambda;
    theta = deg2rad(theta);
    n_eff = n.*sqrt(1 - (n(1)*sin(theta)/n).^2); % 有效折射率
    
    delta = k0 .* n_eff * d;
    if pol == 's'
        m = [cos(delta), 1i*sin(delta)/n_eff;
             1i*n_eff.*sin(delta), cos(delta)];
    else % p偏振
        m = [cos(delta), 1i*sin(delta)/(n_eff.*n);
             1i*n_eff.*n.*sin(delta), cos(delta)];
    end
end

这里有个小trick：用矩阵分块直接构建传输矩阵，比逐个元素赋值优雅多了。注意p波处理时多乘的那个n，这可是边界条件决定的，手滑写错的话反射谱会亲妈都不认识。

整个系统的传输矩阵是各层矩阵连乘：

M_total = eye(2);
for layer = 2:length(n)-1
    M_layer = TMM_layer(n(layer), thickness(layer), lambda, theta, polarization);
    M_total = M_total * M_layer;
end

这个循环从第二层开始到倒数第二层结束，完美避开首尾的半无限介质层。矩阵连乘顺序别搞反了，光学传输是右乘链式结构，跟穿衣服顺序一样从里到外。

MATLAB仿真传输矩阵法(TMM)计算光学薄膜/一维光子晶体的透射谱/反射谱 可定制s波，p波，增透膜，高反膜，随入射角变化等。

最后算反射率透射率：

Y_in = n(1)*cos(theta); 
Y_out = n(end)*sqrt(1 - (n(1)*sin(theta)/n(end))^2);

r = (M_total(2,1) - Y_out*M_total(1,1) + Y_in*(Y_out*M_total(2,2) - M_total(1,2))) / ...
    (M_total(2,1) + Y_out*M_total(1,1) + Y_in*(Y_out*M_total(2,2) + M_total(1,2)));
R = abs(r).^2;
T = 1 - R;

这堆分式看着头疼？其实本质就是电磁场连续性条件的矩阵版。注意Yin/Yout的处理，不同偏振要对应不同的波阻抗表达式。

举个实战栗子——做个13层的TiO2/SiO2高反膜：

n_high = 2.35; % TiO2
n_low = 1.46;  % SiO2
n_stack = [1.0, repmat([n_high, n_low],1,6), n_high, 1.52];
thickness_stack = [0, 100*ones(1,13), 0]; 

% 跑计算
[R, T] = TMM_calculate(n_stack, thickness_stack, 0, 's'); 

% 画图
plot(lambda, R, 'r', lambda, T, 'b');

看到550nm附近反射率飙到99%以上没？这就是光子禁带的威力。如果换成增透膜结构，把高低折射率顺序倒过来，立马变透射率95%+。

想要玩入射角扫描？套个循环就完事：

angles = 0:80;
spectra = zeros(length(angles), length(lambda));
for i = 1:length(angles)
    [R, ~] = TMM_calculate(n_stack, thickness_stack, angles(i), 'p');
    spectra(i,:) = R;
end
surf(lambda, angles, spectra); % 三维彩图走起

这出来的彩虹色图谱，专业点的叫法是角度分辨反射谱。p偏振在布鲁斯特角附近会出个反射率凹陷，这是判别偏振态的活标本。

代码里藏着几个坑得提醒：

1. 波长别用linspace，直接冒号运算符更高效
2. 有效折射率计算要带n(1)的平方项，很多人漏这个导致大角度算崩
3. 矩阵乘法用MATLAB原生*号，别手贱写for循环乘——效率差十倍不止

最后甩个冷知识：TMM算光子晶体时，记得把结构周期数设到20层以上，否则禁带特征不明显。但层数超过50之后计算量会指数爆炸，这时候得换FDTD之类的时域算法了。