二维光子晶体带隙绘制程序_平面波展开法(最终版)
1.主程序
%This is a simple demo for Photonic Crystals simulation %10 points is considered. %by Gao Haikuo %date:20170411clear; clc; global NG G f Nkpoints eigenValue modeset kCorner global epsa epsb epssys a b1 b2 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]; Nkpoints=10; %每个方向上取的点数, modeset=2;% 1:'TE' 2 'TM' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %定义晶格的参数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% epsa = 8.9; %inner epsb = 1; %outer Pf = 0.1257; %Pf = Ac/Au 填充率,可根据需要自行设定Au =a^2; %二维格子原胞面积Rc = (Pf *Au/pi)^(1/2); %介质柱截面半径Ac = pi*(Rc)^2; %介质柱横截面积kCorner=(2*pi/a)*[epssys 0;1/2 0;1/2 1/2]; %T X M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%get gap[G,f]=getGAndf(Pf,Rc);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%get gapeigenValue=getFrequency(kCorner);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%get gapgap=getGap();%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%绘draw banddrawBand(gap);
本程序又可以分为以下几个步骤:
- 求解G和f
- 求解本征频率
- 求解光子带隙
- 绘图
2.求解G和f
function [G,f]=getGAndf(Pf,Rc) global NG G f Nkpoints eigenValue kCorner modeset a global epsa epsb epssys b1 b2 NrSquare = 10; NG =(2*NrSquare+1)^2; % NG is the number of the G value G = zeros(NG,2); i = 1; for l = -NrSquare:NrSquarefor m = -NrSquare:NrSquareG(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;
3.求解本征频率
function eigenValue=getFrequency(kCorner) global NG Nkpoints modeset n_kCorner=size(kCorner,1); kCorner_ex=[kCorner;kCorner(1,:)]; eigenValue=zeros(Nkpoints,NG,n_kCorner,2); for i_mode=modesetfor i_kCorner=1:n_kCornereigenValue(:,:,i_kCorner,i_mode)=...getEigValue(kCorner_ex(i_kCorner,:),...kCorner_ex(i_kCorner+1,:),i_mode);end end
该程序调用了求解一个布里渊边界本征频率的子程序
function [eigValue]=getEigValue(k_begin,k_end,mode) %mode :1 for TE 2 for TM global NG G f Nkpoints THETA=zeros(NG,NG); %待解的TE波矩阵 stepsize=0:1/(Nkpoints-1):1; %每个方向上的步长 TX_TE_eig = zeros(Nkpoints,NG); for n=1:Nkpoints %scan the 10 points along the directionfprintf(['\n k-point:',int2str(n),'of',int2str(Nkpoints),'.\n']);step = stepsize(n)*(k_end-k_begin)+k_begin; % get the kfor i=1:(NG-1) % Gfor j=(i+1):NG % G' kGi = step+G(i,:); %k+GkGj = step+G(j,:); %K+G' switch modecase 1 %TE modeTHETA(i,j)=f(i,j)*dot(kGi,kGj); %(K+G)(K+G')f(G-G')case 2 %TH modeTHETA(i,j)=f(i,j)*norm(kGi)*norm(kGj);endTHETA(j,i)=conj(THETA(i,j)); endend for i=1:NGkGi = step+G(i,:);THETA(i,i)=f(i,i)*norm(kGi)*norm(kGi); endeigValue(n,:)=sort(sqrt(eig(THETA))).'; end
4.求解光子晶体带隙
function gap=getGap() global eigenValue kCorner modeset NG a gap=[]; band=[]; n_kCorner=size(kCorner,1); kCorner_ex=[kCorner;kCorner(1,:)]; for mode=modeset for i_kCorner=1:n_kCornerfor k=1:NGsubLine=eigenValue(:,k,i_kCorner,mode)*a/(2*pi);L=min(subLine);H=max(subLine);band=mergeBand(band,L,H);endend end if size(band,1)>=2gap=[band(1:end-1,2),band(2:end,1)]; end
里面调用了区间合并的函数mergeBand
function band=mergeBand(band,L,H) flag=1; i_L=0;%最小值所在行数 flag_L_in=0;%是否在之间 i_H=0; flag_H_in=0; i=1; flag_H_scan=1; if isempty(band)band=[L,H];return; end while flagif L<band(i,1)i_L=i;flag=0;flag_L_in=0;elseif L<=band(i,2)i_L=i;flag=0;flag_L_in=1;elsei=i+1;if i>length(band(:,1))flag=0;flag_H_scan=0;band=[band;L,H];endend end flag=1; if flag_H_scanwhile flagif H<band(i,1)i_H=i;flag=0;flag_H_in=0;elseif H<=band(i,2)i_H=i;flag=0;flag_H_in=1;elsei=i+1;if i>length(band(:,1))flag=0;i_H=i;endendend%mergeif i_L==i_H if L<band(i_L,1) && H>=band(i_H,1)band(i_L,1)=L;elseif H<band(i_L,1) band=[band(1:i_L-1,:);[L,H];band(i_L:end,:)];end elseif i_H>length(band(:,1))band(i_L,1)=min([band(i_L,1),L]);band(i_L,2)=H;if i_L+1<=length(band(:,1))band(i_L+1:end,:)=[];endelseif H>=band(i_H,1) band(i_L,1)=min([band(i_L,1),L]);band(i_L,2)=max([band(i_H,2),H]);band(i_L+1:i_H,:)=[];elseband(i_L,1)=min([band(i_L,1),L]);band(i_L,2)=H; if i_L+1<=i_H-1band(i_L+1:i_H-1,:)=[];endend end end
5.最后是绘制图形的程序
function drawBand(gap) global NG G f Nkpoints eigenValue kCorner modeset a n_kCorner=size(kCorner,1); figure (1) hold on; colorSet=['r','b']; for mode=modesetfor i_kCorner=1:n_kCornerfor k=1:NGh=plot((i_kCorner-1)*(Nkpoints-1) : i_kCorner*(Nkpoints-1),...eigenValue(:,k,i_kCorner,mode)*a/(2*pi),colorSet(mode));set(h, 'linesmoothing', 'on');endend end grid on; xlabel('K-Space'); yLabel('Frequency(\omegaa/2\piC)'); xmax=n_kCorner*(Nkpoints-1); axis([0 xmax 0 0.8]); set(gca,'XTick',0:(Nkpoints-1):xmax); xtixlabel = strvcat('T','X','M','T'); set(gca,'XTickLabel',xtixlabel); %draw gap for i=1:size(gap,1)fill([0, xmax,xmax,0],[gap(i,1),gap(i,1),gap(i,2),gap(i,2)],'r'); end hold off;
转载于:https://www.cnblogs.com/Iknowyou/p/6723749.html
二维光子晶体带隙绘制程序_平面波展开法(最终版)相关推荐
- matlab晶体能带,matlab平面波展开法的二维光子晶体能带研究+程序
摘 要 :二维光子晶体可以作为对光子传输控制的新型材料.本文主要通过平面波展开法对二维光子晶体进行数值计算及其性质分析.首先我们介绍了二维光子晶体的基础概念.结构.介电性能等特性.然后基于麦克斯韦方 ...
- 如何创建二维数组 微信小程序_微信小程序遍历二维数组
在微信小程序中遍历二维数组,代码如下 data 中二维数组结构如下 data: { familys:[ { familyName:'贾家', users: [ {name:'贾宝玉'}, {name: ...
- 『小程序开发』关于微信小程序扫普通链接二维码打开小程序的具体配置流程...
前言: 对于扫普通链接二维码打开小程序的功能详解,官方api已经可以说是接近手把手的教学,咱们这里不做累述,直接上图走起...官方接入指南 功能介绍 扫二维码登录小程序...^_^ 限制 1.对于普通 ...
- 绘制恒定频率下的二维光子晶体光子带图谱
绘制恒定频率下的二维光子晶体光子带图谱 在光子晶体研究领域中,光子带图谱是非常重要的工具.它可以展示出不同频率下光的传播状态,并且为设计光子晶体的光子学性质提供依据.本篇文章将介绍如何使用Matlab ...
- matlab三维选取二维,基于Matlab绘制二维和三维图形以及其他图形控制函数的使用方法...
Matlab绘图 强大的绘图功能是Matlab的特点之一,Matlab提供了一系列的绘图函数,用户不需要过多的考虑绘图的细节,只需要给出一些基本参数就能得到所需图形,这类函数称为高层绘图函数.此外,M ...
- 二维光子晶体禁带的遗传优化算法实现
二维光子晶体禁带的遗传优化算法MATLAB源代码 光子晶体中因周期性结构而存在的频率禁带称为光子禁带,光子禁带的存在是光子晶体具有广泛应用前景的重要原因. 禁带越大,可控光的频带也越宽,因此如何设计合 ...
- uni-app微信小程序扫普通二维码分享小程序
这里需要扫普通二维码分享的话就需要先产生二维码了 文档:https://github.com/yingye/weapp-qrcode 1.绘制二维码 我这里使用的是资源是weapp.qrcode.es ...
- matlab建成二维数组,matlab绘制二维数组
hist 累计图 rose 极座标累计图 stairs 阶梯图 stem 针状图 fill 实心图 feather 羽毛图 compass 罗盘图 quiver 向量场图 Matlab 如何画出一个二 ...
- matlab二维谐振子,基于有限差分法求解的二维谐振子的MATLAB程序如下。哪位大神能帮我做个注明啊,完全看不懂啊,,急...
基于有限差分法求解的二维谐振子的MATLAB程序如下.哪位大神能帮我做个注明啊,完全看不懂啊,,急0 ____丿呆呆丶2017.04.15浏览20次分享举报 tic clc clear L=20; W ...
最新文章
- 关于.NET技术体系的思维导图
- Android中JSON解析细解及实例
- 国内外软件开发上的差距与分析
- apache目录 vscode_VsCode搭建Java开发环境(Spring Boot项目创建、运行、调试)
- python打代码运行图形_利用aardio给python编写图形界面
- QT学习小结之信号与槽
- C#WinForm的TextBox 按TAB键让光标按照指定顺序走
- 让你做个《五子棋》怎么存储棋盘上的棋子信息?
- cad立体图怎么旋转看图_CAD趣事之对CAD图纸进行旋转,360°无死角查看的方法-dwg文件查看器...
- Visio 2003 Professional
- React Native集成阿里云推送----广播推送
- 调频连续波雷达(FMCW)测距/测速原理
- 【数据结构与算法】试卷 1(含答案)
- 小型的代码管理仓库Gitea安装指南
- B站下载姿势合集——亲测
- JAVA集合05_Collection.toMap()应用、三个重载方法、解决重复key问题
- centos7搭建hexo博客步骤
- 模糊数学 5、模糊综合评判
- 边缘设备、系统及计算杂谈(12)——k8s学习之二
- Linux访问Windows7共享文件夹