beam2d_no_polarization.m

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

clear;
% % The function can be used to calculate the pattern of interference

% % Just for the case of one beam propagates through a prism which can devides and combines the beam

% the parameter
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

vSourceK = [0, 0, 1];                                             %the wave vector

%vBeamAngleX = [0,2/3,4/3] * pi;% angle from x-axis, the projection of beams on X-Y plane

%--------------- the case of the angles between each beam is equal -------------------------
vBeamAngleX = zeros(1,iBeamNum);
for n = 1:iBeamNum-1
    vBeamAngleX(n+1) = vBeamAngleX(n) + 2*pi/iBeamNum;
end
%----------------------------------------------------

PlaneAngle = 38.9;                                                % angle between the plane and the source-beam
vPlaneAngle = (PlaneAngle * pi / 180).*ones(1,iBeamNum);

BeamPower = 1;                                    % the amplitude
vBeamPower = BeamPower.*ones(1,iBeamNum);                                           

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

%vBeamAngleZ = 40.3/180*pi.*[1,1,1];      %  the change of the angle between beam and Z axis

for k = 1 : iBeamNum
%------------ the wave vector of each beam --------------------
    vK(k,:) = [sin(vBeamAngleZ(k)) * cos(vBeamAngleX(k)),...
               sin(vBeamAngleZ(k)) * sin(vBeamAngleX(k)),...
               cos(vBeamAngleZ(k))];
end

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

n = 500;
% 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 = 0;
        vP=[y,x,z];
        

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

% Draw the pattern

pcolor (mInten);
%figure
%imshow(mInten,[]);
axis equal;