clear;clc;close?all;
sz_X=41;sz_Y=23;
sz_M=sz_X-2;sz_N=sz_Y-2;
h=0.254;
%?sz_X=81;sz_Y=45;
%?sz_M=sz_X-2;sz_N=sz_Y-2;
%?h=0.127
%?sz_X=21;sz_Y=12;
%?sz_M=sz_X-2;sz_N=sz_Y-2;
%?h=0.508;
%TM波計(jì)算
DM=diag（4*ones（sz_M1）0）+diag（-ones（sz_M-11）1）+diag（-ones（sz_M-11）-1）;
for?i=1:sz_N
????D1（::i）=DM;
end
?sz_D1?=?size（D1）;
c?=?mat2cell（reshape（D1sz_D1（1）[]）sz_D1（1）repmat（sz_D1（2）1sz_D1（3）））;
?K1=feval（@blkdiagc{:}）;
?KM=K1+diag（-ones（1sz_M*sz_N-sz_M）sz_M）+diag（-ones（1sz_M*sz_N-sz_M）-sz_M）;
?[v_TMd_TM]=eig（KM）;
?[dx_TMdy_TM]=find（d_TM>0）;
?l_dxTM=length（dx_TM）;
?for?i=1:l_dxTM
?????d2_TM（i）=d_TM（dx_TM（i）dy_TM（i））;
?end
?for?i=1:l_dxTM-1
?????minIndex=i;
?????for?j=i+1:l_dxTM
?????????if（d2_TM（j）?????????????minIndex=j;
?????????end
?????end
?????temp=d2_TM（i）;
?????d2_TM（i）=d2_TM（minIndex）;
?????d2_TM（minIndex）=temp;
?end
?kc_TM=sqrt（d2_TM）/h;
?%TE波計(jì)算
?DE=diag（4*ones（sz_X1）0）+diag（-ones（sz_X-11）1）+diag（-ones（sz_X-11）-1）;
?DE（12）=-2;DE（sz_Xsz_X-1）=-2;
for?i=1:sz_Y
????D2（::i）=DE;
end
?sz_D2?=?size（D2）;
c?=?mat2cell（reshape（D2sz_D2（1）[]）sz_D2（1）repmat（sz_D2（2）1sz_D2（3）））;
?K1=feval（@blkdiagc{:}）;
?KE=K1+diag（-1*ones（1sz_X*sz_Y-sz_X）sz_X）+diag（-1*ones（1sz_X*sz_Y-sz_X）-sz_X）;
?for?i=1:sz_X
?????KE（ii+sz_X）=-2;
?????KE（i+sz_X*（sz_Y-1）i+sz_X*（sz_Y-2））=-2;?%KE矩陣
?end
?[v_TEd_TE]=eig（KE）;
?[dx_TEdy_TE]=find（d_TE>1e-6）;
?l_dxTE=length（dx_TE）;
?for?i=1:l_dxTE
?????d2_TE（i）=d_TE（dx_TE（i）dy_TE（i））;
?end
?for?i=1:l_dxTE-1
?????minIndex=i;
?????for?j=i+1:l_dxTE
?????????if（d2_TE（j）