develop.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;
lambda= 5.32;
% 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 = [1, 0, 0];
vSourceK = [0, 0, 1]*(2*pi/lambda);
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
Sourcepower=[3,0,0];
sk=[1,0,0];
pk=[0,1,0];
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 
    if k < iBeamNum 
        beam(k,:) = [sin(vBeamAngleZ(k)) * cos(vBeamAngleX(k)),...
                  sin(vBeamAngleZ(k)) * sin(vBeamAngleX(k)),...
                  cos(vBeamAngleZ(k))];
%beam1=[0,1,0];
%beam2=[1*sin(vIncidentAngle)*cos(pi/6),1*sin(vIncidentAngle)*sin(pi/6),1*sin(vIncidentAngle)];
%beam3=[-1*sin(vIncidentAngle)*cos(pi/6),-1*sin(vIncidentAngle)*sin(pi/6),1*sin(vIncidentAngle)];
beam4=[1.5,1.5,0];
n = 100;
for m = 0 : n
    x = 1.5 * m;
    for j = 0 : n
        y = 1.5 * j;
        for p = 0 : 80 
            z = 1.5 * p;
        vP=[y,x,z];
        dSum = 0;
       for k = 1 : iBeamNum
          for l = 1 : iBeamNum 
              dSum=dSum+abs(beam(k,:)*beam(l,:)') * cos ((beam(k,:) - beam(l,:)) * vP');
              
              
          end
      end
        mInten (m + 1, j + 1, p + 1) = dSum;
    end
end
end
% % Draw the pattern
% 
% pcolor (mInten);
% axis equal;

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