%%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
num_eigenvals=10;              %Number of Eigenvalues to keep
nmax=9;                        %Num of plane waves in expansion
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'

%f==filling fraction=pi*ro^2/Ac
f=.906;          %filling fraction
lat_a=.1;       %Lattice spacing for frequency units
ro=.0418;        %Cylinder Radius

%%pa,pb are densities of materials
pa=7.27e3;             %kg/m3  Steel
pb=1.29;              %kg/m3  Air

delp=(pb/pa-1)/(f*pb/pa+1-f);

%%pcla,pCtb
va=6086;             %m/s     Steel
vb=340;              %m/s     Air
pcla=(va^2)*pa;      %dyn/cm2 Steel
pclb=(vb^2)*pb;         %dyn/cm2 Air
delt=(pclb/pcla-1)/(f*pclb/pcla+1-f);

n1inc_arr=repmat(n1inc,1,2*nmax+1);
n2inc_arr=round(n2inc/(2*nmax+1));

[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=4*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+delp.*fG_Gp.*k_p_g_dot_kp_p_g.*(1-kro_del_ggp);
        N=kro_del_ggp+fG_Gp.*delt.*(1-kro_del_ggp);

        A=N^-1*M;
        [a,b]=eig(A);
        g=find(b~=0);
        eigs=(sort(b(g)));
        eigs(1:num_eigenvals)
        k=repmat(kxt,1,num_eigenvals);
        min(b(g));
    
        
         if onoff==1
            omega=[((340)/lat_a)*sqrt(eigs(1:num_eigenvals))'];
            karr=[k];
	    if kyind==1
             omegafull=[((340)/lat_a)*sqrt(eigs(1:num_eigenvals))'];
             karrfull=[k];
	    end		
            onoff=0;  
        else
          omega=[omega,((340)/lat_a)*sqrt(eigs(1:num_eigenvals))'];
          karr=[karr,k];
	      omegafull=cat(1,omegafull,...
		      ((340)/lat_a)*sqrt(eigs(1:num_eigenvals))');
          karrfull=cat(1,karrfull,k);
        end        
    end
    
    if kyind==1
        ymax=max(omega);
        subplot(1,3,1)
        b1=plot(karr,omega,'*');
        ax1=gca; 
        po=get(gca,'position');
        set(gca,'YLim',[0,ymax]); 
        set(gca,'XLim',[0,.5]);
        set(gca,'nextplot','add');	
	    left_side=po(1);
        set(gca,'xtick',[]);
        width=po(3);
        right_side=left_side+width;
        set(gca,'xticklabel',{});
        ylabel('Frequency (kHz)');
    
  elseif kyind==2
        subplot(1,3,2)
        b2=plot(karr,omega,'*');
        po=get(gca,'position');
        po(1)=right_side;
 	    ax2=gca;
	    set(gca,'nextplot','add');
        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;
        set(gca,'xticklabel',{});
        
        
 elseif kyind==3
        subplot(1,3,3)
        b3=plot(karr,omega,'*');
        po=get(gca,'position')
        po(1)=right_side;
    	ax3=gca;
    	set(gca,'nextplot','add');
        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;
        set(gca,'xticklabel',{});
        text(-1,-.05*ymax,'M');
        text(-.5,-.05*ymax,'\Gamma');
        text(0,-.05*ymax,'X');
        text(.5,-.05*ymax,'M');
        text(-.3,-.1*ymax,'Reduced Wave Vector','HorizontalAlignment','center');
	    text(-.3,1.075*ymax,['Dispersion Relation for Triangular Lattice of Steel Cylinders in Air'],...
                             'HorizontalAlignment','Center');
        text(-.3,1.025*ymax,['f=',num2str(f),'            a=',num2str(lat_a), 'm           r_{o}=',...
                              num2str(ro), '           N=',num2str(nmax)],...
                              'HorizontalAlignment','center');
	%Plot the Bandgaps
	for evals=2:num_eigenvals
	  if min(omegafull(:,evals))-.025*ymax > max(omegafull(:,evals-1))
	    top=min(omegafull(:,evals))-.01*ymax;
		bottom=max(omegafull(:,evals-1))+.01*ymax;
	  	xband=[0 0 .5 .5 0];
		yband=[bottom top top bottom bottom];
		bg1=area(xband,yband);
		bg2=area(xband,yband,'Parent',ax2);
		bg3=area(xband,yband,'Parent',ax1);
        bottom
        top
      end
  end
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