2d-fdtd.m

二维光子晶体禁带的MATLAB源码,光子晶体中因周期性结构而存在的频率禁带称为光子禁带

clear;
% % The function can be used to calculate the pattern of interference
% % Just for one beam case

% the param
i = sqrt (-1);
dNglass = 1.58;                                   % the incident index of glass;

% input the configure of prism ------------------------------------------------------------

iBeamNum = 3;                                         % the number of the prism's plane
dAlpha = 0 * pi / 180;                               % the polarization of the source beam
vSourcePolarization = [cos(dAlpha), sin(dAlpha), 0];
vSourceK = [0, 0, 1];
vBeamAngleX = [0, 120, 240] * pi / 180;          % angle from x-axis
vPlaneAngle = [50.9, 50.9, 50.9] * pi / 180;    % angle between the plane and the source-beam
vBeamPower = [1, 1, 1];                            % the amplitude

%------------------------------------------------------------------------------------------
vRefractionAngle = asin (sin (pi / 2 - vPlaneAngle) / dNglass); % the refraction angle
vIncidentAngle = pi / 2 - vPlaneAngle; 
vBeamAngleZ = vIncidentAngle - vRefractionAngle; % angle from z-axis

for k = 1 : iBeamNum
        vK(k,:) = [sin(vBeamAngleZ(k)) * cos(vBeamAngleX(k)),...
                  sin(vBeamAngleZ(k)) * sin(vBeamAngleX(k)),...
                  cos(vBeamAngleZ(k))];
   
% considered the polarization of the beams (For Method 2)-----------------------------------
% In the air
        vSair(k,:) = cross (vK(k,:), vSourceK);
        vSair(k,:) = vSair(k,:) / sqrt (vSair(k,:) * vSair(k,:)');
        vPair(k,:) = cross (vSair(k,:), vSourceK);
        vPair(k,:) = vPair(k,:) / sqrt (vPair(k,:) * vPair(k,:)');
        dPowerP(k) = vPair(k,:) * vSourcePolarization' * vBeamPower(k);
        dPowerS(k) = vSair(k,:) * vSourcePolarization' * vBeamPower(k);
   
% In the glass
        vSglass(k,:) = vSair(k,:);
        vPglass(k,:) = cross (vSglass(k,:), vK(k,:));
        vPglass(k,:) = vPglass(k,:) / sqrt (vPglass(k,:) * vPglass(k,:)');
        dPowerP(k) = dPowerP(k) * 2 * cos (vIncidentAngle(k)) / (dNglass * cos(vIncidentAngle(k)) + cos(vRefractionAngle(k)));
        dPowerS(k) = dPowerS(k) * 2 * cos (vIncidentAngle(k)) / (cos(vIncidentAngle(k)) + dNglass * cos(vRefractionAngle(k)));

% the polarization of the refraction beams
        vE(k,:) = dPowerP(k) * vPglass(k,:) + dPowerS(k) * vSglass(k,:);
%-----------------------------------------------------------------------------------------
end

% Calculate the pattern -------------------------------------------------------------------

n = 100;
% mInten(n+1,n+1)=zeros(n+1,n+1);

for m = 0 : n
    x = 1.5 * m;

    for j = 0 : n
        y = 1.5 * j;
        for p = 0 : 10 
            z = 1.5 * p;
        vP=[y,x,z];
        

%% calculate the intensity -----------------------------------------------------------------
% % method 1
%       vSum = 0;
%       for k = 1 : iBeamNum
%           vSum = vSum + vBeamPower(k) * exp (i * vK(k,:) * vP');
%       end
% 
%         if (vSum * vSum' > 8)
%           mInten (m + 1, j + 1) = vSum * vSum';
%         end
% % method 1 

% method 2
       dSum = 0;
       for k = 1 : iBeamNum
          for l = 1 : iBeamNum
             dSum = dSum + abs (vE(k,:) * vE(l,:)') * cos ((vK(k,:) - vK(l,:)) * vP');
          end
       end
%        if dSum > 2
           mInten (m + 1, j + 1, p + 1) = dSum;
%        else
%            mInten (m + 1, j + 1) = 0;
%        end
% method 2
       
     end
end
end
% % Draw the pattern
% 
% pcolor (mInten);
% axis equal;

% draw the 3D graphic
mTemp = smooth3(mInten);
hiso = patch(isosurface(mTemp, 1),...
    'FaceColor','blue','EdgeColor','none');
hcap = patch(isocaps(mTemp, 1),...
    'FaceColor','blue','EdgeColor','none');
isonormals(mTemp,hiso)
view(3); axis vis3d tight
camlight left; lighting phong