clear,close all
na=3.1693;
nb=3.3884;
nc=3.1693;
nd=1;
ne=3.1693;
nf=1;
ng=1;
n=[na nb nc nd ne nf ng]
lengh=1.3;
k0=2*pi/lengh;
%%%%%%%%%%%%%%%%%%%%% 網格劃分 %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% 非等間距劃分 橫向 31點、縱向81點
deltx=0.2;               % 橫向 31 點 
delty1=0.24;            % 縱向 81點
delty2=0.04;
delty3=0.015;
delty4=0.01;
delty5=0.01;
delty6=0.14;
delty7=0.01;
delty8=0.095
             
y1=1/delty1;
y2=1/delty2;
y3=1/delty3;
y4=1/delty4;
y5=1/delty5;
y6=1/delty6;
y7=1/delty7;
y8=1/delty8;
x=1/deltx;
for i=1:31              % 給每個格點的折射率賦值
    for j=1:81
        if j<=16
            N(j,i)=na;
        elseif j>16&j<=36
            N(j,i)=nb;
        elseif j>36&j<=46
            N(j,i)=nc;
        elseif j>46&j<=66&i>1&i<=11
            N(j,i)=nd;
        elseif j>46&j<=66&i>11&i<=21
            N(j,i)=ne;
        elseif j>46&j<=66&i>21&i<=31
            N(j,i)=nf;
        else N(j,i)=ng
        end
    end
end
 
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 構造本征方程的矩陣 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CEX=zeros(31*81,31*81);        % 創建一個矩陣，就是系數矩陣
M=81
k=1
for m=1:31
   for n=1:81
                        
% x=(n-1)*M+m 格點數與坐標的轉換關系 即N(n,m)=N((n-1)*M+m)
% 周圍四個點的轉換

a2=(m-1)*M+(n-1);     % m,n-1
a1=((m-1)-1)*M+n;     % m-1,n
a3=(m-1)*M+n;         % m,n
a5=((m+1)-1)*M+n;     % m+1,n
a4=(m-1)*M+(n+1);     % m,n+1
    if (n-1)<=0
        a2=0;
          end
    if (n+1)>81
        a4=0;
          end
    if (m-1)<=0
        a1=0;
          end   
    if (m+1)>31
        a5=0;
          end

if a1~=0    
Tered=2*N(a1)/(N(a1)+N(a3));
else Tered=0;
end
if a5~=0
Teadd=2*N(a5)/(N(a5)+N(a3));
else Teadd=0;
end

if n<=11
    y=y1;
elseif n>11&n<=21
    y=y2;
elseif n>21&n<=31
    y=y3;
elseif n>31&n<=41
    y=y4;
elseif n>41&n<=51
    y=y5;
elseif n>51&n<=61
    y=y6;
elseif n>61&n<=71
    y=y7;
else
    y=y8;
end
Cx=[x^2*Tered y^2 -2*y^2-(4-Teadd-Tered)*x^2+k0^2*na^2 y^2 x^2*Teadd ];
 
%%%%%%%%%%%給矩陣賦值    
    if a1~=0
        CEX(k,a1)=Cx(1);
         end    
    if a2~=0
         CEX(k,a2)=Cx(2);
         end    
    if a3~=0
             CEX(k,a3)=Cx(3);
         end    
    if a4~=0
                 CEX(k,a4)=Cx(4);
         end    
     if a5~=0
                     CEX(k,a5)=Cx(5);
         end
        
k=k+1;
end
end


%%%%%%%%%%%%%%%%%%%%%%%%  求本征值和本征向量 %%%%%%%%%%%%%%%%%%%%%%%%%%
[V,D]=eig(CEX);      % V 代表本征向量構成的矩陣，D 代表本征值構成的對角矩陣
%  取出本征值，求出其中的最大數，作為所求的準TE模的本征值,也即傳播常數
B=diag(D);           
b=B(1);
for i=1:2511
    y=B(i)-b;
    if y>0
        b=B(i);
        k=i
    end
end
neff=sqrt(b/k0^2)
E=V(:,k);
k=1
for i=1:31
for j=1:81
EX(i,j)=abs(E(k));
k=k+1
end
end
% 構造本征向量矩陣，矩陣中的點的值對應與分割網格上的電場的值
%%%%%%% 繪圖 %%%%%%%%%





%%%%%%%%%  %%%%%%%%%
%Cex*Ex=b^2*Ex