Projet Cristal Photonique

main.m

%The program for the computation of the PhC photonic
%band structure for 2D PhC.
%Parameters of the structure are defined by the PhC
%period, elements radius, and by the permittivities
%of elements and background material.
%Input parameters: PhC period, radius of an
%element, permittivities
%Output data: The band structure of the 2D PhC.
clear all; close all; clf;
%The variable a defines the period of the structure.
%It influences on results in case of the frequency
%normalization.
a=1e-6; % in meters
%The variable r contains elements radius. Here it is
%defined as a part of the period.
r=0.9*a; % in meters
%The variable eps1 contains the information about
%the relative background material permittivity.
eps1=9;
%The variable eps2 contains the information about
%permittivity of the elements composing PhC. In our
%case permittivities are set to form the structure
%perforated membrane
eps2=1;
%The variable precis defines the number of k-vector
%points between high symmetry points
precis=5;
%The variable nG defines the number of plane waves.
%Since the number of plane waves cannot be
%arbitrary value and should be zero-symmetric, the
%total number of plane waves may be determined
%as (nG*2-1)^2
nG=4;
%The variable precisStruct defines contains the
%number of nodes of discretization mesh inside in
%unit cell discretization mesh elements.
precisStruct=30;
%Variable latticeType allow you to choose a type of lattice
%0 for square
%1 for hexagonal
latticeType=1;
%The following loop carries out the definition
%of the unit cell. The definition is being
%made in discreet form by setting values of inversed
%dielectric function to mesh nodes.
nx=1;
for countX=-a/2:a/precisStruct:a/2
    ny=1;
    for countY=-a/2:a/precisStruct:a/2
        %The following condition allows to define the circle
        %with of radius r
        if(sqrt(countX^2+countY^2)<r)
            %Setting the value of the inversed dielectric
            %function to the mesh node
            struct(nx,ny)=1/eps2;
            %Saving the node coordinate
            xSet(nx)=countX;
            ySet(ny)=countY;
        else
            struct(nx,ny)=1/eps1;
            xSet(nx)=countX;
            ySet(ny)=countY;
        end
        ny=ny+1;
    end
    nx=nx+1;
end
%The computation of the area of the mesh cell. It is
%necessary for further computation of the Fourier
%expansion coefficients.
dS=(a/precisStruct)^2;
%Forming 2D arrays of nods coordinates
xMesh=meshgrid(xSet(1:length(xSet)-1));
yMesh=meshgrid(ySet(1:length(ySet)-1))';
%Transforming values of inversed dielectric function
%for the convenience of further computation
structMesh=struct(1:length(xSet)-1,...
    1:length(ySet)-1)*dS/(max(xSet)-min(xSet))^2;
%Defining the k-path within the Brillouin zone. The
%k-path includes all high symmetry points and precis
%points between them

switch latticeType
    case 0 %Square lattice
        kx(1:precis+1)=0:pi/a/precis:pi/a; %Path GX
        kx(precis+2:precis+precis+1)=pi/a; %Path XM
        kx(precis+2+precis:precis+precis+1+precis)=...
            pi/a-pi/a/precis:-pi/a/precis:0; %Path MG
        
        ky(1:precis+1)=zeros(1,precis+1); %Path GX
        ky(precis+2:precis+precis+1)=...
            pi/a/precis:pi/a/precis:pi/a; %Path XM
        ky(precis+2+precis:precis+precis+1+precis)=...
            pi/a-pi/a/precis:-pi/a/precis:0; %Path MG
    case 1 %Hexagonal lattice
        %Defining the coordinates of the high-symmetry point lattice in
        %reciprocal coordinates
        % Gamma
        Gx = 0;
        Gy = 0;
        
        % M
        Mx = pi/sqrt(a);
        My = -pi/(sqrt(3)*a);
        
        % K
        Kx = 4*pi/(3*a);
        Ky = 0;
        
        %Path GM
        kx(1:precis+1)=linspace(Gx,Mx,precis+1); 
        ky(1:precis+1)=linspace(Gy,My,precis+1); 
        
        %Path MK
        kx(precis+1:2*precis)=linspace(Mx,Kx,precis); 
        ky(precis+1:2*precis)=linspace(My,Ky,precis); 
        
        %Path KG
        kx(2*precis-1:3*precis-2)=linspace(Kx,Gx,precis); 
        ky(2*precis-1:3*precis-2)=linspace(Ky,Gy,precis); 
end
%After the following loop, the variable numG will
%contain real number of plane waves used in the
%Fourier expansion
numG=1;
%The following loop forms the set of reciprocal
%lattice vectors.
for Gx=-nG*2*pi/a:2*pi/a:nG*2*pi/a
    for Gy=-nG*2*pi/a:2*pi/a:nG*2*pi/a
        G(numG,1)= Gx;
        G(numG,2)= Gy;
        numG=numG+1;
    end
end
%The next loop computes the Fourier expansion
%coefficients which will be used for matrix
%differential operator computation.
for countG=1:numG-1
    for countG1=1:numG-1
        CN2D_N(countG,countG1)=sum(sum(structMesh.*...
            exp(1i*((G(countG,1)-G(countG1,1))*...
            xMesh+(G(countG,2)-G(countG1,2))*yMesh))));
    end
end
%The next loop computes matrix differential operator
%in case of TE mode. The computation
%is carried out for each of earlier defined
%wave vectors.
for countG=1:numG-1
    for countG1=1:numG-1
        for countK=1:length(kx)
            M(countK,countG,countG1)=...
                CN2D_N(countG,countG1)*((kx(countK)+G(countG,1))*...
                (kx(countK)+G(countG1,1))+...
                (ky(countK)+G(countG,2))*(ky(countK)+G(countG1,2)));
        end
    end
end
%The computation of eigen-states is also carried
%out for all wave vectors in the k-path.
for countK=1:length(kx)
    %Taking the matrix differential operator for current
    %wave vector.
    MM(:,:)=M(countK,:,:);
    %Computing the eigen-vectors and eigen-states of
    %the matrix
    [D, V]=eig(MM);
    %Transforming matrix eigen-states to the form of
    %normalized frequency.
    dispe(:,countK)=sqrt(V*ones(length(V),1))*a/2/pi;
end
%Plotting the band structure
%Creating the output field
figure;
%Creating axes for the band structure output.
ax1=axes;
%Setting the option
%drawing without cleaning the plot
hold on;
%Plotting the first 8 bands
for u=1:8
    plot(abs(dispe(u,:)),'r','LineWidth',2);
    %If there are the PBG, mark it with blue rectangle
    if(min(dispe(u+1,:))>max(dispe(u,:)))
        rectangle('Position',[1,max(dispe(u,:)),...
            length(kx)-1,min(dispe(u+1,:))-...
            max(dispe(u,:))],'FaceColor','b',...
            'EdgeColor','b');
    end
end
%Signing Labeling the axes
switch latticeType
    case 0 %Square lattice
        set(ax1,'xtick',...
                [1 precis+1 2*precis+1 3*precis+1]);
        set(ax1,'xticklabel',['G';'X';'M';'G']);
    case 1 %Hexagonal lattice
        set(ax1,'xtick',...
                [1 precis 2*precis-1 3*precis-2]);
        set(ax1,'xticklabel',['G';'M';'K';'G']);
end
ylabel('Frequency \omega a/2\pi c','FontSize',16);
xlabel('Wavevector','FontSize',16);
xlim([1 3*precis-2])
set(ax1,'XGrid','on');