pn_node.m

％ 用于模拟PN结中载流子的分布、电场和电势 ％ 总长为 L 的半导体截面上

% ---------------------------------------------------------------------------------------------------------------
% 用于模拟PN结中载流子的分布、电场和电势
% 总长为 L 的半导体截面上，左侧为 N 区，掺杂浓度 NA；右侧为 P 区，掺杂浓度为 ND
% 假设截面积足够大
% ---------------------------------------------------------------------------------------------------------------

L = 100;                                                    % 总长度
OnesL = ones(1,L);                                          % 计算用的临时变量

MatSumLeft = eye(L);                                        % 用于计算电场的矩阵，左下角为-1，右上角为1，对角线为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;                                                % 假设的 Ni2 ＝ N×P 恒定数值，本征浓度 Ni ＝ 100

NA = 2000;                                                  % 初始状态右侧 P 区空穴浓度，约等于掺杂浓度
ND = 5000;                                                  % 初始状态左侧 N 区电子浓度，约等于掺杂浓度

CN_Diffuse = 0.01;                                          % 扩散系数
CP_Diffuse = 0.01;                                          % 扩散系数 

CN_Drift = 0.0000001;                                       % 漂移系数
CP_Drift = 0.0000001;                                       % 漂移系数

N = [ones(1,L/2)*ND ones(1,L/2)*Ni2/NA];                    % 初始状态，按 P*N = Ni2 得到
P = [ones(1,L/2)*Ni2/ND ones(1,L/2)*NA];                    % 初始状态，按 P*N = Ni2 得到

Grid = N - P;                                               % 晶格电荷，真正的掺杂浓度

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

for i=1:10000,                                              % 步进循环
    
    % 计算电场，并根据电场计算漂移量-----------------------------------------------------------------------------

    E = ( Grid + P - N ) * MatSumE;                         % 无限大截面假设下，电场正比于电荷
                                                     
    LeftE = (E < 0) .* abs(E);                              % 计算用的辅助判断矢量
    RightE = (E > 0) .* abs(E);                             % 计算用的辅助判断矢量

    dP_Drift_Left = LeftE  * CP_Drift;                      % 空穴受到向左的电场而产生的漂移量（比例）
    Temp = (dP_Drift_Left - 1) <= 0;                        % 做限幅运算，避免漂移量大于原有空穴量
    dP_Drift_Left = dP_Drift_Left .* Temp + OnesL - Temp;   % 
    dP_Drift_Left(1) = 0;                                   % 最左侧的位置处，空穴不能向左继续漂移
    
    dP_Drift_Right = RightE  * CP_Drift;                    % 空穴受到向右的电场而产生的漂移量（比例）
    Temp = (dP_Drift_Right - 1) <= 0;                       % 做限幅运算，避免漂移量大于原有空穴量
    dP_Drift_Right = dP_Drift_Right .* Temp + OnesL - Temp; % 
    dP_Drift_Right(L) = 0;                                  % 最右侧的位置处，空穴不能向右继续漂移

    dP_Drift_Reduce = dP_Drift_Left + dP_Drift_Right;       % 任意位置处总共减少的电荷量（比例）

    dP_Drift_Left = dP_Drift_Left .* P;                     % 任意位置处向左漂移的空穴量（绝对数值）
    dP_Drift_Right = dP_Drift_Right .* P;                   % 任意位置处向右漂移的空穴量（绝对数值）
    dP_Drift_Reduce = dP_Drift_Reduce .* P;                 % 任意位置处总共减少的电荷量（绝对数值）
    
    
    dN_Drift_Left = RightE  * CN_Drift;                     % 以下是电子的运算，完全类似于空穴部分
    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;   % 计算调整量
    dN_Drift = [0,dN_Drift_Right(1:L-1)] + [dN_Drift_Left(2:L),0] -  dN_Drift_Reduce;   % 计算调整量

    N = N + dN_Drift;                                       % 计算调整后的电子数量
    P = P + dP_Drift;                                       % 计算调整后的空穴数量
    
    % 计算扩散量-------------------------------------------------------------------------------------------------

    
    dN_Diffuse = [N(1) N(1:L-1)] + [N(2:L) N(L)] - 2*N;     % 两侧的梯度差决定了电子的增减量       
    dP_Diffuse = [P(1) P(1:L-1)] + [P(2:L) P(L)] - 2*P;     % 两侧的梯度差决定了空穴的增减量

    N = N + dN_Diffuse * CN_Diffuse;                        % 计算调整后的电子数量
    P = P + dP_Diffuse * CP_Diffuse;                        % 计算调整后的空穴数量

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

    d_PN = sqrt(Ni2 ./ (P .* N));                           % 根据 Ni2 / （P×N） 计算修正量是 大于1 或 小于 1
    
    P_Larger = P > N;                                       % 计算用的辅助判断矢量
    N_Larger = OnesL - P_Larger;                            % 计算用的辅助判断矩阵
    
    dP_Modify = N_Larger .* (d_PN - 1) .* P;                % 如果 P < N ，则令 P = P *  sqrt(Ni2 / （P×N))
    P = P + dP_Modify;                                      %
    N = N + dP_Modify;                                      % 并对 N 做相应的调整
    
    dN_Modify = P_Larger .* (d_PN - 1) .* N;                % 如果 N < P ，则令 N = N *  sqrt(Ni2 / （P×N))
    N = N + dN_Modify;                                      % 
    P = P + dN_Modify;                                      % 并对 P 做相应的调整
    
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;                        % 绘制 电场曲线
subplot(3,1,3);plot(1:L,E*MatSumRight); grid on;            % 绘制 电势曲线