% Berechnung des Verlaufs der Ladungsträgerdichten in einem npn-Transistor

clear all;
close all;

% Konstanten
eps0 = 8.854e-12;
k    = 8.262e-5;
e    = 1.602e-19;

% Halbleiter-Parameter
epsr = 12;
T    = 300;
NE   = 1e20;
NB   = 1e17;
NC   = 1e15;
ni   = 1.45e10;
Ln   = 50e-6;
Lp   = 50e-6;

% Weiten der Bahngebiete und RLZen
xje = 2e-6;
xjc = 4e-6;
xk  = 8e-6;

% UT, UD bestimmen
UT  = k * T;
UDe = UT * log(NB * NE / ni^2);
UDc = UT * log(NB * NC / ni^2);

% Raster festlegen
deltax = 1e-4 * xk;

% für thd. Glg. und AP bestimmen
for i = 0:1
  UBE = i * 0.8;
  UBC = i * -1;

  % Weiten bestimmen
  wbe = sqrt(2 * epsr * eps0 * (UDe - UBE) / e * NE / NB * 1 / (NE + NB) * 1e-6);
  wbc = sqrt(2 * epsr * eps0 * (UDc - UBC) / e * NB / NC * 1 / (NC + NB) * 1e-6);
  we  = xje - wbe * NB / NE;
  x_e = xje + wbe;
  x_c = xjc - wbc * NC / NB;
  wb  = x_c - x_e;
  wc  = xk - xjc - wbc;

  % runden
  wbe = round(wbe / deltax) * deltax;
  wbc = round(wbc / deltax) * deltax;
  we  = round(we / deltax) * deltax;
  wb  = round(wb / deltax) * deltax;
  wc  = round(wc / deltax) * deltax;

  % Bereiche des Transistors
  xe    = 0:deltax:we;
  xbe   = 0:deltax:wbe;
  xb    = 0:deltax:wb;
  xbc   = 0:deltax:wbc;
  xc    = 0:deltax:wc;

  % Potential berechnen
  phibe = e * NB * 1e6  / (2 * epsr * eps0) * (wbe - xbe).^2;
  phibc = -e * NC * 1e6 / (2 * epsr * eps0) * (xbc - wbc).^2 + UDc - UBC;

  % Ladungsträgerdichten bestimmen
  ne  = NE * ones(1, length(xe));
  pe0 = ni^2 / NE * ones(1, length(xe));
  pe  = pe0 + pe0 * (exp(UBE / UT) - 1) .* sinh(xe / Lp) / sinh(we / Lp);
  nbe = NE * exp((phibe - UDe + UBE) / UT);
  pbe = NB * exp(-phibe / UT);
  nb0 = ni^2 / NB * ones(1, length(xb));
  nb  = nb0 + nb0 * (exp(UBE / UT) - 1) .* sinh((wb - xb + eps) / Ln) / sinh(wb / Ln)...
            + nb0 * (exp(UBC / UT) - 1) .* sinh(xb / Ln) / sinh(wb / Ln);
  pb  = NB * ones(1, length(xb));
  nbc = NC * exp((phibc - UDc + UBC) / UT);
  pbc = NB * exp(-phibc / UT);
  nc  = NC * ones(1, length(xc));
  pc0 = ni^2 / NC * ones(1, length(xc));
  pc  = abs(pc0 + pc0 * (exp(UBC / UT) - 1) .* sinh((wc - xc + eps) / Lp) / sinh(wc / Lp));

  x = [xe xbe+we xb+we+wbe xbc+we+wbe+wb xc+we+wbe+wb+wbc];
  n = [ne nbe nb nbc nc];
  p = [pe pbe pb pbc pc];

  figure();
  semilogy(x, n, '-', x, p, '--');
  grid on;
  legend('n(x)', 'p(x)');
  ylabel('Ladungstraegerdichte / cm^{-3}');
  xlabel('x / m');
  title(['Ladungstraegerverteilung in RLZ bei U_{BE} = ' num2str(UBE) ' V, U_{BC} = '...
        num2str(UBC) ' V']);
end
