matlab仿真二维光子晶体最简程序

本程序为初学者使用,只考虑MT方向

下面的程序为matlab代码

只考虑MT方向

%This is a simple demo for Photonic Crystals simulation
%This demo is for TE wave only, so only h wave is considered.
%for TM direction only,10 points is considered.
%---------------------------------------M
%|                                         /    |
%|                                   /          |
%|                             /                |
%|                      --------------------|X
%|                     T                       |
%|                                              |
%|                                              |
%---------------------------------------
%equation :sum_{G‘,k}(K+G)(K+G‘)f(G-G‘)hz(k+G‘)=(omega/c)^2*hz(k+G)
%G‘ can considerd as the index of column, and G as index of rows
%[(K+G1)(K+G1)f(G1-G1)   (K+G1)(K+G2)f(G1-G2)  ][hz(G1)]=(omega/c)^2[hz(G1)]
%[(K+G2)(K+G1)f(G2-G1)   (K+G2)(K+G2)f(G2-G2)  ][hz(G2)]            [hz(G2)]
%or:   THETA_TE*Hz=(omega/c)^2*Hz
%by Gao Haikuo
%date:20170411

clear; clc; epssys=1.0e-6; %设定一个最小量,避免系统截断误差或除0错误

%this is the lattice vector and the reciprocal lattice vector
a=1; a1=a*[1 0]; a2=a*[0 1];
b1=2*pi/a*[1 0];b2=2*pi/a*[0 1];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%定义晶格的参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsa = 1; %介质柱的介电常数
epsb = 13; %背景的介电常数
Pf = 0.7; %Pf = Ac/Au 填充率,可根据需要自行设定
Au =a^2; %二维格子原胞面积
Rc = (Pf *Au/pi)^(1/2); %介质柱截面半径
Ac = pi*(Rc)^2; %介质柱横截面积

%construct the G  list
NrSquare = 10;
NG =(2*NrSquare+1)^2;  % NG is the number of the G value
G = zeros(NG,2);
i = 1;
for l = -NrSquare:NrSquare
    for m = -NrSquare:NrSquare
        G(i,:)=l*b1+m*b2;
        i = i+1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成k空间中的f(Gi-Gj)的值,i,j 从1到NG。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=zeros(NG,NG);
for i=1:NG
    for j=1:NG
        Gij=norm(G(i,:)-G(j,:));
        if (Gij < epssys)
            f(i,j)=(1/epsa)*Pf+(1/epsb)*(1-Pf);
        else
            f(i,j)=(1/epsa-1/epsb)*Pf*2*besselj(1,Gij*Rc)/(Gij*Rc);
        end;
    end;
end;
T=(2*pi/a)*[epssys 0];
M=(2*pi/a)*[1/2 1/2]; %????????
X=(2*pi/a)*[1/2 0]; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%对于简约布里渊区边界上的每个k,求解其特征频率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
THETA_TE=zeros(NG,NG); %待解的TE波矩阵
Nkpoints=10; %每个方向上取的点数,
stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长

TX_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵
XM_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵
MT_TE_eig = zeros(Nkpoints,NG); %沿MT方向的TE波的待解的特征频率矩阵

for n=1:Nkpoints %scan the 10 points along the TM direction
    fprintf([‘\n k-point:‘,int2str(n),‘of‘,int2str(Nkpoints),‘.\n‘]);
    MT_step = stepsize(n)*(T-M)+M;  % get the k
    %先求非对角线上的元素
    for i=1:(NG-1)   % G
        for j=(i+1):NG % G‘
            kGi = TX_step+G(i,:); %k+G
            kGj = TX_step+G(j,:); %K+G‘
            THETA_TE(i,j)=f(i,j)*dot(kGi,kGj); %(K+G)(K+G‘)f(G-G‘)
            THETA_TE(j,i)=conj(THETA_TE(i,j));
        end
    end
    %再求对角线上的元素
    for i=1:NG
        kGi = TX_step+G(i,:);
        THETA_TE(i,i)=f(i,i)*norm(kGi)*norm(kGi);
    end
    MT_TE_eig(n,:)=sort(sqrt(eig(THETA_TE))).‘;
end

%draw
kaxis = 0;
TXaxis = kaxis:norm(T-X)/(Nkpoints-1):(kaxis+norm(T-X));
kaxis = kaxis + norm(T-X);
XMaxis = kaxis:norm(X-M)/(Nkpoints-1):(kaxis+norm(X-M));
kaxis = kaxis + norm(X-M);
MTaxis = kaxis:norm(M-T)/(Nkpoints-1):(kaxis+norm(M-T));
kaxis = kaxis + norm(M-T); 

