%%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 is for an acoustic wave equation with pressure instead of the
%displacement as prescribed in App. Phys. Lett. 70 3218
%The indices of the reciprocal lattice are nx, ny
num_eigenvals=10;                %Number of Eigenvalues to keep
nmax=9;
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=0.9;          %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);

nx_arr=repmat(n1,1,2*nmax+1);
ny_arr=round(n2/(2*nmax+1));

[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;
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;
%fG_Gp=(2*pi*ro^2/lat_a^2)*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+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)*real(sqrt(eigs(1:num_eigenvals)))'];
            karr=[k];
	    if kyind==1
             omegafull=[((340)/lat_a)*real(sqrt(eigs(1:num_eigenvals)))'];
             karrfull=[k];
	    end		
            onoff=0;  
        else
          omega=[omega,((340)/lat_a)*real(sqrt(eigs(1:num_eigenvals)))'];
          karr=[karr,k];
	  omegafull=cat(1,omegafull,...
		      ((340)/lat_a)*real(sqrt(eigs(1:num_eigenvals)))');
          karrfull=cat(1,karrfull,k);
        end        
        
      end
    
      if kyind==1
        myfig=figure;
        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 Square 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)],...
                              'HorizontalAlignment','center');
	
    %Plot the Bandgaps
	for evals=2:num_eigenvals
	  if min(real(omegafull(:,evals)))-.025*ymax > max(real(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

%Determine operating system
osroot=matlabroot;
if osroot(1)=='C'
   filename=['C:\Documents and Settings\Lance\Desktop\School\Qual\Square_Disp',...
                '_f_',num2str(f),'_nmax_',num2str(nmax),'_a_',num2str(lat_a),'.tif']
      saveas(myfig,filename,'tif');   
else
   filename=['/afs/slac/u/ki/lances/matlab/figures/Square_Disp',...
                '_f_',num2str(f),'_nmax_',num2str(nmax),'_a_',num2str(lat_a),'.tif']
      saveas(myfig,filename,'tif');
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

   
