clear all
syms k w p;
n1=1.0;
n2=3.2;
d1=1; % 850nm
d2=1;
d=d1+d2;
eps1=n1^2;
eps2=n2^2;
k_0=eps1*d1/d + eps2*d2/d;
kappa=-i*( (eps1-eps2)*(exp(i*2*pi*d1*p/d)-1) )/(2*pi*p);
k_p1=subs(kappa,1);
k_n1=subs(kappa,-1);
G0=subs(2*pi*p/d,1);
c=3*10^8;
perme=1/(c^2);
A=[ (k-G0)^2-w^2*perme*k_0 ,-w^2*perme*k_n1 , 0;
-w^2*perme*k_p1 ,k^2-w^2*perme*k_0 ,-w^2*perme*k_n1;
0 ,-w^2*perme*k_p1 ,(k+G0)^2-w^2*perme*k_0 ];
y=solve(det(A),w);
k1= 0:0.02:pi/d;
y1=subs(y(1),k,k1);
plot(k1/G0,y1/(G0*c),'R-o');
hold on
% y2=subs(y(2),k,k1);
% plot(k1/G0,y2/(G0*c),'k');
y3=subs(y(3),k,k1);
plot(k1/G0,y3/(G0*c),'b-*');
hold on
%y4=subs(y(4),k,k1);
%plot(k1/G0,y4/(G0*c),'R*');
y5=subs(y(5),k,k1);
plot(k1/G0,y5/(G0*c),'k>-');
%y6=subs(y(6),k,k1);
%plot(k1/G0,y6/(G0*c),'bo-');
grid on;
%title({'\epsilon_1/\epsilon_2=1/13; d_1/d_2=3/2'},'Color','b','fontsize',15);
title({'n1/n2=1.0/3.2; d_1/d_2=1/1'},'Color','b','fontsize',15);
xlabel('Wave vetor (kd/2\pi)','fontsize',12);
ylabel('Frequency (\omegad/2\piC)','fontsize',12);
%---------------------計算Vp Vg------------------
kz2=pi/(2*d)
vp_11=subs(y(3),k,kz2)/kz2
t1=diff(y(3))
vg_11=subs(t1,k,kz2)
vp_12=subs(y(5),k,kz2)/kz2
t2=diff(y(5))
vg_12=subs(t2,k,kz2)
%---------------------------------------
kz3=pi/d
vp_31=subs(y(3),k,kz3)/kz3
vg_31=subs(t1,k,kz3)
vp_32=subs(y(5),k,kz3)/kz3
vg_32=subs(t2,k,kz3)
%---------------------------------------
kz4=3*pi/(2*d)
vp_41=subs(y(3),k,kz4)/kz4
vg_41=subs(t1,k,kz4)
vp_42=subs(y(5),k,kz4)/kz4
vg_42=subs(t2,k,kz4)
