%%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 nx, ny
n1=(-10:10);
n2=(-220:220);
%Form the matrix indices gg' 
%   g=row index                 nx,ny
%   g'=column index             nx',ny'

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

%%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
pctb=2.79e11;        %dyn/cm2
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

kx=1;
ky=.5;

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
        subplot(1,3,1)
        b=plot(karr,omega,'*');
        po=get(gca,'position');
        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,2.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,2.5]);
        left_side=po(1);
        width=po(3);
        right_side=left_side+width;
    end
    
end