te.m

利用平面波展开发计算二维光子晶体的能带结构

function TE
filling_ratio=0.6;
planewave_number=3;
dc_background=1;
dc_dielectric=13.6;
twopi=2*pi;
d_lattice_constant=1;
%triangular
lattice_1(1)=d_lattice_constant;
lattice_1(2)=0;
lattice_2(1)=d_lattice_constant/2;
lattice_2(2)=(d_lattice_constant/2)*3^0.5;
%calculating the wscell area of this two dimensional pc
	wscell=lattice_1(1)*lattice_2(2)-lattice_2(1)*lattice_1(2);
%calculating the radius of the dielectric circular rod
	R=(filling_ratio*wscell/pi)^0.5;
%calculating the reciprocal lattice vector
  	lattice_f1(1)=twopi*(lattice_2(2))/wscell;
	lattice_f1(2)=twopi*(-lattice_2(1))/wscell;
	lattice_f2(1)=twopi*(-lattice_1(2))/wscell;
	lattice_f2(2)=twopi*(lattice_1(1))/wscell;
%defining the bloch vector of main crystal line
% main crystal line for triangular lattice
	Bloch_vector_k_crystalline(2,1)=0;
	Bloch_vector_k_crystalline(2,2)=0;
	Bloch_vector_k_crystalline(1,1)=0;
	Bloch_vector_k_crystalline(1,2)=0.5*lattice_f2(2);
	Bloch_vector_k_crystalline(3,1)=Bloch_vector_k_crystalline(1,2)*tan(pi/6);
	Bloch_vector_k_crystalline(3,2)=Bloch_vector_k_crystalline(1,2);
	Bloch_vector_k_crystalline(4,1)=0;
	Bloch_vector_k_crystalline(4,2)=0.5*lattice_f2(2);

    planewave_setting_number=0;
    
    %setting the planewavelist
    for i=-planewave_number:1:planewave_number
  		for j=-planewave_number:1:planewave_number
			planewave_setting_number=planewave_setting_number+1;
			for s=1:1:2
				planewavelist(planewave_setting_number,s)=i*lattice_f1(s)+j*lattice_f2(s);
			end  
        end 
	end   
%allocate the memory for coefficient_matrix_B
for i=1:1:((2*planewave_number+1)^2)
    for j=1:1:((2*planewave_number+1)^2)
        temp_planewavelist(1)=planewavelist(i,1)-planewavelist(j,1);
		temp_planewavelist(2)=planewavelist(i,2)-planewavelist(j,2);
        %the fourier coefficient of dielectric constant of triangular lattice with circular rods
        if(abs(temp_planewavelist(1))==0&abs(temp_planewavelist(2))<1.0E-10)
            % triangular lattice with circular rod
            coefficient_matrix_B(i,j)=dc_background+filling_ratio*(dc_dielectric-dc_background);
        else
            %triangular lattice with circular rod 
		    temp_abs_planewavelist=(temp_planewavelist(1)^2+temp_planewavelist(2)^2)^0.5;
            coefficient_matrix_B(i,j)=2*filling_ratio*(dc_dielectric-dc_background)*besselj(1,temp_abs_planewavelist*R)/(temp_abs_planewavelist*R);
        end
    end
end
%dicating the number of main crystalline 
x=0;
for h=1:1:3
    for s=1:1:21
        %for triangular and square		
		Bloch_vector_k(1)=((s-1.0)/20.0)*(Bloch_vector_k_crystalline(h+1,1)-Bloch_vector_k_crystalline(h,1))+Bloch_vector_k_crystalline(h,1);
		Bloch_vector_k(2)=((s-1.0)/20.0)*(Bloch_vector_k_crystalline(h+1,2)-Bloch_vector_k_crystalline(h,2))+Bloch_vector_k_crystalline(h,2);   	
	%constructing the coefficient matrix A of the equation AX=bBX
			for i=1:1:(2*planewave_number+1)^2
				for j=1:1:(2*planewave_number+1)^2
					if (i==j)
						temp_planewavelist(1)=Bloch_vector_k(1)+planewavelist(i,1);
						temp_planewavelist(2)=Bloch_vector_k(2)+planewavelist(i,2);
                        coefficient_matrix_A(i,j)=temp_planewavelist(1)^2+temp_planewavelist(2)^2;
                    else
                        coefficient_matrix_A(i,j)=0 ;
					end 
				end 
			end          
            %calculating the coefficient matrix B^(-1)*A of the equation Ax=bBx
            coefficient_matrix_equation=inv(coefficient_matrix_B)*coefficient_matrix_A;  
            %solve the equation B^(-1)*Ax=bx
            [V,W]=eig(coefficient_matrix_equation);
            t=0;
            for p=1:1:(2*planewave_number+1)^2
                for q=1:1:(2*planewave_number+1)^2
                    if p==q
                        t=t+1;
                        eigenvalue(t)=W(p,q);
                    end
                end
            end
            %normalizing the eigenvalue by twopi*C/lattice_constant	
	        eigenvalue_flag=(2*planewave_number+1)^2;
			for i=(2*planewave_number+1)^2:-1:1
				eigenvalue_imag=imag(eigenvalue(i));
				if((abs(eigenvalue_imag)==0)&(real(eigenvalue(i))>=0))				
					eigenvalue_temp(eigenvalue_flag)=(real(eigenvalue(i)))^0.5*lattice_1(1)/twopi;
				    eigenvalue_flag=eigenvalue_flag-1;
				end 
			end 
			for i=eigenvalue_flag:-1:1
				eigenvalue_temp(i)=9e+999;
			end 
             
            % for oblique lattice
			Bloch_vector_k_temp(1)=Bloch_vector_k(1)*lattice_1(1)/twopi;
			Bloch_vector_k_temp(2)=Bloch_vector_k(2)*lattice_1(1)/twopi;   
            if s~=21
               x=x+1;
            elseif (s==21&h==3)
                 x=x+1;
             else
                 x=x;
            end
           for z=1:1:(2*planewave_number+1)^2
               u(z,x)=eigenvalue_temp(z);
           end
        
       end
   end
  %sort program
   u1=sort(u);
   dc=[0:1:60];
 for wc=1:1:6
       plot(dc,u1(wc,:),'r');
       
       hold on;
  
end