% 少数キャリア連続の式を体感する
% 望月 1999/10/28

d0=10; %拡散定数
d0L=1; %左の窓が，0=密封 1=開放
d0R=0;

p0=1;

pt0=[	linspace(0,0,50) ...
	linspace(1,1,100) ...
	linspace(0,0,50) ];
% pt00=[ linspace(0,0,20) linspace(1,1,10) linspace(0,0,20)]
% pt0=[pt00 pt00 pt00 pt00 pt00]
len=length(pt0);

e=[	linspace(0,0,5) ...
	linspace(0.5,0.5,95) ...
	linspace(1,1,75) ...
	linspace(1,0,20)];
e=[e linspace(0,0,len-1-length(e)) ];

ea=abs(e);
eplus=(e+ea)/2;
eminus=(e-ea)/2;

g=[	linspace(0,0,60) ...
	linspace(1,1,10) ];
g=[g linspace(0,0,len-length(g)) ];

p=pt0;
pout=p;

i=0;
for ii=50:50:1000
   for i=i+1:ii

      for j=1:d0
        % 拡散{
        p1=[p(1)*d0L p(1:len-1)];
	p2=[p(2:len) p(len)*d0R];
        p=0.5*p+0.25*(p1 + p2); 
        % }

        % 強制的にキャリア濃度設定{
        %p(1)=1+sin(i/100);
        % }
      end

      % 電界(+拡散){
      %p1=p(1:len-1).*eplus/100;
      %p2=p(2:len).*eminus/100;
      %p=p-[p1 0]+[0 p1]+[p2 0]-[0 p2];
      % }

      % 光{
      %p=p+g*1/100; 
      % }

      % 熱平衡{
      %p=p-(p-p0)/1000;
      % }
   end

   pout=[p ; pout];
   waterfall(pout)
   disp(['loop ' num2str(i)])
   pause(1)
end
waterfall(pout)

disp('end.')
disp('PLEASE TYPE  "view([-1 -9 20])" to Change View Point.')
disp('PLEASE TYPE  "edit carrier" to Edit Program.')
disp('//')