ПЗ / 3 / Computing

clear;
%данные
me=9.10*(10^-31);
h=1.05*(10^-34);
l=0.5*(10^-9);
W1=6.409*10^(-19);
W2=0;
W3=1.282*10^(-18);

%находим меньший из барьеров и вычисляем коэффициент b
b=min(W3,W1);
Wmin=min(W3,W1);

%вычисление пси-функций для (не)сим. ямы
%данные для несим. ямы
z1=[0.21 0.78]; %корни графического решения уравнения для несим. ямы
E1=zeros();
delt1=zeros();
A1=zeros();
A11=zeros;
A21=zeros();
B1=zeros();
C1=zeros();
s1=zeros();
k11=zeros();
k31=zeros();
psi1=zeros();
psi21=zeros();

%данные для сим. ямы
z2=[0.19 0.7]; %корни графического решения уравнения для сим. ямы
E2=zeros();
delt2=zeros();
A2=zeros();
A12=zeros;
A22=zeros();
B2=zeros();
C2=zeros();
s2=zeros();
k12=zeros();
k32=zeros();
psi2=zeros();
psi22=zeros();

x=-0.5*l:0.001*(10^-9):1.5*l;

%расчёт параметров для несим. ямы
for i=1:size(z1,2)
%расчёт энергии E(n)
E1(i)=z1(i)*Wmin;
%вычисление s
s1(i)=(sqrt(2*me*(E1(i)-W2)))/(h);
%вычисление k1
k11(i)=(sqrt(2*me*(W1-E1(i))))/(h);
%вычисление k3
k31(i)=(sqrt(2*me*(W3-E1(i))))/(h);
%вычисление δ(n)
delt1(i)=atan(sqrt((E1(i))/(Wmin-E1(i))));
%вычисление константы A
A11(i)=4*k11(i)*k31(i)*s1(i)*(sin(delt1(i)))^2;
A21(i)=2*k31(i)*s1(i)*(sin(delt1(i)))^2-k11(i)*k31(i)*sin(2*l*s1(i)+2*delt1(i))+k11(i)*k31(i)*sin(2*delt1(i))+2*k11(i)*k31(i)*l*s1(i)+2*k11(i)*s1(i)*(sin(l*s1(i)+delt1(i)))^2;
A1(i)=sqrt((A11(i))/(A21(i)));
%вычисление константы B
B1(i)=A1(i)/sin(delt1(i));
%вычисление константы C
C1(i)=(A1(i)*sin(s1(i)*l+delt1(i)))/(sin(delt1(i))*exp(-k31(i)*l));
end

%расчёт параметров для сим. ямы
for i=1:size(z2,2)
%расчёт энергии E(n)
E2(i)=z2(i)*Wmin;
%вычисление s
s2(i)=(sqrt(2*me*(E2(i)-W2)))/(h);
%вычисление k1
k12(i)=(sqrt(2*me*(Wmin-E2(i))))/(h);
%вычисление k3
k32(i)=(sqrt(2*me*(Wmin-E2(i))))/(h);
%вычисление δ(n)
delt2(i)=atan(sqrt((E2(i))/(Wmin-E2(i))));
%вычисление константы A
A12(i)=4*k12(i)*k32(i)*s2(i)*(sin(delt2(i)))^2;
A22(i)=2*k32(i)*s2(i)*(sin(delt2(i)))^2-k12(i)*k32(i)*sin(2*l*s2(i)+2*delt2(i))+k12(i)*k32(i)*sin(2*delt2(i))+2*k12(i)*k32(i)*l*s2(i)+2*k12(i)*s2(i)*(sin(l*s2(i)+delt2(i)))^2;
A2(i)=sqrt((A12(i))/(A22(i)));
%вычисление константы B
B2(i)=A2(i)/sin(delt2(i));
%вычисление константы C
C2(i)=(A2(i)*sin(s2(i)*l+delt2(i)))/(sin(delt2(i))*exp(-k32(i)*l));
end

%расчёт значений пси-функции для несим. ямы
for j=1:size(z1,2)
for i=1:size(x,2)
if x(i)<=0
psi1(i,j)=A1(j)*exp(x(i)*k11(j));
psi21(i,j)=(psi1(i,j))^2;
elseif (0<x(i))&&(x(i)<=l)
psi1(i,j)=B1(j)*sin(x(i)*s1(j)+delt1(j));
psi21(i,j)=(psi1(i,j))^2;
else
psi1(i,j)=C1(j)*exp(-x(i)*k31(j));
psi21(i,j)=(psi1(i,j))^2;
end
end
end