Ntraject = 3;
figure (1)
hold on;
Nk=Nkpoints;
for k=1:NG
    for i=1:Nkpoints
        EigFreq_TE(i+0*Nk) = TX_TE_eig(i,k)/(2*pi/a);
        EigFreq_TE(i+1*Nk) = XM_TE_eig(i,k)/(2*pi/a);
        EigFreq_TE(i+2*Nk) = MT_TE_eig(i,k)/(2*pi/a);
    end
  plot(TXaxis(1:Nk),EigFreq_TE(1+0*Nk:1*Nk),‘r‘,...
XMaxis(1:Nk),EigFreq_TE(1+1*Nk:2*Nk),‘r‘,...
MTaxis(1:Nk),EigFreq_TE(1+2*Nk:3*Nk),‘r‘);
end
时间: 2024-12-27 16:19:36

matlab仿真二维光子晶体最简程序的相关文章

二维光子晶体带隙绘制程序_平面波展开法

1.主程序 %This is a simple demo for Photonic Crystals simulation %10 points is considered. %by Gao Haikuo %date:20170411 clear; clc; global NG G f Nkpoints eigenValue modeset kCorner global epsa epsb epssys a b1 b2 epssys=1.0e-6; %设定一个最小量,避免系统截断误差或除0错误

matlab画二维直方图以及双y轴坐标如何修改另一边y轴的颜色

1.首先讲一下如何用hist画二维直方图 1 x=[-568179 -766698 -935586 -826865 -393971 -771826 -1529945 -1910695 -1694740 -926367 -306998 -844840 -1828334 -2062815 -2297296 -1498824 -411346 -827922 -1826636 -1844777 -1862918 -1881060 -746534 -100479 -845832 -1832756 -194

笔记12:简易的二维码生成解析程序

首先得引用一个文件.直接看代码吧! 1 using System; 2 using System.Collections.Generic; 3 using System.ComponentModel; 4 using System.Drawing; 5 using System.Data; 6 using System.Text; 7 using System.Windows.Forms; 8 using ThoughtWorks.QRCode.Codec; 9 using ThoughtWor

java 获取数组(二维数组)长度实例程序

我们可能知道 js有个length函数,java也有啊length函数 例 如果数组是data[],则data.length 代码如下 复制代码 byte[] phone =new byte[81]; //建立一个byte类型的数组,长度为81 phone[i]!=0中phone[i]! //数组的第i的位置不等0 如: 代码如下 复制代码 byte[] phone =new byte[81]; //建立一个byte类型的数组,长度为81 phone[1]!=0中phone[1]! //数组第二

Android zxing 解析二维码,生成二维码极简demo

zxing 官方的代码很多,看起来很费劲,此demo只抽取了有用的部分,实现了相机预览解码,解析本地二维码,生成二维码三个功能. 简化后的结构如下: 废话少说直接上代码: BaseDecodeHandler: package com.song.zxing.decode; import android.graphics.Bitmap; import android.os.Bundle; import com.google.zxing.BarcodeFormat; import com.google

使用Graham扫描法求二维凸包的一个程序

1 #include "includeall.h" 2 #include "Link.class.h" 3 4 int RightOrLeft(float x1,float y1,float x2,float y2,float x3,float y3)//判断第三个点在前两个点连成的直线的哪个位置,-1 左边,0,直线上,1 右边 5 { 6 float X=(y3-y1)*(x2-x1)/(y2-y1)+x1; 7 if(X<x3) 8 { 9 return

matlab 画二维图与三维图

二维图 ezplot('sin(x)');%默认范围 ezplot('sin(x)',[-4 4]);%自己设定范围 三维图 ezmesh('x*x+y*y');%默认范围

matlab的二维卷积操作(转)

MATLAB的conv2函数实现步骤(conv2(A,B)): 其中,矩阵A和B的尺寸分别为ma*na即mb*nb ① 对矩阵A补零,第一行之前和最后一行之后都补mb-1行,第一列之前和最后一列之后都补nb-1列(注意conv2不支持其他的边界补充选项,函数内部对输入总是补零): ② 将卷积核绕其中心旋转180度: ③ 滑动旋转后的卷积核,将卷积核的中心位于图像矩阵的每一个元素,并求乘积和(即将旋转后的卷积核在A上进行滑动,然后对应位置相乘,最后相加):下面分别是shape=full, same

matlab的二维卷积操作

MATLAB的conv2函数实现步骤(conv2(A,B)): 其中,矩阵A和B的尺寸分别为ma*na即mb*nb ① 对矩阵A补零,第一行之前和最后一行之后都补mb-1行,第一列之前和最后一列之后都补nb-1列(注意conv2不支持其他的边界补充选项,函数内部对输入总是补零): ② 将卷积核绕其中心旋转180度: ③ 滑动旋转后的卷积核,将卷积核的中心位于图像矩阵的每一个元素,并求乘积和(即将旋转后的卷积核在A上进行滑动,然后对应位置相乘,最后相加):下面分别是shape=full, same