% ---------------------------------------------------------------------------------------------------------------
% 用于模擬PN結(jié)中載流子的分布、電場(chǎng)和電勢(shì)
% 總長(zhǎng)為 L 的半導(dǎo)體截面上，左側(cè)為 N 區(qū)，摻雜濃度 NA；右側(cè)為 P 區(qū)，摻雜濃度為 ND
% 假設(shè)截面積足夠大
% ---------------------------------------------------------------------------------------------------------------

L = 100;                                                    % 總長(zhǎng)度
OnesL = ones(1,L);                                          % 計(jì)算用的臨時(shí)變量

MatSumLeft = eye(L);                                        % 用于計(jì)算電場(chǎng)的矩陣，左下角為-1，右上角為1，對(duì)角線為0
for iii=2:L;                                                %
    MatSumLeft(iii,:) = MatSumLeft(iii,:) + MatSumLeft(iii-1,:); %
end;                                                        %
MatSumRight = ones(L) - MatSumLeft;                         %
MatSumLeft = MatSumLeft - eye(L);                           %
MatSumE = MatSumRight - MatSumLeft;                         % 

% ---------------------------------------------------------------------------------------------------------------

Ni2 = 10000;                                                % 假設(shè)的 Ni2 ＝ N×P 恒定數(shù)值，本征濃度 Ni ＝ 100

NA = 2000;                                                  % 初始狀態(tài)右側(cè) P 區(qū)空穴濃度，約等于摻雜濃度
ND = 5000;                                                  % 初始狀態(tài)左側(cè) N 區(qū)電子濃度，約等于摻雜濃度

CN_Diffuse = 0.01;                                          % 擴(kuò)散系數(shù)
CP_Diffuse = 0.01;                                          % 擴(kuò)散系數(shù) 

CN_Drift = 0.0000001;                                       % 漂移系數(shù)
CP_Drift = 0.0000001;                                       % 漂移系數(shù)

N = [ones(1,L/2)*ND ones(1,L/2)*Ni2/NA];                    % 初始狀態(tài)，按 P*N = Ni2 得到
P = [ones(1,L/2)*Ni2/ND ones(1,L/2)*NA];                    % 初始狀態(tài)，按 P*N = Ni2 得到

Grid = N - P;                                               % 晶格電荷，真正的摻雜濃度

% ---------------------------------------------------------------------------------------------------------------

