% carrier2b <==(望月 00y17) dx を 1 にする
% carrier2a <==(望月 00y15) ドリフト電流の取り扱いを更に変える。(x固定)
% carrier2 <==(望月 00y13) ドリフト電流の取り扱いを変える。(xについてゆく)
% carrier (New 望月 99x28) 少数キャリア連続の式を体感する
z01=0; % キャリア濃度の初期条件の分布が, 0〜1=パターン0〜1
z1=1; % 拡散現象が ... 0=無い 1=ある
d0=5; % 拡散定数
d0LL=0; % 左端の濃度を強制的に設定 1=する。0=しない。
d0L=1; % 左の窓が,0=開放 1=密封
d0R=1; % 右の窓 同上
z2=1; % ドリフトが ... 0=無い 1=ある (多少の拡散分を含みます)
z3=0; % キャリアの発生が ... 0=無い 1=ある
z4=0; % キャリアの消滅が ... 0=無い 1=ある
tau=1000; % キャリア消滅の時定数
p00=1; % 熱平衡状態でのキャリア濃度
iiend=500; % シミュレーションの最終時刻
%
if z01==0,
p0=[ linspace(0,0,5) ...
linspace(2,2,45) ...
linspace(0.5,0.5,50) ...
linspace(1,1,100) ...
linspace(1,1.2,20) ...
linspace(1.2,1.2,5) ...
linspace(1.2,1,20) ...
linspace(1,1,155) ]; % 濃度
else
p0=linspace(1,1,30); % 濃度
end
p=p0;
len=length(p);
len1m=len-1;
len2m=len-2;
x0=(0:len1m);
x=x0;
% pt00=[ linspace(0,0,20) linspace(1,1,10) linspace(0,0,20)]
% pt0=[pt00 pt00 pt00 pt00 pt00]
e=[ linspace(0,0,150) ...
linspace(0,1,30) ...
linspace(1,1,90) ...
linspace(1,0,30) ...
linspace(0,0,len1m) ];
e=e(1:len);
g=[ linspace(0,0,340) ...
linspace(1,1,10) ...
linspace(0,0,len1m) ];
g=g(1:len);
%
plot(x,p,'g',x,e*0.1-0.2,'m',x,g*0.1-0.4,'y','LineWidth',2)
pold=p;
fi=1;
i=0;
for ii=20:20:iiend
for i=i+1:ii
dx=diff(x); % 区間の間隔
if z1==1, % 拡散現象
for j=1:d0
p(1)=p(2)*d0L;
p(len)=p(len1m)*d0R;
dp=diff(p)./dx./dx;
ddp=diff(dp)/20;
p=p+[0 ddp 0];
% 強制的にキャリア濃度設定
p(2)=(1-d0LL)*p(2)+d0LL*(1+sin(i/80));
end
end
if z2==1, % 電界
dx0=(x(3:len)-x(1:len2m))/2;
x=x+e/50;
dx2=(x(3:len)-x(1:len2m))/2;
p(2:len1m)=p(2:len1m).*dx0./dx2;
end
% キャリアの発生
p=p+z3*g/5000;
% キャリアの消滅
p=p-z4*(p-p00)/tau;
end
figure(1)
hold on
plot(x,p)
title(['loop ' num2str(i)])
hold off
if z2==1, % ドリフト電流存在時,xの調整
pold=p;
p1to2=0;
dx0=[1 (x0(3:len)-x0(1:len2m))/2]; % 区間の間隔
dx=[1 (x(3:len)-x(1:len2m))/2]; % 区間の間隔
ix2=2; add2=0;
for k=2:len1m
ix1=ix2; add1=add2;
while x(ix2+1)<x0(k+1);
ix2=ix2+1;
end
add2=(x0(k)-x(ix2-1))/(x(ix2)-x(ix2-1));
if ix2==ix1,
p(k)=(dx(ix1)*pold(ix1)*(add2-add1))/dx0(k);
elseif ix2==ix1+1,
p(k)=(dx(ix1)*pold(ix1)*(1-add1) ...
+ dx(ix2)*pold(ix2)*add2)/dx0(k);
elseif ix2==ix1+2,
p(k)=(dx(ix1)*pold(ix1)*(1-add1) ...
+ dx(ix1+1)*pold(ix1+1) ...
+ dx(ix2)*pold(ix2)*add2)/dx0(k);
else
pause
end
end
x=x0;
end
pause(0.1)
end
figure(1)
hold on
plot(x0,p0,'g','LineWidth',2)
plot(x,p,'r')
title(['loop ' num2str(i)])
hold off
|