Is there any way to plot graph between different values of volume fraction "f" and "freq"

clear all

warning off;

epsa=8.9; %material for the atom, alumina

epsb=1; %material of the background

a=1; %lattice constant

R=.2*a; %radius of the 'atom' or the radius of the cylinder

i=sqrt(-1);

f=pi*R^2/a^2; %volume fraction of atom(s) in the cell

NumberofCell=1; %supercell size=1X1

a1=a;

a2=a*i;

b1=2*pi/a/NumberofCell;

b2=2*pi/a/NumberofCell*i;

n=input('please input n:');

display('Fourier transforming.....');

tic;

NumberofPW=(2*n+1)^2;

mind=(-2*n:2*n)'+2*n+1;

mind=mind(:,ones([1,size(mind)]))-2*n-1;

GG=mind'*b1+mind*b2;

%clear mind;

eps21=2*f*(epsa-epsb)*besselj(1,abs(GG).*R)./(abs(GG).*R);

eps21(2*n+1,2*n+1)=epsb+f*(epsa-epsb);

%zz=[0,0]*[a1 a2].';%use 1X1 supercell to verify the algorithm

%5X5 super cell

zz=[

-2 -2;-2 -1;-2 0;-2 1;-2 2;

-1 -2;-1 -1;-1 0;-1 1;-1 2;

0 -2; 0 -1; 0 1; 0,2;

1 -2; 1 -1; 1 0; 1 1; 1 2;

2,-2; 2,-1; 2 0; 2,1; 2,2]*[a1 a2].';

zz=[0,0]*[a1 a2].';%use 1X1 supercell to verify the algorithm

%zz=[0 0; 0 1;1 0;1 1]*[a1 a2].';%2X2 supercell

%zz=[-1 -1;-1 0;-1 1;0 -1;0 0;0 1;1 -1;1 0;1 1]*[a1 a2].';

eps20=zeros(length(eps21));

for x=1:length(zz),

eps20=eps20+exp(i*(real(GG).*real(zz(x))+imag(GG).*imag(zz(x)))).*eps21;

end

ff=pi*R^2*length(zz)/(NumberofCell^2*a^2);

eps20=eps20./NumberofCell^2;

eps20(2*n+1,2*n+1)=epsb+ff*(epsa-epsb);

count=1;

for y=-n:n,

for x=-n:n;

G(count)=x*b1+y*b2;

r(count,:)=[x,y];

count=count+1;

end

end

display('Building eps(G,G) matrix from the FFT matrix');

for x=1:NumberofPW,

for y=x:NumberofPW,

b=r(x,:)-r(y,:)+2*n+1;

eps2(x,y)=eps20(b(2),b(1));

eps2(y,x)=eps2(x,y);

end

end

k1=2*pi/a*0.5.*(0:0.1:1);;

k2=2*pi/a*(0.5+(0.1:0.1:1).*0.5*i);

k3=2*pi/a*(0.5+0.5*i).*(0.9:-0.1:0);

k0=[k1 k2 k3]./NumberofCell;

display('Now calculate the eigen values..');

eps2=inv(eps2);

if max(max(real(eps2))) > 10^7*max(max(imag(eps2)))

display('Your lattice is inversion symmetric');

eps2=real(eps2);

else

display('Imaginary part of FFT is not zero');

stop;

%here we only demonstrate the inversion symmetric case

%%however, nonsymmetric case is also supported

end

options.disp=0;

counter=1;

for ii=1:length(k0),

k=k0(ii);

%k=k0;

M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*eps2; %TE wave

%M=abs(k+G.')*abs(k+G).*eps2; %TM wave

E=sort(abs(eig(M)));

%E=abs(eigs(M,20,'sm',options));%used to calculate only several smallest eigenvalues

freq(:,counter)=sqrt(abs(E)).*a./2./pi;

display(sprintf('calculation of k=%f+%fi is finished',real(k),imag(k)));

counter=counter+1;

end

toc;

[max(freq(1,:)),min(freq(2,:))]

tmpx=1:length(k0);

plot(tmpx,freq,'o-','linewidth',2)

title('TM:Band structure of a 2D square PBG')

xlabel('wave vector')

ylabel('wa/2\pic')

grid on

axis([1 31 0 1])