#!/usr/bin/octave -qf
# consts
 ec = 1.6e-19;                          # (C)
 kB = 1.38e-23;                         # (J/K)
 NA = 6.022e23;
# params
 SM = 1.5;
 SN = 0.3;
 ac = 1.0;
 at = 1.0;
 At = 127 * 1e-4;                       # (m^2)
 Ac = 20 * At;
 Sp = 85 * 1e-3;                        # (m^3/s)
 Jt = 1.4e-3 * 1e4;                     # (A/m^2)
 Tg = 300;                              # (K)
 Mg = 2 * 14.0067 / NA * 1e-3;          # (kg)
 pN = [0.00: 0.0001: 0.03]';            # (Pa)

iF = Jt / ec;                           # ion flux: (ions/m^2s)
F = pN ./ sqrt(2 * pi * kB * Tg * Mg);  # (molecules/m^2s): eq(3)

t1 = 2*at*F ./ (2*at*F + iF*SN);        # eq(4)
iR = 2*ac*F + iF*SN*t1*At/Ac;           # eq(5)
t2 = iR ./ (iR + iF*SM*(1-t1)*At/Ac);

qt = at * F .* (1-t1) * At;             # (molecules/s): eq(6)
qc = ac * F .* (1-t2) * Ac;             # (molecules/s): eq(7)
qp = pN / (kB * Tg) * Sp;               # (molecules/s): eq(2)
q0 = qt + qc + qp;
q0 = q0 / (NA / 22.4e3 / 60);           # (molecules/s) -> (sccm)

RS = iF * (SN*t1 + SM*(1-t1));          # (a.u.)

#
# plot
# It is very much different since octave3
#
setenv("GNUTERM","postscript");
plot(q0, pN)
xlabel("N2-FLOW (sccm)");
ylabel("N2-Pp (Pa)");
print -deps berg-qP.eps
#
plot(q0, RS)
xlabel("N2-FLOW (sccm)");
ylabel("Sputt. Ratio (a.u.)");
print -deps berg-qR.eps

