1.软件版本
matlab2013b
2.本算法理论知识
3.部分源码
clc;
clear;
close all;
figure
p = 3*pi/2;
a = 1;
alpha = -15:0.0005:15;
index = 1;
for i=1:length(alpha)
if abs(alpha(i))<0.001
y(i) = p + 1;
else
y(i) = p*sin(a*alpha(i))/(a*alpha(i)) + cos(a*alpha(i));
end
end
plot(alpha,y,'b');hold on;
plot(-15:0.01:15,1,'r-','LineWidth',4);hold on;
plot(-15:0.01:15,-1,'r-','LineWidth',4);hold on;
title('y = p*sin(a*alpha)/(a*alpha) + cos(a*alpha)');
figure_s();
hold off;
%******注意通过以下代码求出图像中的填充区域*******
%四个区间：
%A:2.2540 - 3.1480
%B:4.7206 - 6.2965
%C:7.4290 - 9.4450
%D:10.3040- 13.2755
% startpoint = 2.2540;
% endpoint = 3.1480;
% N = 10000;
% k = 0;
% alpha2 = solvealpha(startpoint,endpoint,N,k);
% alpha2
%******注意通过以上代码求出图像中的填充区域*******
legend('p*sin(a*alpha)/(a*alpha) + cos(a*alpha)','k');
text(pi, -1,'k=pi');
text(2*pi, 1,'k=2pi');
text(3*pi,-1,'k=3pi');
text(4*pi, 1,'k=4pi');
grid on;
figure
%使用四个区间的平方，作为图2的y轴区间
%使用四个区间的平方，作为图2的y轴区间
%使用四个区间的平方，作为图2的y轴区间
%AAAA
startpoint = 2.2540;
endpoint = 3.1480;
N = 10000;
k = -pi:(2*pi)/20:pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
plot(-pi:(2*pi)/80:pi,alpha2(length(alpha2)));hold on;
clear k
clear alpha2
%BBBB
startpoint = 4.7206;
endpoint = 6.2965;
N = 10000;
k = -2*pi:(2*pi)/20:-1*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
k = pi:(2*pi)/20:2*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
plot(-1*pi:(2*pi)/80:1*pi,alpha2(1));hold on;
plot(-2*pi:(2*pi)/80:2*pi,alpha2(length(alpha2)));hold on;
clear k
clear alpha2
%CCCC
startpoint = 7.4290;
endpoint = 9.4450;
N = 10000;
k = -3*pi:(2*pi)/20:-2*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
k = 2*pi:(2*pi)/20:3*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
plot(-2*pi:(2*pi)/80:2*pi,alpha2(1));hold on;
plot(-3*pi:(2*pi)/80:3*pi,alpha2(length(alpha2)));hold on;
clear k
clear alpha2
%DDDD
startpoint = 10.3040;
endpoint = 13.2755;
N = 10000;
k = -4*pi:(2*pi)/20:-3*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
k = 3*pi:(2*pi)/20:4*pi;
for i = 1:length(k)
alpha2(i) = solvealpha(startpoint,endpoint,N,k(i));
alpha2(i) = alpha2(i)*alpha2(i);
end
plot(k,alpha2,'k','LineWidth',4);hold on;
plot(-3*pi:(2*pi)/80:3*pi,alpha2(1));hold on;
plot(-4*pi:(2*pi)/80:4*pi,alpha2(length(alpha2)));hold on;
clear k
clear alpha2
%画虚线
plot(pi, 0:4.7206^2,'r');hold on;
plot(2*pi,0:7.4290^2,'r');hold on;
plot(3*pi,0:10.3040^2,'r');hold on;
plot(4*pi,0:13.2755^2,'r');hold on;
text(-1, 145,'允带');
text(-1, 96,'禁带');
text(-1, 75,'允带');
text(-1, 45,'禁带');
text(-1, 30,'允带');
text(-1, 15,'禁带');
4.仿真结论
5.参考文献
[1]王忆锋, 唐利斌. Kronig-Penney能带模型的MATLAB分析与计算[J]. 激光与红外, 2010, 40(7):5.A16-07
