clear close all maxw=2.5e15; minw=1e11; Nw=2000; w=linspace(minw,maxw,Nw); kj=1; for j=1:Nw; lam0=2*pi*3e14./w(j); % vacuum wavelength in micron %% output and input materials nout=1.5; nin=1; k0out=2*pi*nout./lam0; k0in=2*pi*nin./lam0; Np=20; thetain=pi/4; ky=k0in*sin(thetain); thetaout=asin(ky./k0out) %% unit cell (AB)%%%%%%%% nA=2; % refractive index of material A dA=0.25; % thickness of layer A in micron k0A=2*pi*nA./lam0; % modulus of the wavevector in layer A ( micron^(-1)) nB=1.428; % refractive index of material B dB=0.35; % thickness of layer B in micron k0B=2*pi*nB./lam0; % modulus of the wavevector in layer B ( micron^(-1)) a=dA+dB; % unit cell thickness kxA=sqrt(k0A^2-ky^2); % x component of the wavevector in layer A ( micron^(-1)) kxB=sqrt(k0B^2-ky^2); % x component of the wavevector in layer B (micron^(-1)) kxin=sqrt(k0in^2-ky^2); kxout=sqrt(k0out^2-ky^2); %% Transfer matrix of the unit cell %%%% M=zeros(2,2); PB= [ exp(-i.*kxB.*dB) 0; 0 exp(i.*kxB.*dB)]; IBA=0.5*[ 1+(kxB/kxA) 1-(kxB/kxA); 1-(kxB/kxA) 1+(kxB/kxA)]; PA= [ exp(-i.*kxA.*dA) 0; 0 exp(i.*kxA.*dA)]; IAB=0.5*[ 1+(kxA/kxB) 1-(kxA/kxB); 1-(kxA/kxB) 1+(kxA/kxB)]; M=IAB*PA*IBA*PB; %% finite structure Iout=0.5*[ 1+(kxout/kxB) 1-(kxout/kxB); 1-(kxout/kxB) 1+(kxout/kxB)]; Iout1=0.5*[ 1+(kxout/kxA) 1-(kxout/kxA); 1-(kxout/kxA) 1+(kxout/kxA)]; Iin=0.5*[ 1+(kxA/kxin) 1-(kxA/kxin); 1-(kxA/kxin) 1+(kxA/kxin)]; Iin1=0.5*[ 1+(kxB/kxin) 1-(kxB/kxin); 1-(kxB/kxin) 1+(kxB/kxin)]; NM=(Iin*PA*IBA*PB*(M^(Np-1))*Iout); %NM=M^Np; %% eigen values %% Ka=acos(0.5*(M(1,1)+M(2,2)))./a; %K=0; %KY=0; %if abs(imag(Ka))<2; tt(kj)=0; if abs(real(Ka))>0; if abs(real(Ka))