%%This script will attempt to find the band gaps in a phononic crystal by
%%carrying out a Plane Wave Expansion of the various quantities

%The indices of the reciprocal lattice are n1, n2
nmax=10;
n1inc=(-nmax:nmax);            %G=(2pi/a)(nx x + ny y)
n2inc=(-floor((2*nmax+1)^2/2):floor((2*nmax+1)^2/2));

%The reciprcoal lattice of the triangular lattice is given by  G
%           G=(n1+n2)x+sqrt(3)(n1-n2)y


%Form the matrix indices gg' 
%   g=row index                 nx,ny
%   g'=column index             nx',ny'

%%pa,pb are densities of materials
pa=8.936;            %g/cm3  Ni
pb=2.697;            %g/cm3  Al
pa=7.67;             %g/cm3  Steel
pb=1.29;             %g/cm3  Air

f=.906;               %f==filling fraction=pi*ro^2/Ac
delp=(pa/pb-1)/(f*pa/pb+1-f);

%%pCta,pCtb
pcta=7.54e11;        %dyn/cm2 Ni
pctb=2.79e11;        %dyn/cm2 Al
pcta=4.6097e11;      %dyn/cm2 Steel
pctb=4.3860e9;       %dyn/cm2 Air
delt=(pcta/pctb-1)/(f*pcta/pctb+1-f);

n1inc_arr=repmat(n1inc,1,21);
n2inc_arr=round(n2inc/21);

[n1,n1p]=meshgrid(n1inc_arr,n1inc_arr);
[n2,n2p]=meshgrid(n2inc_arr,n2inc_arr);

%Set up the G vectors for each nx,ny
%G=(n1+n2).^2+sqrt(3)*(n1-n2).^2;
%Gp=(n1p+n2p).^2+sqrt(3)*(n1p-n2p)^2;
%To deal with divide by zero, we can set up 
fill=ones(1,length(n1inc_arr));
fillt=diag(fill);
G_Gp=sqrt(2*pi*f/sqrt(3)*(((n1-n1p)+(n2-n2p)).^2+3*((n1-n1p)-(n2-n2p)).^2))+fillt;

%fG=2*f*besselj(1,G)./G;
%fGp=2*f*besselj(1,Gp)./Gp;
fG_Gp=2*f*besselj(1,G_Gp)./G_Gp;
%Form wavevector kx,ky

mode=1;
%%Mode=1 will calculate the band structure and plot it along the
%%Gamma-M-directions

if mode==1 

    for kyind=1:3
    
    onoff=1;
    for kx=.05:.05:.5
    
        if kyind==1
            kxt=kx;
            kx=.55-kx;
            ky=sqrt(3)*kx;
        elseif kyind==2
            kxt=kx;
            ky=0;
        elseif kyind==3
            kxt=kx;
            ky=kx;
            kx=.5;
        end
        
        mag_k_p_g=((kx+n1+n2).^2+(ky+sqrt(3).*(n1-n2)).^2);
        k_p_g_dot_kp_p_g=(kx+n1+n2).*(kx+n1p+n2p)+...
                         (ky+sqrt(3)*(n1-n2)).*(ky+sqrt(3)*(n1'-n2'));
        kro_del_ggp=(n1p==n1) & (n2p==n2);

        M=mag_k_p_g.*kro_del_ggp+delt.*fG_Gp.*k_p_g_dot_kp_p_g.*(1-kro_del_ggp);
        N=kro_del_ggp+fG_Gp.*delp.*(1-kro_del_ggp);

        A=N^-1*M;
        [a,b]=eig(A);
        g=find(b~=0);
        eigs=(sort(b(g)));
        eigs(1:5);
        k=repmat(kxt,1,5)
        min(b(g));
    
        if onoff==1
            omega=[eigs(1:5)'];
            karr=[k];
            onoff=0;  
        else
            omega=[omega,eigs(1:5)'];
            karr=[karr,k];
        end        
        
    end
    
    if kyind==1
        ymax=max(omega);
        subplot(1,3,1)
        b=plot(karr,omega,'*');
        po=get(gca,'position');
        set(gca,'YLim',[0,ymax]);
        set(gca,'XLim',[0,.5]);
        left_side=po(1);
        width=po(3);
        right_side=left_side+width;
        
    elseif kyind==2
        subplot(1,3,2)
        b=plot(karr,omega,'*');
        po=get(gca,'position');
        po(1)=right_side;
        set(gca,'position',po);
        set(gca,'YLim',[0,ymax]); 
        set(gca,'XLim',[0,.5]);
        set(gca,'ytick',[]);
        left_side=po(1);
        width=po(3);
        right_side=left_side+width;
        
    elseif kyind==3
        subplot(1,3,3)
        b=plot(karr,omega,'*');
        po=get(gca,'position')
        po(1)=right_side;
        set(gca,'position',po);
        set(gca,'ytick',[]);
        set(gca,'XLim',[0,.5]);
        set(gca,'YLim',[0,ymax]);
        left_side=po(1);
        width=po(3);
        right_side=left_side+width;
    end
    
end



%%Mode=2 will plot the equifrequency surfaces 
elseif mode==2
    %Prepare Matrices to hold eigenfrequencies
    kxa=0:0.05:.5;
    kya=0:0.05:.5;
    [kxm,kym]=meshgrid(kxa,kya);
    omega=zeros(length(kxa),length(kya))
    onoff=1;
    
    for i=1:length(kxa)
        kx=kxa(i);
        for j=1:length(kya)
          ky=kya(j);  
          mag_k_p_g=((kx+nx).^2+(ky+ny).^2);
          k_p_g_dot_kp_p_g=(kx+nx).*(kx+nxp)+(ky+ny).*(ky+nyp);
          kro_del_ggp=(nxp==nx) & (nyp==ny);

          M=mag_k_p_g.*kro_del_ggp+delt.*fG_Gp.*k_p_g_dot_kp_p_g.*(1-kro_del_ggp);
          N=kro_del_ggp+fG_Gp.*delp.*(1-kro_del_ggp);

          A=N^-1*M;
          [a,b]=eig(A);
          g=find(b~=0);
          eigs=(sort(b(g)));
          min(b(g));
    
          omega(i,j)=eigs(2);
          end        
        
    end                     %end of ky loop
  end                       %end of kx loop
end                         %end of If mode==2

end                         %end of file
end