%% 1D photonic band off axis prop. isofrequency curves
clear
close all
maxw=2.5e15;
minw=1e11;
Nw=20;
w=linspace(minw,maxw,Nw);

for j=1:Nw;
kj=1;
lam0=2*pi*3e14./w(j); % vacuum wavelength in micron

%% unit cell  (AB)%%%%%%%%

nA=1.5;             % 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=3;     % refractive index of material B
dB=0.2;    % 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



kymin=-k0B;   %  wavevector component
kymax=k0B; %  max
kymed=k0A;
Nk=400*j;
kyv=[ linspace((kymin),(-kymed),Nk) linspace((-kymed),(kymed),Nk) linspace((kymed),(kymax),Nk)];
Nkk=length(kyv);
for k=1:Nkk
    ky=kyv(k);
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))

%% 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;

%% eigen values %%
Ka=acos(0.5*(M(1,1)+M(2,2)))./a;
%K=0;
%KY=0;
%if abs(imag(Ka))<2;
    if abs(real(Ka))>0;
        if abs(real(Ka))<pi./a;
K(kj)=real(Ka);
freq(kj)=w(j);
KY(kj)=ky;

K(kj+1)=-Ka;
freq(kj+1)=w(j);
KY(kj+1)=ky;

kj=kj+2;
end;
end;

end;
if kj>1
   % y=interp1(KY,K,kyv);
   
%% plot results   
hold off
AA=num2str(w(j),'% 5.2e');
TEX=['w=' AA]
plot(K*a./pi,KY.*a./pi,'.','MarkerSize',3)
hold on
plot(K*a./pi + 2,KY.*a./pi,'r.','MarkerSize',3)
plot(K*a./pi -2,KY.*a./pi,'r.','MarkerSize',3)

axis([-2 2, -2 2]);
legend(TEX)
grid on
xlabel('Kx (units of \pi/a)')
ylabel('Ky (units of \pi/a)')
pause
end
clear K KY kj
end;