%расчёт значений пси-функции для сим. ямы
for j=1:size(z2,2)
for i=1:size(x,2)
if x(i)<=0
psi2(i,j)=A2(j)*exp(x(i)*k12(j));
psi22(i,j)=(psi2(i,j))^2;
elseif (0<x(i))&&(x(i)<=l)
psi2(i,j)=B2(j)*sin(x(i)*s2(j)+delt2(j));
psi22(i,j)=(psi2(i,j))^2;
else
psi2(i,j)=C2(j)*exp(-x(i)*k32(j));
psi22(i,j)=(psi2(i,j))^2;
end
end
end

set(0,'DefaultAxesFontSize',18,'DefaultAxesFontName','Times New Roman');
set(0,'DefaultTextFontSize',18,'DefaultTextFontName','Times New Roman');

if (mod(size(z1,2),2)==0)
a1=size(z1,2)/2+1;
else
a1=(size(z1,2)-1)/2+1;
end

if (mod(size(z2,2),2)==0)
a2=size(z2,2)/2+1;
else
a2=(size(z2,2)-1)/2+1;
end

%графики для несим. ямы
figure(1) %график E(n)
plot(1:size(z1,2),E1,'black','Linewidth',2);
title("Для несим. ямы");
xlabel('\it{}n\rm');
ylabel('\it{}E\rm(\it{}n\rm) [Дж]');

figure(2) %графики Ψ(x) при n=1..?
for j=1:size(z1,2)
subplot(a1,a1,j)
plot(x,psi1(:,j),'black','Linewidth',2);
title(['Для несим. ямы. \it{}n\rm=',num2str(j)])
xlabel('\it{}x\rm [м]');
ylabel('\it{}Ψ\rm(\it{}x\rm)');
end

figure(3) %графики |Ψ|^2(x) при n=1..?
for j=1:size(z1,2)
subplot(a1,a1,j)
plot(x,psi21(:,j),'black','Linewidth',2);
title(['Для несим. ямы. \it{}n\rm=',num2str(j)])
xlabel('\it{}x\rm [м]');
ylabel('\it{}|Ψ|\rm^2(\it{}x\rm)');
end

%графики для сим. ямы
figure(4) %график E(n)
plot(1:size(z2,2),E2,'black','Linewidth',2);
title("Для сим. ямы");
xlabel('\it{}n\rm');
ylabel('\it{}E\rm(\it{}n\rm) [Дж]');

figure(5) %графики Ψ(x) при n=1..?
for j=1:size(z2,2)
subplot(a2,a2,j)
plot(x,psi2(:,j),'black','Linewidth',2);
title(['Для сим. ямы. \it{}n\rm=',num2str(j)])
xlabel('\it{}x\rm [м]');
ylabel('\it{}Ψ\rm(\it{}x\rm)');
end

figure(6) %графики |Ψ|^2(x) при n=1..?
for j=1:size(z2,2)
subplot(a2,a2,j)
plot(x,psi22(:,j),'black','Linewidth',2);
title(['Для сим. ямы. \it{}n\rm=',num2str(j)])
xlabel('\it{}x\rm [м]');
ylabel('\it{}|Ψ|\rm^2(\it{}x\rm)');
end

%сим. яма с уровнями E(n)
W=zeros();
E=zeros();
for j=1:size(z2,2)
for i=1:size(x,2)
if x(i)<=0
x(i);
W(i,j)=W1;
E(i,j)=E2(j);
elseif (0<x(i))&&(x(i)<=l)
x(i);
W(i,j)=W2;
E(i,j)=E2(j);
else
W(i,j)=W1;
E(i,j)=E2(j);
end
end
end

%несим. яма с уровнями E(n)
w=zeros();
e=zeros();
for j=1:size(z1,2)
for i=1:size(x,2)
if x(i)<=0
x(i);
w(i,j)=W3;
e(i,j)=E1(j);
elseif (0<x(i))&&(x(i)<=l)
x(i);
w(i,j)=W2;
e(i,j)=E1(j);
else
w(i,j)=W1;
e(i,j)=E1(j);
end
end
end

figure(7)
plot(x,W,'black',x,w,'blue',x,E,'red',x,e,'green','LineWidth',2);
title('(Не)сим. яма с уровнями \it{}E\rm(\it{}n\rm)');
xlabel('\it{}x\rm [м]');
ylabel('\it{}W\rm(\it{}x\rm) [Дж]');
%legend('Несим. яма','Сим. яма','Энерг. уровни сим. ямы','Энерг. уровни несим. ямы');