matlab光栅仿真

neff=3.5; k0=3000; L=20e-3; p=0.1e-3; q=0.4e-3; M=20; dz=p/M; d=0.5e-3; pe=221e-9; wl=1533e-9; l=L/500; z=zeros(1,20); B=zeros(2,2,20); B1=zeros(2,2,40); B2=zeros(2,2,40); C=zeros(1,300); z=0:0.1e-4:0.2e-3; for x=1:300 wl=wl+1e-10; cc=2*pi*neff/wl; dc=cc-pi/pe; dcc=2*dc; for t=1:40 if mod(t,2)==1 for m=1:20 k=k0*exp(-4*(m*dz-p/2)^2/(p^2)); sm=sqrt((k^2-dc^2)); T1=[cosh(sm*dz)+(i*dc/sm)*sinh(sm*dz)]*exp(-i*dc*dz); T2=i*(k/sm)*sinh(sm*dz)*exp(-i*dc*(z(m)+z(m+1))); T3=-(i*k/sm)*sinh(sm*dz)*exp(i*dc*(z(m)+z(m+1))); T4=[cosh(sm*dz)-(i*dc/sm)*sinh(sm*dz)]*exp(i*dc*dz); if m==1 B(:,:,m)=[T1 T2;T3 T4]; else B(:,:,m)=B(:,:,m-1)*[T1,T2;T3,T4]; end end B1(:,:,t)=B(:,:,20); else B1(:,:,t)=[exp(cc*q*i),0;0,exp(-cc*q*i)]; end if t==1 B2(:,:,t)=B1(:,:,t); else B2(:,:,t)=B2(:,:,t-1)*B1(:,:,t); end end A=B2(:,:,40); %C(x)=(abs(-A(2,1)*A(1,2)/A(2,2)+A(1,1)))^2;%投射谱 C(x)=(abs(A(2,1)/A(2,2)))^2;%反射谱 end x=1:300; %subplot(1,2,1); %plot(1550e-9+t*0.1e-9,C(t)); %subplot(1,2,2); plot(1533e-9+x*(1e-10),C(x));