一个可以算相对介电常数和相对磁导率对角各向异性的FDFD（频域有限差分算法） matlab代码

注意这个代码不是我自己写的，是取自《Electromagnetic and Photonic Simulation for the Beginner Finite-Difference Frequency-Domain in MATLAB》这书的第七章，这本书的PDF在谷歌上很容易下载，并且作者在书中贴出了全书的代码链接。（ps：作者 Raymond C. Rumpf真的太无私了 很了不起的分享精神）

但是这个代码我自己改了一下可以适用于对角各向异性，也就是相对介电常数和相对磁导率可以是对角矩阵，且对角元可以都不等；同时我和comsol计算结果对比了一下，无论是TE还是TM模式的实、虚部都和仿真软件符合得非常好。

仿真对比的comsol文件和.dat数据已上传在蓝奏云，读者可根据合理要求从我这里获得。

% Chapter7_photonicbanddiagram.m% INITIALIZE MATLAB
close all;
clc;
clear all;CRe=dlmread('ComsolRe.dat');
CIm=dlmread('ComsolIm.dat');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DASHBOARD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PHOTONIC CRYSTAL PARAMETERS
a      = 1;%晶格常数
%r      = 0.48*a;
erhole = diag([3-15*1j,5+1j,9+8*1j]);%柱子相对介电常数
urhole = diag([1-1j,4+5*1j,7-5*1j]);%柱子相对磁导率erfill = diag([2-1j,3-1j,4-1j]);%背景相对介电常数
urfill = diag([11-2*1j,10-2*1j,9-2*1j]);%背景相对磁导率len=3*a/5;%正方形柱子边长
xp=[len/2;len/2;-len/2;-len/2];
yp=[len/2;-len/2;-len/2;len/2];
% FDFD PARAMETERS
Nx     = 100;
Ny     = 100;
%NBETA  = 100;%BZ离散点数
NBANDS = 12;
%wnmax  = 0.6;
num_point=10;%BZ一条边上离散点数kxky0=kxky(num_point,a);
NBETA=length(kxky0);
BZpar=0:1/num_point:3;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CALCULATE OPTIMIZED GRID
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% YEE GRID
dx = a/Nx;
dy = a/Ny;% 2X GRID
Nx2 = 2*Nx;             dx2 = dx/2;
Ny2 = 2*Ny;             dy2 = dy/2;% CALCULATE 2X MESHGRID
xa2 = [1-0.5:Nx2-0.5]*dx2;      xa2 = xa2 - mean(xa2);%坐标原点移到中间
ya2 = [1-0.5:Ny2-0.5]*dy2;      ya2 = ya2 - mean(ya2);%坐标原点移到中间
[Y2,X2] = meshgrid(ya2,xa2);%实际坐标 ,先Y2 再X2的顺序可能是为了和后面的数组的维度匹配
%mean函数作用于等差数组上  等价于0.5*(等差数组中最小的数+等差数组中最大的数)
%Y2 X2大小都是Nx2*Ny2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BUILD UNIT CELL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BUILD UNIT CELL
sf=@(xx,yy)shape_f(xx,yy,xp,yp);
[uu,vv]=size(X2);
ER2=zeros(uu,vv);
for aa=1:uufor bb=1:vvER2(aa,bb)=sf(X2(aa,bb),Y2(aa,bb));%逻辑值矩阵end
end%ER2 = (X2.^2 + Y2.^2) <= r^2;%逻辑值矩阵，圆内为1
%ER2 = erfill + (erhole - erfill)*ER2;
ER2xx = erfill(1,1)*(1-ER2) + erhole(1,1)*ER2;%实际相对介电常数分布矩阵
ER2yy = erfill(2,2)*(1-ER2) + erhole(2,2)*ER2;%实际相对介电常数分布矩阵
ER2zz = erfill(3,3)*(1-ER2) + erhole(3,3)*ER2;%实际相对介电常数分布矩阵UR2xx =urfill(1,1)*(1-ER2) + urhole(1,1)*ER2;%实际相对磁导率分布矩阵
UR2yy =urfill(2,2)*(1-ER2) + urhole(2,2)*ER2;%实际相对磁导率分布矩阵
UR2zz =urfill(3,3)*(1-ER2) + urhole(3,3)*ER2;%实际相对磁导率分布矩阵% EXTRACT YEE GRID MATERIAL ARRAYS
ERxx = ER2xx(2:2:Nx2,1:2:Ny2);
ERyy = ER2yy(1:2:Nx2,2:2:Ny2);
ERzz = ER2zz(1:2:Nx2,1:2:Ny2);
URxx = UR2xx(1:2:Nx2,2:2:Ny2);
URyy = UR2yy(2:2:Nx2,1:2:Ny2);
URzz = UR2zz(2:2:Nx2,2:2:Ny2);% SHOW UNIT CELL
figure
title('Unit Cell')
subplot(2,3,1)
imagesc(xa2,ya2,imag(ER2xx.'));%加xa2 ya2是为了标出实际坐标
%imagesc(real(ER2.'));
axis equal tight;
colorbar;
title('er_{xx}');subplot(2,3,2)
imagesc(xa2,ya2,imag(ER2yy.'));%加xa2 ya2是为了标出实际坐标
%imagesc(real(ER2.'));
axis equal tight;
colorbar;
title('er_{yy}');subplot(2,3,3)
imagesc(xa2,ya2,imag(ER2zz.'));%加xa2 ya2是为了标出实际坐标
%imagesc(real(ER2.'));
axis equal tight;
colorbar;
title('er_{zz}');subplot(2,3,4)
imagesc(xa2,ya2,imag(UR2xx.'));%加xa2 ya2是为了标出实际坐标
%imagesc(real(ER2.'));
axis equal tight;
colorbar;
title('ur_{xx}');subplot(2,3,5)
imagesc(xa2,ya2,imag(UR2yy.'));%加xa2 ya2是为了标出实际坐标
%imagesc(real(ER2.'));
axis equal tight;
colorbar;
title('ur_{yy}');subplot(2,3,6)
imagesc(xa2,ya2,imag(UR2zz.'));%加xa2 ya2是为了标出实际坐标
%imagesc(real(ER2.'));
axis equal tight;
colorbar;
title('ur_{zz}');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CALCUALTE LIST OF BLOCH WAVE VECTORS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RECIPROCAL LATTICE VECTORS
% T1 = (2*pi/a) * [ 1 ; 0 ];
% T2 = (2*pi/a) * [ 0 ; 1 ];
% 
% % KEY POINTS OF SYMMETRY
% G = [ 0 ; 0 ];
% X = 0.5*T1;
% M = 0.5*T1 + 0.5*T2;% CHOOSE PATH AROUND IBZ
%KP = [ G X M G ];
KL = { 'M' 'X' '\Gamma' 'M' };% % DETERMINE LENGTH OF IBZ PERIMETER
% NKP  = length(KP);
% LIBZ = 0;
% for m = 1 : NKP-1
%     LIBZ = LIBZ + norm(KP(:,m+1) - KP(:,m));
% end
% 
% % GENERATE LIST OF POINTS AROUND IBZ
% dibz  = LIBZ/NBETA;
% BETA  = KP(:,1);
% KT    = 1;
% NBETA = 1;
% for m = 1 : NKP-1
%     dK      = KP(:,m+1) - KP(:,m);
%     N       = ceil(norm(dK)/dibz);
%     BETA    = [ BETA , KP(:,m) + dK*[1:N]/N ];
%     NBETA   = NBETA + N;
%     KT(m+1) = NBETA;
% end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PERFORM FDFD ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DIAGONALIZE MATERIAL TENSORS
ERxx = diag(sparse(ERxx(:)));%对角稀疏矩阵
ERyy = diag(sparse(ERyy(:)));
ERzz = diag(sparse(ERzz(:)));
URxx = diag(sparse(URxx(:)));
URyy = diag(sparse(URyy(:)));
URzz = diag(sparse(URzz(:)));% INITIALIZE BAND DATA
WNTE = zeros(NBANDS,NBETA);
WNTM = zeros(NBANDS,NBETA);
HermiTE=zeros(1,NBETA);
HermiTM=zeros(1,NBETA);%
% MAIN LOOP -- ITERATE OVER IBZ
%
for nbeta = 1 : length(kxky0)Mark=nbeta/NBETA% Get Next Block Wave Vectorbeta =kxky0(:,nbeta);% Build Derivative MatricesNS  = [Nx Ny];RES = [dx dy];BC  = [1 1];[DEX,DEY,DHX,DHY] = yeeder2d(NS,RES,BC,beta);%求导矩阵的大小都是M*M ,M=Nx*Ny% TM Mode AnalysisA = - DHX/URyy*DEX - DHY/URxx*DEY;B = ERzz;HermiTM(1,nbeta)=ishermitian(inv(B)*A);D = eigs(A,B,NBANDS,0);%0参数表示最接近0的本征值D = sort(D);WNTM(:,nbeta) = D(1:NBANDS);% TE Mode AnalysisA = - DEX/ERyy*DHX - DEY/ERxx*DHY;B = URzz;HermiTE(1,nbeta)=ishermitian(inv(B)*A);D = eigs(A,B,NBANDS,0);D = sort(D);WNTE(:,nbeta) = D(1:NBANDS);end% NORMALIZE THE FREQUENCIES
ReWNTE = a/(2*pi)*real(sqrt(WNTE));%wa/(2*pi*c)
ReWNTM = a/(2*pi)*real(sqrt(WNTM));%wa/(2*pi*c)
ImWNTE = a/(2*pi)*imag(sqrt(WNTE));%wa/(2*pi*c)
ImWNTM = a/(2*pi)*imag(sqrt(WNTM));%wa/(2*pi*c)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DRAW PROFESSIONAL BAND DIAGRAM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
hold on
plot(HermiTE,'b.')
plot(HermiTM,'ro')
hold off
legend('TE','TM')
title('ishermitian')
% PREPARE FIGURE WINDOW% DRAW BAND DATA
figure
subplot(1,2,1)
hold on
plot(BZpar,ReWNTM,'.g');
scatter(CRe(:,1),CRe(:,2),'bo')
hold off
% SET VIEW
xlim([min(BZpar) max(BZpar)]);
%ylim([0 0.6]);
set(gca,'XTick',[0 1 2 3],'XTickLabel',KL);
xlabel('Bloch Wave Vector $\vec{\beta}$',...'Interpreter','LaTex');
ylabel('Frequency $\omega_{\textrm{n}} = a/\lambda_0$',...'Interpreter','LaTex');subplot(1,2,2)
hold on
plot(BZpar,ImWNTM,'.k');
scatter(CIm(:,1),CIm(:,2),'bo')
hold off
title('TM')figure
subplot(1,2,1)
hold on;
plot(BZpar,ReWNTE,'.r');
%scatter(CRe(:,1),CRe(:,2),'bo')
hold off;
% SET VIEW
xlim([min(BZpar) max(BZpar)]);
%ylim([0 0.6]);
set(gca,'XTick',[0 1 2 3],'XTickLabel',KL);
xlabel('Bloch Wave Vector $\vec{\beta}$',...'Interpreter','LaTex');
ylabel('Frequency $\omega_{\textrm{n}} = a/\lambda_0$',...'Interpreter','LaTex');subplot(1,2,2)
hold on
plot(BZpar,ImWNTE,'.k');
%scatter(CIm(:,1),CIm(:,2),'bo')
hold off;
title('TE')function [op]=kxky(num_point,a)
x0=0:1/num_point:3;
op=zeros(2,length(x0));
for jj=1:length(x0)x=x0(jj);if x>=0 && x<1op(:,jj)=[pi/a; (1-x)*(pi/a)];elseif x>=1 && x<2op(:,jj)=[(2-x)*(pi/a); 0];elseif x>=2 && x <=3op(:,jj)=[(x-2)*(pi/a); (x-2)*(pi/a)];endendendfunction [out]=shape_f(xx,yy,xp,yp)
%xx、yy是输入的要判断的点，
%xp 、yp是多边形的顶点(顺时针或逆时针，首尾相连)
%xp、yp是向量
%如果在多边形内或边上，out为1，否则为0
out=double(inpolygon(xx,yy,xp,yp));endfunction [DEX,DEY,DHX,DHY] = yeeder2d(NS,RES,BC,kinc)
% YEEDER2D      Derivative Matrices on a 2D Yee Grid
%
% [DEX,DEY,DHX,DHY] = yeeder2d(NS,RES,BC,kinc);
%
% INPUT ARGUMENTS
% =================
% NS    [Nx Ny] Grid Size
% RES   [dx dy] Grid Resolution
% BC    [xbc ybc] Boundary Conditions
%         0: Dirichlet boundary conditions
%         1: Periodic boundary conditions
% kinc  [kx ky] Incident Wave Vector
%       This argument is only needed for PBCs.
%
% Note: For normalized grids use k0*RES and kinc/k0
%
% OUTPUT ARGUMENTS
% =================
% DEX   Derivative Matrix wrt x for Electric Fields
% DEY   Derivative Matrix wrt to y for Electric Fields
% DHX   Derivative Matrix wrt to x for Magnetic Fields
% DHY   Derivative Matrix wrt to y for Magnetic Fields%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% HANDLE INPUT ARGUMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EXTRACT GRID PARAMETERS
Nx = NS(1);     dx = RES(1);
Ny = NS(2);     dy = RES(2);% DEFAULT KINC
if ~exist('kinc')kinc = [0 0];
end% DETERMINE MATRIX SIZE 确定矩阵大小
M = Nx*Ny;% ZERO MATRIX
Z = sparse(M,M);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BUILD DEX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HANDLE IF SIZE IS 1 CELL
if Nx==1DEX = -1i*kinc(1)*speye(M,M); %为什么？% HANDLE ALL OTHER CASES
else% Center Diagonald0 = -ones(M,1);%这里必须为一列数组% Upper Diagonald1 = ones(M,1);%这里必须为一列数组d1(Nx+1:Nx:M) = 0;% Build Derivative Matrix with Dirichlet BCsDEX = spdiags([d0 d1]/dx,[0 1],Z);%把数值填到Z里 其他部分不变% Incorporate Periodic Boundary Conditionsif BC(1)==1d1 = zeros(M,1);%这里必须为一列数组d1(1:Nx:M) = exp(-1i*kinc(1)*Nx*dx)/dx;DEX = spdiags(d1,1-Nx,DEX);endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BUILD DEY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HANDLE IF SIZE IS 1 CELL
if Ny==1DEY = -1i*kinc(2)*speye(M,M);% HANDLE ALL OTHER CASES
else% Center Diagonald0 = -ones(M,1);% Upper Diagonald1 = ones(M,1);% Build Derivative Matrix with Dirichlet BCsDEY = spdiags([d0 d1]/dy,[0 Nx],Z);% Incorporate Periodic Boundary Conditionsif BC(2)==1d1 = (exp(-1i*kinc(2)*Ny*dy)/dy)*ones(M,1);DEY = spdiags(d1,Nx-M,DEY);endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BUILD DHX AND DHY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DHX = -DEX';
DHY = -DEY';
end

不得不说Rumpf教授的FDFD代码是我目前见过的算得和comsol软件吻合得最好的代码！