for i=1:10000,                                              % 步進(jìn)循環(huán)
    
    % 計(jì)算電場(chǎng)，并根據(jù)電場(chǎng)計(jì)算漂移量-----------------------------------------------------------------------------

    E = ( Grid + P - N ) * MatSumE;                         % 無(wú)限大截面假設(shè)下，電場(chǎng)正比于電荷
                                                     
    LeftE = (E < 0) .* abs(E);                              % 計(jì)算用的輔助判斷矢量
    RightE = (E > 0) .* abs(E);                             % 計(jì)算用的輔助判斷矢量

    dP_Drift_Left = LeftE  * CP_Drift;                      % 空穴受到向左的電場(chǎng)而產(chǎn)生的漂移量（比例）
    Temp = (dP_Drift_Left - 1) <= 0;                        % 做限幅運(yùn)算，避免漂移量大于原有空穴量
    dP_Drift_Left = dP_Drift_Left .* Temp + OnesL - Temp;   % 
    dP_Drift_Left(1) = 0;                                   % 最左側(cè)的位置處，空穴不能向左繼續(xù)漂移
    
    dP_Drift_Right = RightE  * CP_Drift;                    % 空穴受到向右的電場(chǎng)而產(chǎn)生的漂移量（比例）
    Temp = (dP_Drift_Right - 1) <= 0;                       % 做限幅運(yùn)算，避免漂移量大于原有空穴量
    dP_Drift_Right = dP_Drift_Right .* Temp + OnesL - Temp; % 
    dP_Drift_Right(L) = 0;                                  % 最右側(cè)的位置處，空穴不能向右繼續(xù)漂移

    dP_Drift_Reduce = dP_Drift_Left + dP_Drift_Right;       % 任意位置處總共減少的電荷量（比例）

    dP_Drift_Left = dP_Drift_Left .* P;                     % 任意位置處向左漂移的空穴量（絕對(duì)數(shù)值）
    dP_Drift_Right = dP_Drift_Right .* P;                   % 任意位置處向右漂移的空穴量（絕對(duì)數(shù)值）
    dP_Drift_Reduce = dP_Drift_Reduce .* P;                 % 任意位置處總共減少的電荷量（絕對(duì)數(shù)值）
    
    
    dN_Drift_Left = RightE  * CN_Drift;                     % 以下是電子的運(yùn)算，完全類(lèi)似于空穴部分
    Temp = (dN_Drift_Left - 1) <= 0;                        % ...
    dN_Drift_Left = dN_Drift_Left .* Temp + OnesL - Temp;   % ...
    dN_Drift_Left(1) = 0;                                   % ...
        
    dN_Drift_Right = LeftE .*  CN_Drift;                    % ...
    Temp = (dN_Drift_Right - 1) <= 0;                       % ...
    dN_Drift_Right = dN_Drift_Right .* Temp + OnesL - Temp; % ...
    dN_Drift_Right(L) = 0;                                  % ...

    dN_Drift_Reduce = dN_Drift_Left + dN_Drift_Right;       % ...
    
    dN_Drift_Left = dN_Drift_Left .* N;                     % ...
    dN_Drift_Right = dN_Drift_Right .* N;                   % ...
    dN_Drift_Reduce = dN_Drift_Reduce .* N;                 % ...
   
    dP_Drift = [0,dP_Drift_Right(1:L-1)] + [dP_Drift_Left(2:L),0] -  dP_Drift_Reduce;   % 計(jì)算調(diào)整量
    dN_Drift = [0,dN_Drift_Right(1:L-1)] + [dN_Drift_Left(2:L),0] -  dN_Drift_Reduce;   % 計(jì)算調(diào)整量

    N = N + dN_Drift;                                       % 計(jì)算調(diào)整后的電子數(shù)量
    P = P + dP_Drift;                                       % 計(jì)算調(diào)整后的空穴數(shù)量
    
    % 計(jì)算擴(kuò)散量-------------------------------------------------------------------------------------------------

    
    dN_Diffuse = [N(1) N(1:L-1)] + [N(2:L) N(L)] - 2*N;     % 兩側(cè)的梯度差決定了電子的增減量       
    dP_Diffuse = [P(1) P(1:L-1)] + [P(2:L) P(L)] - 2*P;     % 兩側(cè)的梯度差決定了空穴的增減量

    N = N + dN_Diffuse * CN_Diffuse;                        % 計(jì)算調(diào)整后的電子數(shù)量
    P = P + dP_Diffuse * CP_Diffuse;                        % 計(jì)算調(diào)整后的空穴數(shù)量

    % 維持 P×N = Ni2--------------------------------------------------------------------------------------------

    d_PN = sqrt(Ni2 ./ (P .* N));                           % 根據(jù) Ni2 / （P×N） 計(jì)算修正量是 大于1 或 小于 1
    
    P_Larger = P > N;                                       % 計(jì)算用的輔助判斷矢量
    N_Larger = OnesL - P_Larger;                            % 計(jì)算用的輔助判斷矩陣
    
    dP_Modify = N_Larger .* (d_PN - 1) .* P;                % 如果 P < N ，則令 P = P *  sqrt(Ni2 / （P×N))
    P = P + dP_Modify;                                      %
    N = N + dP_Modify;                                      % 并對(duì) N 做相應(yīng)的調(diào)整
    
    dN_Modify = P_Larger .* (d_PN - 1) .* N;                % 如果 N < P ，則令 N = N *  sqrt(Ni2 / （P×N))
    N = N + dN_Modify;                                      % 
    P = P + dN_Modify;                                      % 并對(duì) P 做相應(yīng)的調(diào)整
    
end;

% ---------------------------------------------------------------------------------------------------------------

subplot(3,1,1);plot(1:L,N,'r-',1:L,P,'b-'); grid on;        % 繪制 P 和 N 的分布曲線
subplot(3,1,2);plot(1:L,E); grid on;                        % 繪制 電場(chǎng)曲線
subplot(3,1,3);plot(1:L,E*MatSumRight); grid on;            % 繪制 電勢(shì)曲線