%%This script will attempt to find the band gaps in a phononic crystal by
%%carrying out a Plane Wave Expansion of the various quantities
%
%It solves for the displacement according to the prescription in 
%Phys Rev Lett 71 2022

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


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

%f==filling fraction=pi*ro^2/Ac
f=.25;

%%pa,pb are densities of materials
pa=8.936;            %g/cm3  Ni
pb=2.697;            %g/cm3  Al
delp=(pa/pb-1)/(f*pa/pb+1-f);

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

nx_arr=repmat(n1,1,21);
ny_arr=round(n2/21);

[nx,nxp]=meshgrid(nx_arr,nx_arr);
[ny,nyp]=meshgrid(ny_arr,ny_arr);

%Set up the G vectors for each nx,ny
%G=nx.^2+ny.^2;
%Gp=nxp.^2+nyp.^2;
%To deal with divide by zero, we can set up 
fill=ones(1,length(nx_arr));
fillt=diag(fill);
G_Gp=sqrt(4*pi*f*((nx-nxp).^2+(ny-nyp).^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=kx;
        elseif kyind==2
            kxt=kx;
            ky=0;
        elseif kyind==3
            kxt=kx;
            ky=.5;
        end
        
        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)));
        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);
        set(gca,'ytick',[]);
        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,'YLim',[0,ymax]); 
        set(gca,'XLim',[0,.5]);
        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