onebeam3d.m

理论模拟计算干涉法生成二维及三维光子晶体的结构和各个参数的调节对结果的影响

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 = [20, 20, 20] * 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 + 1
    if k < iBeamNum + 1
        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,:);
%-----------------------------------------------------------------------------------------
    else
        vE(k,:) = vSourcePolarization;
        vK(k,:) = vSourceK;
    end
end

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

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

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

    for j = 0 : n
        y = 0.5 * j;
        z = 4;
        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 + 1
          for l = 1 : iBeamNum + 1
             dSum = dSum + abs (vE(k,:) * vE(l,:)') * cos ((vK(k,:) - vK(l,:)) * vP');
          end
       end
%        if dSum > 2
           mInten (m + 1, j + 1) = dSum;
%        else
%            mInten (m + 1, j + 1) = 0;
%        end
% method 2
       
     end
end

% Draw the pattern

pcolor (mInten);
axis equal;