CONTENT:
program A_tools; {Hannes.Tammet@ut.ee 201200205}
{$APPTYPE CONSOLE}
uses windows; // delete this line when outside of Windows or Delphi!
const title = ' A_tools 20120121';
{An overview of included procedures is available in web:
http://ael.physic.ut.ee/tammet/a_tools.htm
"A toolbox for air ion and aerosol calculations"}
procedure initconsole; // delete this procedure when outside of Delphi!
var i : integer;
crd : coord;
rct : small_rect;
consolehandle : thandle;
begin
consolehandle := GetStdHandle (STD_OUTPUT_HANDLE);
rct.left := 0; rct.right := 79;
rct.top := 0; rct.bottom := 50;
SetConsoleWindowInfo (consolehandle, true, rct);
SetConsoleTextAttribute (consolehandle, 240);
setconsoletitle (title);
for i := 1 to 1200 do write (' ');
crd.x := 0; crd.y := 0;
SetConsoleCursorPosition (consolehandle, crd);
end;
function log10 (x : real) : real; begin log10 := ln (x) / ln (10) end;
function power (x, k : real) : real; begin power := exp (k * ln (x)); end;
{A. PARTICLE SIZE AND MOBILITY IN AIR}
function electrical_mobility_air ( {Tammet95m}
Millibar,
Celsius,
ParticleDensity {g cm-3, for cluster ions typically 2.08},
ParticleCharge {e, for cluster ions 1},
MassDiameter {nm} : real) : real;
function Omega11 (x : real) : real; {*(1,1)*(T*) for (*-4) potential}
var p, q : real; {and elastic-specular collisions}
begin
if x > 1 then Omega11 := 1 + 0.106 / x + 0.263 / exp ((4/3) * ln (x))
else begin p := sqrt (x); q := sqrt (p);
Omega11 := 1.4691 / p - 0.341 / q + 0.181 * x * q + 0.059 end;
end;
const GasMass = 28.96; {amu} Polarizability = 0.00171; {nm3}
a = 1.165; b = 0.48; c = 1.001; {the slip factor coefficients}
ExtraDistance = 0.115 {nm}; TransitionDiameter = 2.48 {nm};
var Temperature, GasDiameter, Viscosity, FreePath, DipolEffect,
DeltaTemperature, CheckMark, ParticleMass, CollisionDistance,
Kn, Omega, s, x, y, r1, r2 : real;
begin
if MassDiameter < 0.2 then {emergency exit}
begin electrical_mobility_air := 1e99; exit; end;
Temperature := Celsius + 273.15;
r1 := Temperature / 296.15; r2 := 406.55 / (Temperature + 110.4);
Viscosity {microPa s} := 18.3245 * r1 * sqrt (r1) * r2;
FreePath {nm} := 67.3 * sqr (r1) * (1013 / millibar) * r2;
ParticleMass {amu} := 315.3 * ParticleDensity * exp (3 * ln (MassDiameter));
DeltaTemperature := Temperature;
repeat
CheckMark := DeltaTemperature;
GasDiameter {nm} := 0.3036 * (1 + exp (0.8 * ln (44 / DeltaTemperature)));
CollisionDistance {nm} := MassDiameter / 2 + ExtraDistance +
GasDiameter / 2;
DipolEffect := 8355 * sqr (ParticleCharge) * Polarizability /
sqr (sqr (CollisionDistance));
DeltaTemperature := Temperature + DipolEffect;
until abs (CheckMark - DeltaTemperature) < 0.01;
if ParticleCharge = 0 then Omega := 1
else Omega := Omega11 (Temperature / DipolEffect);
Kn := FreePath / CollisionDistance;
if Kn < 0.03 {underflow safe} then y := 0 else y := exp (- c / Kn);
x := (273.15 / DeltaTemperature) *
exp (3 * ln (TransitionDiameter / MassDiameter));
if x > 30 {overflow safe} then s := 1
else if x > 0.001
then s := 1 + exp (x) * sqr (x / (exp (x) - 1)) * (2.25 / (a + b) - 1)
else {underflow safe} s := 1 + (2.25 / (a + b) - 1);
electrical_mobility_air := 1.6022 * ParticleCharge *
((2.25 / (a + b)) / (Omega + s - 1)) *
sqrt (1 + GasMass / ParticleMass) *
(1 + Kn * (a + b * y)) /
(6 * PI * Viscosity * CollisionDistance);
end;
function mass_diameter_air ( {Tammet95m}
millibar,
Celsius,
density {g cm-3, for cluster ions typically 2},
charge {e, for cluster ions 1},
Z {cm2 V-1 s-1} : real) : real;
{Uses external function "electrical_mobility_air"}
var y0, dx, d1, d2 : real;
function y (d : real) : real;
begin // nnn := nnn + 1;
y := 1.6022 * charge /
electrical_mobility_air (millibar, Celsius, density, charge, d);
end;
begin
y0 := 1.6022 * charge / Z;
d1 := ((0.6 + sqrt (0.36 + 200 * 150 / y0)) / (150 / y0) - 0.35) /
(1.1 + 0.001 * z * z * z * charge * charge);
d2 := d1;
while y (d1) > y0 do d1 := d1 / 1.2;
while y (d2) < y0 do d2 := d2 * 1.2;
repeat dx := (d1 + d2) / 2;
if y (dx) > y0 then d2 := dx else d1 := dx;
until abs (d2 - d1) / d1 < 1e-6;
mass_diameter_air := (d2 + d1) / 2;
end;
function electrical_mobility_1_air (millibar, Celsius, d {nm} : real) : real;
{cm2 V-1 s-1, singly charged particles of density 2 g/cm3 in air,
uses external function "electrical_mobility_air"}
begin electrical_mobility_1_air :=
electrical_mobility_air (millibar, Celsius, 2, 1, d);
end;
function mass_diameter_1_air (millibar, Celsius, Z {cm2 V-1 s-1} : real) : real;
{nm, singly charged particles of density 2 g/cm3 in air,
uses external function "mass_diameter_air"}
begin mass_diameter_1_air :=
mass_diameter_air (millibar, Celsius, 2, 1, Z);
end;
function Z_approx_air (millibar, Celsius, d {nm} : real) : real;
{cm2 V-1 s-1, mathematical approximation of mobility Tammet95m}
const a = 1.165; b = 0.48; c = 1.001; {Jung et al., 2011}
var Viscosity, FreePath, Kn, r1, r2, y, dd : real;
begin
y := d + 0.001 * Celsius - 1.2;
dd := d + 0.6 - 0.22 * exp (-1.8 * y * y) - 0.033 / d;
r1 := (Celsius + 273.15) / 296.15;
r2 := 406.55 / (Celsius + 383.55);
Viscosity {microPa s} := 18.3245 * r1 * sqrt (r1) * r2;
FreePath {nm} := 67.3 * sqr (r1) * (1013 / millibar) * r2;
Kn := 2 * FreePath / dd;
if Kn < 0.03 {underflow safe} then y := 0 else y := exp (- c / Kn);
Z_approx_air :=
1.6022 * (1 + Kn * (a + b * y)) / (3 * PI * Viscosity * dd);
end;
function d_approx_air (millibar, Celsius, Z {cm2V-1s-1} : real) : real;
{Uses external function "Z_approx_air". If failed then result = 0}
var c, d, test : real;
n : integer;
begin
c := 100; n := 0;
repeat
n := n + 1;
d := (0.6 + sqrt (0.36 + 200 * c * Z)) / (c * Z) - 0.4;
test := Z_approx_air (millibar, Celsius, d);
c := (1.2 / (d + 0.4) + 200 / sqr (d + 0.4)) / test;
until (abs (test / Z - 1) < 1e-6) or (n = 99);
if n < 99 then d_approx_air := d else d_approx_air := 0;
end;
function plain_Millikan_mobility (millibar, Celsius, d {nm} : real) : real;
{cm2 V-1 s-1, according to ISO15900 without diameter correction,
the correction 0.3-0.6 can be added to the input value of d}
const a = 1.165; b = 0.48; c = 1.001; {Jung et al, 2011}
var Kelvin, Viscosity, FreePath, Kn, r1, r2, y : real;
begin
Kelvin := Celsius + 273.15;
r1 := Kelvin / 296.15; r2 := 406.55 / (Kelvin + 110.4);
Viscosity {microPa s} := 18.3245 * r1 * sqrt (r1) * r2;
FreePath {nm} := 67.3 * sqr (r1) * (1013 / millibar) * r2;
Kn := 2 * FreePath / (d);
if Kn < 0.03 {underflow safe} then y := 0 else y := exp (- c / Kn);
plain_millikan_mobility := 1.6022 * (1 + Kn * (a + b * y)) /
(3 * PI * Viscosity * d);
end;
function plain_Millikan_diameter (millibar, Celsius, Z {cm2V-1s-1} : real) : real;
{Uses external function "plain_Millikan_mobility". If failed then result = 0.
The diameter correction 0.3-0.6 can be added to the returned value of d}
var c, d, test : real;
n : integer;
begin
c := 100; n := 0;
repeat
n := n + 1;
d := (0.6 + sqrt (0.36 + 200 * c * Z)) / (c * Z) - 0.4;
test := plain_Millikan_mobility (millibar, Celsius, d);
c := (1.2 / (d + 0.4) + 200 / sqr (d + 0.4)) / test;
until (abs (test / Z - 1) < 1e-6) or (n = 99);
if n < 99 then plain_Millikan_diameter := d else plain_Millikan_diameter := 0;
end;
{B. ATTACHMENT, COAGULATION AND CHARGE DISTRIBUTION}
function attachment_coefficient {Ion-to-particle, ions are singly charged}
{cm3/s} (q, {number of charges on particle, attracting -, repelling +}
di, {diameter of ion : nm}
gcm3, {ion density : g cm-3, typically 2.08}
dp, {diameter of particle : nm}
T, {temperature : K}
p {pressure : mb} : real) : real;
{Uses external function "electrical_mobility_air"
Comment: diameters of atmospheric ions of mobility 1.36 and 1.56 cm2/Vs,
are in standard conditions 0.79 and 0.70 nm}
var d, bi, bp, pi2Di, pi2Dp, x, f, y, bb : real;
begin
bi := electrical_mobility_air (p, T - 273.15, gcm3, 1, di)/1.6022;
bp := electrical_mobility_air (p, T - 273.15, gcm3, abs (q), dp)/1.6022;
pi2Di := 8.674e-8 * T * bi; {2*pi*D : SI}
pi2Dp := 8.674e-8 * T * bp; {2*pi*D : SI}
d := dp + di + 0.5; {nm}
x := 2 * q * (1.671e4 / T) / d;
y := 1 - 1 / (1 + 0.5 * q * (q - 1) + d / 23);
if x < -33 then f := -x
else if abs (x) < 1e-6 then f := 1
else if x > 33 then f := 0
else f := x / (exp (x) - 1);
bb := 1e-3 * y * d * f * sqrt (sqr (pi2Dp) + sqr (pi2Di));
if (abs (q + 1) < 0.001) and ((bb < (1.6e-6 * bi)) or (dp < 5))
then bb := 1.6e-6 * bi;
attachment_coefficient {cm3/s} := bb;
end;
function coagulation_coefficient {Particle-to-particle, Sahni approximation}
{cm3/s} (d1, {diameter of small neutral particle : nm}
d2, {diameter of large charged or neutral particle : nm}
gcm3, {small particle density : g cm-3, typically 1.5-2}
q, {large particle charge : e}
aa, {small particle polarizability : angstrom^3}
T, {air temperature : K}
mb {air pressure : mb} : real) : real;
{Uses external function "electrical_mobility"}
{Comment: polarizability in angstrom^3 is often estimated as equal to
the number of atoms in the cluster or as r^3 for particles}
const {that could be edited or moved to the list of parameters}
h = 0.115; {extra distance : nm}
d0 = 2.5; {critical diameter of quantum rebound : nm}
T0 = 300; {extra temperature of quantum rebound: K}
var x, y, d, b, m1, m2, g, p, a, fe, pp, qq : real;
begin
x := 8.674e-23 * T; {2*pi*k*T : SI}
d := 1e-9 * (d1 + d2 + 2 * h); {2 * delta : SI}
m1 := 5.236e-25 * gcm3 * d1 * d1 * d1; {SI}
m2 := 5.236e-25 * gcm3 * d2 * d2 * d2; {SI}
b := 1e15 * {B1 + B2 : SI}
(electrical_mobility_air (mb, T, gcm3, 0, d1) +
electrical_mobility_air (mb, T, gcm3, 0, d2)) / 1.6022;
g := 2 * (b / d) * sqrt (x * m1 * m2 / (m1 + m2)); {gamma = Kc / Kk}
if d0 = 0 then p := 1 else begin
y := (273 / (T + T0)) * exp (3 * ln (d0 / d2));
p := sqr (y) * exp (y) / sqr (exp (y) - 1); {sticking probability}
end;
if (q = 0) or (aa = 0) then fe := 1 else begin
a := aa * 1e-30; {polarizability : m^3}
y := 119700 * T * sqr (sqr (d / 2)) / (a * sqr (q));
{y = T* = kT/Upol, 119700 = 8 * pi * eps0 * k / e^2}
if y > 1 then fe := 1 + 0.106 / y + 0.263 / exp ((4/3) * ln (y))
else begin pp := sqrt (y); qq := sqrt (pp);
fe := 1.4691 / pp - 0.341 / qq + 0.182 * y * qq + 0.059;
end;
end;
coagulation_coefficient {cm3/s} := 1e6 * fe * x * b * d /
(1 + g - 0.299 * g / (exp (1.1 * ln (g)) + 0.64) + g * (1 - p) / p);
end;
function Boltzmann_probability {of charge qe on a particle of diameter d}
(q {number of elementary charges}: integer;
d {particle diameter : nm},
T {air temperature : K} : real) : real;
const k = 1.38e-23; {Boltzmann constant}
e = 1.6e-19; {elementary charge}
eps = 8.85e-12; {electric constant}
var i : integer;
a, x, s : real;
begin
a := e * e / (4 * pi * eps * d * 1e-9 * k * T);
i := 0; s := 0;
repeat i := i + 1;
x := exp (-i * i * a);
s := s + x;
until x < 1e-9;
Boltzmann_probability := exp (-q * q * a) / ( 1 + 2 * s);
end;
function Wiedensohler_probability {of charge qe on a particle of diameter dp}
(q {number of elementary charges}: integer;
d {particle diameter : nm}: real) : real;
{According to Wiedensohler et al. (2010) AMTD 3, 5521-5587,
valid for q in [-2..2]}
const Wiedensohler : array [0..5, -2..2] of real =
((-26.3328, -2.3197, -0.0003, -2.3484, -44.4756),
(35.9044, 0.6175, -0.1014, 0.6044, 79.3772),
(-21.4608, 0.6201, 0.3073, 0.4800, -62.8900),
(7.0867, -0.1105, -0.3372, 0.0013, 26.4492),
(-1.3088, -0.1260, 0.1023, -0.1553, -5.7480),
(0.1051, 0.0297, -0.0105, 0.0320, 0.5049));
var i : integer;
logd, sum : real;
begin
if (q < -2) or (q > 2) then begin
writeln ('Illegal arguments in charge_probability_W, press ENTER!');
readln; halt;
end;
logd := log10 (d);
sum := Wiedensohler [0, q];
for i := 1 to 5 do sum := sum + Wiedensohler [i, q] * power (logd, i);
Wiedensohler_probability := power (10, sum);
end;
function charge_probability {of charge qe on a particle of diameter d}
(q0 {number of elementary charges}: integer;
dp {particle diameter : nm},
di {ion diameter : nm},
T {air temperature : K} : real) : real;
{Uses external function "attachment_coefficient", presumtions:
+- symmetry, ion density 2 g cm-3, air pressure 1013 mb}
var j, q : integer;
Nq, Nj, c, sum : real;
begin
q := abs (q0);
Nj := 1; c := 1; sum := 1; j := 0; Nq := 0;
while (j <= q) or (c > 1e-12) do begin
if j = q then Nq := Nj;
c := attachment_coefficient (j, di, 2, dp, T, 1013) /
attachment_coefficient (-j - 1, di, 2, dp, T, 1013);
Nj := Nj * c;
sum := sum + 2 * Nj;
j := j + 1;
end;
charge_probability := Nq / sum;
end;
procedure ion_equilibrium {in free air without nucleation}
{arguments:} (Z_pos {mobility of positive small ions : cm2 V-1 s-1},
Z_neg {mobility of negative small ions : cm2 V-1 s-1},
i_rate {ionization rate cm-3s-1},
T {air temperature : K},
n_particle {number concentration of background aerosol particles : cm-3},
d_particle {mean diameter of background aerosol particles : nm} : real;
{results:} var n_pos {concentration of positive small ions, cm-3},
n_neg {concentration of negative small ions, cm-3} : real);
{charges carried by background aerosol particles = n_neg - n_pos}
const a = 1.5e-6; {recombinationcoefficient : cm3 s-1}
eps = 1e-3; {tolerance}
var dia_concentration, correct,
pi2D_pos, pi2D_neg,
sink_p, sink_n, sink_a,
nip, nin, dip, din,
dt, dq: real;
begin {internal calculations in cm and s}
dq := 1.671e-3 / T; {Coulomb length : cm}
pi2D_pos := 5.414e-4 * T * Z_pos; {2*pi*D : cm2 s-1}
pi2D_neg := 5.414e-4 * T * Z_neg; {2*pi*D : cm2 s-1}
dia_concentration := 1e-7 * (d_particle - 1.5) * n_particle; {cm-2}
correct := (d_particle + 9) / (d_particle + 23);
{first approximation}
sink_a := 0.5 * (pi2D_pos + pi2D_neg) * dia_concentration;
nip := (sqrt (sink_a * sink_a + 4 * a * i_rate) - sink_a) / (2 * a);
nin := nip;
{iteration}
dt := 1 / (2 * sink_a + (1 + i_rate) / 200);
repeat
sink_p := pi2D_pos * (dia_concentration - correct * dq * (nin - nip));
sink_n := pi2D_neg * (dia_concentration + correct * dq * (nin - nip));
dip := (i_rate - a * nip * nin - sink_p * nip) * dt; nip := nip + dip;
din := (i_rate - a * nip * nin - sink_n * nin) * dt; nin := nin + din;
until (abs (dip) < eps) and (abs (din) < eps);
n_pos := nip; n_neg := nin;
end;
{C. PARTICLE GROWTH AND SINK}
function Knudsen_growthrate {nm h-1}
(dg, {diameter of a growth unit : nm}
gcm3, {density of growth unit substance: g cm-3}
ng, {number concentration of growth units: cm-3}
T {air temperature : K} : real) : real;
{Nucleating particles are expected to grow by condensation of a substance,
which particles are called the growth units. Growth units are molecules
or clusters. Plain Knudsen growth rate assumes the free molecule regime
when growth rate does not depend on the size of the growing particle}
const k = 1.38e-23; {Boltzmann constant}
var d, m, g : real;
begin
d := 1e-9 * dg; {internal calculations in SI}
m := (pi / 6) * d * d * d * (gcm3 * 1000);
g := (pi / 12) * d * d * d * (1e6 * ng) * sqrt (8 * k * T / (pi * m)); {SI}
Knudsen_growthrate := (1e9 * 3600) * g;
end;
function growthfactor0 {dimensionless}
(d1, {diameter of small neutral particle : nm}
d2, {diameter of large charged or neutral particle : nm}
gcm3, {small particle density : g cm-3, typically 1.5-2}
q, {large particle charge : e}
aa, {small particle polarizability : angstrom^3}
T, {air temperature : K}
mb {air pressure : mb} : real) : real;
{Uses external function "electrical_mobility_air"
Presumption: condensing substance is expected not to evaporate
Growthfactor1 = (Growth rate) / (Plain Knudsen growth rate), see
"Tools for air ion and aerosol calculations" and (Tammet, Kulmala, 2005)
Applicable for growth units of d < 0.7 nm, e.g. sulphuric acid
This function is an impementation of the method of effective polarizability.
Polarizability aa in angstrom^3 is often estimated as equal to the number
of atoms in a cluster or as equal to r^3 for large particles. Polarizability
of a molecule of sulphuric acid is extra high: about 149 angstrom^3.
Alternative methods are used in functions growthfactor1 and growthfactor2}
const {that could be edited or moved to the list of parameters}
h = 0.115; {extra distance : nm}
d0 = 2.5; {critical diameter of quantum rebound : nm}
T0 = 300; {extra temperature of quantum rebound: K}
var x, y, d, b, m1, m2, g, p, a, fe, pp, qq, k, u : real;
begin
x := 8.674e-23 * T; {2*pi*k*T : SI}
d := 1e-9 * (d1 + d2 + 2 * h); {2 * delta : SI}
m1 := 5.236e-25 * gcm3 * d1 * d1 * d1; {SI}
m2 := 5.236e-25 * gcm3 * d2 * d2 * d2; {SI}
b := 6.2414e14 * {B1 + B2 : SI, 0.62414 = 1/1.6022}
(electrical_mobility_air (mb, T, gcm3, 0, d1) +
electrical_mobility_air (mb, T, gcm3, 0, d2));
g := 2 * (b / d) * sqrt (x * m1 * m2 / (m1 + m2)); {gamma = Kc / Kk}
y := (273 / (T + T0)) * exp (3 * ln (d0 / d2)); {quantum effect}
p := sqr (y) * exp (y) / sqr (exp (y) - 1); {sticking probability}
if (q = 0) or (aa = 0) then fe := 1 else begin
a := aa * 1e-30; {polarizability : m^3}
y := 119700 * T * sqr (sqr (d / 2)) / (a * sqr (q));
{y = T* = kT/Upol, 119700 = 8 * pi * eps0 * k / e^2}
if y > 1 then fe := 1 + 0.106 / y + 0.263 / exp ((4/3) * ln (y))
else begin pp := sqrt (y); qq := sqrt (pp);
fe := 1.4691 / pp - 0.341 / qq + 0.182 * y * qq + 0.059;
end;
end;
k {m3/s} := fe * x * b * d /
(1 + g - 0.299 * g / (exp (1.1 * ln (g)) + 0.64) + g * (1 - p) / p);
u := sqrt (3.516e-23 * T / m1);
growthfactor0 := 1.273e18 * k / (d2 * d2 * u);
end;
function growthfactor1 {dimensionless}
(d1, {diameter of small neutral particle : nm}
d2, {diameter of large charged or neutral particle : nm}
gcm3, {small particle density : g cm-3, typically 1.5-2}
q, {large particle charge : e}
T, {air temperature : K}
mb, {air pressure : mb}
yua, {first dipole enhancement coefficient}
yub {second dipole enhancement coefficient} : real) : real;
{Uses external function "electrical_mobility_air"
Presumption: condensing substance is expected not to evaporate
Growthfactor1 = (Growth rate) / (Plain Knudsen growth rate), see
"Tools for air ion and aerosol calculations" and (Tammet, Kulmala, 2005)
Applicable for growth units of d < 0.7 nm, e.g. sulphuric acid
This function is an impementation of the method of Nadykto and yu.
yua = sqr(f1 - 1) / (f2 - 1) and yub = ln((f1 - 1) / (f2 - 1)),
where f1 is the Nadykto-Yu dipole enhanchement factor for d = 1 nm
and f2 is the Nadykto-Yu dipole enhanchement factor for d = 2 nm.
Nadykto and Yu expected f1 = 4.35 and f2 = 1.8 for sulphuric acid
at temperature 298 K, in this case yua = 14 and yub = 1.43
Alternative methods are used in functions growthfactor0 and growthfactor2}
const {that could be edited or moved to the list of parameters}
h = 0.115; {extra distance : nm}
d0 = 2.5; {critical diameter of quantum rebound : nm}
T0 = 300; {extra temperature of quantum rebound: K}
var x, y, d, b, m1, m2, g, p, k, u : real;
begin
x := 8.674e-23 * T; {2*pi*k*T : SI}
d := 1e-9 * (d1 + d2 + 2 * h); {2 * delta : SI}
m1 := 5.236e-25 * gcm3 * d1 * d1 * d1; {SI}
m2 := 5.236e-25 * gcm3 * d2 * d2 * d2; {SI}
b := 6.2414e14 * {B1 + B2 : SI, 0.62414 = 1/1.6022}
(electrical_mobility_air (mb, T, gcm3, 0, d1) +
electrical_mobility_air (mb, T, gcm3, 0, d2));
g := 2 * (b / d) * sqrt (x * m1 * m2 / (m1 + m2)); {gamma = Kc / Kk}
y := (273 / (T + T0)) * exp (3 * ln (d0 / d2)); {quantum effect}
p := sqr (y) * exp (y) / sqr (exp (y) - 1); {sticking probability}
k {m3/s} := x * b * d /
(1 + g - 0.299 * g / (exp (1.1 * ln (g)) + 0.64) + g * (1 - p) / p);
u := sqrt (3.516e-23 * T / m1);
y := 1.273e18 * k / (d2 * d2 * u);
growthfactor1 := (1 + yua * exp (-yub * d2)) * y;
end;
function growthfactor2 {dimensionless}
(d1, {diameter of small neutral particle : nm}
d2, {diameter of large charged or neutral particle : nm}
gcm3, {small particle density : g cm-3}
q, {large particle charge : e}
aa, {small particle polarizability : angstrom^3}
T, {air temperature : K}
mb {air pressure : mb} : real) : real;
{Uses external function "electrical_mobility_air"
Presumption: condensing substance is expected to evaporate back
into air according to the nano-Koehler law
Growthfactor2 = (Growth rate) / (Plain Knudsen growth rate), see
"Tools for air ion and aerosol calculations" and (Tammet, Kulmala, 2005)
Applicable for growth units of d > 0.7 nm, e.g. organic compounds
Alternative methods are used in functions growthfactor0 and growthfactor2}
const {that could be edited or moved to the list of parameters}
h = 0.115; {extra distance : nm}
d0 = 3.0; {critical diameter of nano-Koehler model : nm}
p = 2; {power of the nano-Koehler model}
var x, y, d, b, m1, m2, g, a, fe, pp, qq, k, kk, u : real;
begin
x := 8.674e-23 * T; {2*pi*k*T : SI}
d := 1e-9 * (d1 + d2 + 2 * h); {2 * delta : SI}
m1 := 5.236e-25 * gcm3 * d1 * d1 * d1; {SI}
m2 := 5.236e-25 * gcm3 * d2 * d2 * d2; {SI}
b := 6.2414e14 * {B1 + B2 : SI, 0.62414 = 1/1.6022}
(electrical_mobility_air (mb, T, gcm3, 0, d1) +
electrical_mobility_air (mb, T, gcm3, 0, d2));
g := 2 * (b / d) * sqrt (x * m1 * m2 / (m1 + m2)); {gamma = Kc / Kk}
if (q = 0) or (aa = 0) then fe := 1 else begin
a := aa * 1e-30; {polarizability : m^3}
y := 119700 * T * sqr (sqr (d / 2)) / (a * sqr (q));
{y = T* = kT/Upol, 119700 = 8 * pi * eps0 * k / e^2}
if y > 1 then fe := 1 + 0.106 / y + 0.263 / exp ((4/3) * ln (y))
else begin pp := sqrt (y); qq := sqrt (pp);
fe := 1.4691 / pp - 0.341 / qq + 0.182 * y * qq + 0.059;
end;
end;
k {m3/s} := fe * x * b * d /
(1 + g - 0.299 * g / (exp (1.1 * ln (g)) + 0.64));
u := sqrt (3.516e-23 * T / m1);
if d2 < d0 then kk := 0 {nano-Koehler coefficient}
else kk := 1 - exp ( p * ln (d0 / d2));
growthfactor2 := 1.273e18 * k * kk / (d2 * d2 * u);
end;
function ion_needlesink {s-1}
(Z {ion mobility : cm2 V-1 s-1},
dneedle {mm, about 0.9 for Pinus Sylvestris},
L {m-2, total length of needles in 1 m^3 of the canopy},
wind {m s-1, inside of the canopy},
T {air temperature : K},
mb {air pressure : mb} : real) : real;
var nyy, Re, Sc, kTB : real;
begin
kTB := 8.62e-9 * T * Z;
nyy := 5.5e-7 * exp (1.8 * ln (T)) / mb;
Re := (wind * dneedle / nyy) / 2000;
Sc := nyy / kTB;
ion_needlesink := kTB * L * (0.3 +
0.62 * sqrt (Re) * exp (0.33333 * ln (Sc)) *
exp (0.8 * ln (1 + exp (0.625 * ln (Re / 282000)) )) /
sqrt (sqrt (1 + exp (0.66667 * ln (0.4 / Sc)) )) );
end;
function particle_needlesink {s-1}
(dparticle {nm},
gcm3 {particle density : g cm-3},
dneedle {mm, about 0.9 for Pinus Sylvestris},
L {m-2, total length of needles in 1 m^3 of the canopy},
wind {m s-1, inside of the canopy},
T {air temperature : K},
mb {air pressure : mb} : real) : real;
{Uses external function "electrical_mobility_air"}
var nyy, Re, Sc, kTB : real;
begin // 8.517e-9 = 1.38e-8 / 1.6022
kTB := 8.517e-9 * T * electrical_mobility_air (mb, T, gcm3, 0, dparticle);
nyy := 5.5e-7 * exp (1.8 * ln (T)) / mb;
Re := (wind * dneedle / nyy) / 2000;
Sc := nyy / kTB;
particle_needlesink := kTB * L * (0.3 +
0.62 * sqrt (Re) * exp (0.33333 * ln (Sc)) *
exp (0.8 * ln (1 + exp (0.625 * ln (Re / 282000)) )) /
sqrt (sqrt (1 + exp (0.66667 * ln (0.4 / Sc)) )) );
end;
function ion_aerosolsink {s-1}
(Z {ion mobility : cm2 V-1 s-1},
p {polarity of ions +1 or -1},
n_particle {cm-3, number concentration of background aerosol particles},
d_particle {nm, mean diameter of background aerosol particles},
q_particle {e, algebraic mean charge number of background aerosol particles},
T {air temperature : K} : real) : real;
var pi2D, dq, c, s : real;
begin
pi2D := 5.414e-8 * T * Z; {2*pi*D : SI}
dq := 1.671e4 / T; {Coulomb length : nm}
c := (d_particle + 9) / (d_particle + 23);
s := pi2D * ((d_particle - 1.5) - c * p * q_particle * dq) * 1e-9 *
n_particle * 1e6;
if s > 0 then ion_aerosolsink := s else ion_aerosolsink := 0;
end;
{D. PARTICLE SIZE AND MECHANICAL MOBILITY IN GASES}
function mechanical_mobility { air nitrogen }
{velocity/force} (GasMass {amu}, { 28.96 28.02 }
{ (m/s) / fN } Polarizability {nm3}, { 0.00171 0.00174 }
{ 1e12 CGS } VisCon1 {nm}, { 0.3036 0.2996 }
{JAS26, 1995} VisCon2 {K}, { 44 40 }
{pp. 459-475} VisCon3, { 0.8 0.7 }
{const a, b, c} Pressure {mb},
{corrected 2012} Temperature {K},
ParticleDensity {g cm-3, for cluster ions typically 2.08},
ParticleCharge {e, for cluster ions 1},
MassDiameter {nm} : real) : real;
{NB: Multiply with 1.6022 to get electrical mobility of
a single charged cluster or particle in cm2V-1s-1}
function Omega11 (x : real) : real; {*(1,1)*(T*) for (*-4) potential}
var p, q : real; {and elastic-specular collisions}
begin
if x > 1 then Omega11 := 1 + 0.106 / x + 0.263 / exp ((4/3) * ln (x))
else begin p := sqrt (x); q := sqrt (p);
Omega11 := 1.4691 / p - 0.341 / q + 0.181 * x * q + 0.059 end;
end;
const a = 1.165; b = 0.48; c = 1.001; {the slip factor coefficients}
// a = 1.2; b = 0.5; c = 1; {old version}
ExtraDistance = 0.115 {nm}; TransitionDiameter = 2.48 {nm};
var GasDiameter, MeanVelocity, Viscosity, FreePath, DipolEffect,
DeltaTemperature, CheckMark, ParticleMass, CollisionDistance,
Kn, Omega, s, x, y : real;
begin
if MassDiameter < 0.2 then // emergency exit
begin mechanical_mobility := 1e99; exit; end;
Viscosity {microPa s} := 0.02713 * sqrt (GasMass * Temperature) /
sqr (VisCon1 * (1 + exp (VisCon3 * ln (VisCon2 / Temperature))));
MeanVelocity {m/s} := 145.5 * sqrt (Temperature / GasMass);
FreePath {nm} := (166251 * Viscosity * Temperature) /
(GasMass * Pressure * MeanVelocity);
ParticleMass {amu} := 315.3 * ParticleDensity *
exp (3 * ln (MassDiameter));
DeltaTemperature := Temperature;
repeat
CheckMark := DeltaTemperature;
GasDiameter {nm} := VisCon1 *
(1 + exp (VisCon3 * ln (VisCon2 / DeltaTemperature)));
CollisionDistance {nm} := MassDiameter / 2 + ExtraDistance +
GasDiameter / 2;
DipolEffect := 8355 * sqr (ParticleCharge) * Polarizability /
sqr (sqr (CollisionDistance));
{NB: erratum (ParticleCharge) in JAS is corrected}
DeltaTemperature := Temperature + DipolEffect;
until abs (CheckMark - DeltaTemperature) < 0.01;
if ParticleCharge = 0 then Omega := 1
else Omega := Omega11 (Temperature / DipolEffect);
Kn := FreePath / CollisionDistance;
if Kn < 0.03 {underflow safe} then y := 0 else y := exp (- c / Kn);
{NB! erratum in JAS (y := 1) corrected!}
x := (273.15 / DeltaTemperature) *
exp (3 * ln (TransitionDiameter / MassDiameter));
if x > 30 {overflow safe} then s := 1
else if x > 0.001
then s := 1 + exp (x) * sqr (x / (exp (x) - 1)) * (2.25 / (a + b) - 1)
else {underflow safe} s := 1 + (2.25 / (a + b) - 1);
mechanical_mobility := ((2.25 / (a + b)) / (Omega + s - 1)) *
sqrt (1 + GasMass / ParticleMass) *
(1 + Kn * (a + b * y)) /
(6 * PI * Viscosity * CollisionDistance);
end;
function mechanical_diameter {nm} { air nitrogen }
(GasMass {amu}, { 28.96 28.02 }
Polarizability {nm3}, { 0.00171 0.00174 }
VisCon1 {nm}, { 0.3036 0.2996 }
{JAS26, 1995} VisCon2 {K}, { 44 40 }
{pp. 459-475} VisCon3, { 0.8 0.7 }
{updated 2012} Pressure {mb},
Temperature {K},
ParticleDensity {g cm-3, for cluster ions typically 2.08},
ParticleCharge {e, for cluster ions 1},
Mechmobility {(m/s) / fN} : real) : real;
{Uses external function "mechanical_mobility"
If failed then result = 0}
var c, d, test : real;
n : integer;
begin
c := 300; n := 0;
repeat
n := n + 1;
d := (0.6 + sqrt (0.36 + 200 * c * MechMobility)) /
(c * MechMobility) - 0.3;
test := mechanical_mobility (GasMass, Polarizability,
VisCon1, VisCon2, VisCon3,
Pressure, Temperature,
ParticleDensity, ParticleCharge, d);
c := (1.2 / (d + 0.3) + 200 / sqr (d + 0.3)) / test;
until (abs (test / MechMobility - 1) < 0.0001) or (n = 99);
if n < 99 then mechanical_diameter := d else mechanical_diameter := 0;
end;
{==============================================================================}
procedure test_mobility_air;
const d : array [1..10] of real = (0.5, 0.7, 1, 1.5, 2, 3, 5, 10, 20, 50);
var i : integer;
p, t, di,
z1, z2, z3, z4 : real;
begin
write ('Write air pressure (mb) = '); readln (p);
write ('Write air temperature (C) = '); readln (t);
writeln ('Electrical mobility according 4 functions:');
writeln (' 1.6*mech electric approx Millikan');
writeln ('d:nm Z1 Z2 Z3 ',
' Z4 Z1/Z2 Z3/Z2 Z4/Z2');
for i := 1 to 10 do begin
di := d [i];
z1 := 1.6022 * mechanical_mobility (28.96, 0.00171, 0.3036, 44, 0.8,
p, t + 273.25, 2, 1, di);
z2 := electrical_mobility_1_air (p, t, di);
z3 := Z_approx_air (p, t, di);
z4 := plain_Millikan_mobility (p, t, di + 0.3);
writeln (di:4:1, z1:11:5, z2:10:5, z3:10:5, z4:10:5,
z1/z2:7:3, z3/z2:7:3, z4/z2:7:3);
end;
writeln (' Z4 = plain_Millikan_mobility (pressure, temperature, d + 0.3 nm)');
end;
procedure test_diameter_air;
const z : array [1..10] of real = (0.001, 0.005, 0.02, 0.1,
0.2, 0.4, 0.7, 1, 1.5, 2);
var i : integer;
p, t, zi,
d1, d2, d3, d4 : real;
begin
write ('Write air pressure (mb) = '); readln (p);
write ('Write air temperature (C) = '); readln (t);
writeln ('Particle diameter according 4 functions:');
writeln (' mechdia massdia1 approx Millikan');
writeln (' Z d1 d2 d3 ',
' d4 d1/d2 d3/d2 d4/d2');
for i := 1 to 10 do begin
zi := z [i];
d1 := mechanical_diameter (28.96, 0.00171, 0.3036, 44, 0.8,
p, t + 273.25, 2, 1, zi / 1.6022);
d2 := mass_diameter_1_air (p, t, zi);
d3 := d_approx_air (p, t, zi);
d4 := plain_Millikan_diameter (p, t, zi) - 0.3;
writeln (zi:4:3, d1:10:3, d2:10:3, d3:10:3, d4:10:3,
d1/d2:7:3, d3/d2:7:3, d4/d2:7:3);
end;
writeln (' d4 = plain_Millikan_diameter (pressure, temperature, Z) - 0.3 nm');
end;
procedure test_attachment; {Compare results with Table by Hoppel and Frick}
const dh : array [1..6] of real = (1, 2, 4, 10, 100, 1000); {d : nm}
r = 2.08;
var i : integer;
p, tc, tk, di : real;
begin
write ('Write air pressure (mb) = '); readln (p);
write ('Write air temperature (C) = '); readln (tc); tk := tc + 273.15;
write ('Write ion diameter (~0.75 nm) = '); readln (di);
writeln (' 1E6 * attachment coefficient (cm3/s):');
writeln (' d -3 -2 -1 0 1 2');
for i := 1 to 6 do writeln (dh [i] :4:0,
1e6 * attachment_coefficient (-3, di, r, dh [i], tk, p):9:3,
1e6 * attachment_coefficient (-2, di, r, dh [i], tk, p):9:3,
1e6 * attachment_coefficient (-1, di, r, dh [i], tk, p):9:3,
1e6 * attachment_coefficient ( 0, di, r, dh [i], tk, p):9:3,
1e6 * attachment_coefficient ( 1, di, r, dh [i], tk, p):9:3,
1e6 * attachment_coefficient ( 2, di, r, dh [i], tk, p):9:3);
writeln ('Compare results with the table by Hoppel and Frick');
end;
procedure test_coagulation; {Compare results with Table 27 by Fuchs}
{NB: Fuchs uses Ks for selfcoagulation dn/dt = -Ks(d) n n.
Function coagulation_coefficient issues Ks for mutual coagulation
dn1/dt = -K(d1, d2) n1 n2. Let n1 = n2 = n/2. Now
d(n)/dt = d(n/2)/dt + d(n/2)/dt = -2 K(d, d) (n/2) (n/2) =
= -K(d, d) n n/2. Thus Ks(d) = K(d, d) / 2}
const fuchs : array [1..10] of real =
(1, 2, 5, 10, 20, 30, 100, 200, 500, 1000);
var i : integer;
d, y, ks, k0, p1, p2,
p3, p4, p5, p6, p7 : real;
begin
writeln ('Comparison with Table 27 by Fuchs');
writeln ('NB: Fuchs uses Ks for selfcoagulation dn/dt = -Ks(d) n n.');
writeln (' Function issues coefficient for mutual coagulation');
writeln (' dn1/dt = -K(d1, d2) n1 n2. Let n1 = n2 = n/2. Now');
writeln (' d(n)/dt = d(n/2)/dt + d(n/2)/dt = -2 K(d, d) (n/2) (n/2) =');
writeln (' = -K(d, d) n n/2. Thus Ks(d) = K(d, d) / 2');
writeln;
writeln (' r 1e10*k0 1e10*ks');
for i := 1 to 10 do begin
d := 2 * fuchs [i];
ks := coagulation_coefficient (d, d, 2, 0, 0, 296, 1013) / 2;
k0 := 4 * pi * d * 1e-7 * 1.38e-16 * 296 * 1e12 * mechanical_mobility
(28.96, 0.00171, 0.3036, 44, 0.8, 1013, 296, 2, 0, d); {CGS}
writeln (fuchs [i] :4:0, 1e10 * k0 :9:3, 1e10 * ks :9:3);
end;
writeln;
writeln ('Individual test:');
writeln ('coagulation_coefficient (d1, d2, gcm3, q2, aa, T, p)');
writeln ('Enter values of arguments (d = 0 ==> exit)');
write (' Diameter of neutral small particle (nm) = '); readln (p1);
if p1 = 0 then exit;
write (' Diameter of neutral or charged large particle (nm) = ');
readln (p2);
write (' Density of small particle (g cm-3) = '); readln (p3);
write (' Charge of large particle (e) = '); readln (p4);
write (' Polarizability of small particle (angstrom^3) = ');readln (p5);
write (' Air temperature (Kelvin) = '); readln (p6);
write (' Air pressure (mb) = '); readln (p7);
y := coagulation_coefficient (p1, p2, p3, p4, p5, p6, p7);
writeln ('Coagulation coefficient = ', y:12, ' cm3/s');
end;
procedure test_distribution;
const dh : array [1..7] of real = (2, 5, 10, 20, 50, 100, 500); {d : nm}
var i, q : integer;
d : real;
begin
writeln (' Assumptions: T = 293 K, ion diameter = 0.8 nm');
writeln (' <====== Wiedensohler ======> <== Advanced => <= Boltzmann =>');
writeln ('d:nm q=-2 q=-1 q= 0 q=+1 q=+2 q=0 q=1 q=2 q=0 q=1 q=2');
for i := 1 to 7 do begin
d := dh [i]; write (d:4:0);
for q := -2 to 2 do write (Wiedensohler_probability (q, d):6:3);
write (' ');
for q := 0 to 2 do write (charge_probability (q, d, 0.8, 293):6:3);
write (' ');
for q := 0 to 2 do write (Boltzmann_probability (q, d, 293):6:3);
writeln;
end;
end;
procedure test_equilibrium;
var p1, p2, p3, p4, p5, p6, p7, p8 : real;
begin
writeln ('Equilibrium of ions in free air without nucleation');
writeln ('ion_equilibrium (Z_pos, Z_neg, i_rate, T,');
writeln (' n_particle, d_particle, n_pos, n_neg)');
writeln ('Enter values of arguments:');
write (' Mobility of positive ions (cm2V-1s-1) = '); readln (p1);
write (' Mobility of negative ions (cm2V-1s-1) = '); readln (p2);
write (' Ionization rate (cm-3s-1) = '); readln (p3);
write (' Air temperature (Kelvin) = '); readln (p4);
write (' Number concentration of particles (cm-3) = '); readln (p5);
write (' Diameter of particles (nm) = '); readln (p6);
ion_equilibrium (p1, p2, p3, p4, p5, p6, p7, p8);
writeln ('Concentration of positive ions = ', p7:0:1, ' cm-3');
writeln ('Concentration of negative ions = ', p8:0:1, ' cm-3');
writeln ('Space charge on background aerosol = ', p8-p7:0:1, ' cm-3');
end;
procedure test_Knudsengrowth;
var p1, p2, p3, p4, y : real;
begin
writeln ('Plain Knudsen growthrate expecting the free molecule regime');
writeln ('Knudsen_growthrate (d1, gcm3, n, T)');
writeln ('Enter values of arguments:');
write (' Diameter of growth units (nm) = '); readln (p1);
write (' Density of growth units (g cm-3) = '); readln (p2);
write (' Number concentration of growth units (cm-3) = '); readln (p3);
write (' Air temperature (Kelvin) = '); readln (p4);
y := Knudsen_growthrate (p1, p2, p3, p4);
writeln ('Knudsen_growthrate = ', y:12, ' nm/h');
end;
procedure test_growthfactors;
var p1, p2, p3, p4, p5, p6, p7, p8, p9, y : real;
begin
writeln ('Growthfactor = (Growth rate) / (Plain Knudsen growth rate),');
writeln ('non-evaporating condensing substance tyically of d < 0.7 nm.');
writeln ('Please learn the accompanying document and source code.');
writeln ('growthfactor0 (d1, d2, gcm3, q2, aa, T, p)');
writeln ('growthfactor1 (d1, d2, gcm3, q2, aa, T, p, yua, yub)');
writeln ('growthfactor2 (d1, d2, gcm3, q2, aa, T, p)');
writeln ('Enter values of arguments for growthfactor0:');
write (' Diameter of neutral growth unit (nm) = '); readln (p1);
write (' Diameter of neutral or charged large particle (nm) = ');
readln (p2);
write (' Density of growth unit (g cm-3) = '); readln (p3);
write (' Number of elementary charges on large particle = ');
readln (p4);
write (' Polarizability of growth unit (angstrom^3) = '); readln (p5);
write (' Air temperature (Kelvin) = '); readln (p6);
write (' Air pressure (mb) = '); readln (p7);
write (' First dipole enhancement coefficient = '); readln (p8);
write (' Second dipole enhancement coefficient = '); readln (p9);
y := growthfactor0 (p1, p2, p3, p4, p5, p6, p7);
writeln ('Gr0 / Gr_Knudsen = ', y:0:3);
y := growthfactor1 (p1, p2, p3, p4, p6, p7, p8, p9);
writeln ('Gr1 / Gr_Knudsen = ', y:0:3);
y := growthfactor2 (p1, p2, p3, p4, p5, p6, p7);
writeln ('Gr2 / Gr_Knudsen = ', y:0:3);
end;
procedure test_ion_needlesink;
var p1, p2, p3, p4, p5, p6, y : real;
begin
writeln ('Sink of small ions on cylindical needles e.g. in forest');
writeln ('ion_needlesink (Z, dneedle, L, wind, T, p)');
writeln ('Enter values of arguments:');
write (' Mobility of ions (cm2V-1s-1) = '); readln (p1);
write (' Diameter of needles (mm) = '); readln (p2);
write (' Total length of needles in 1 m^3 (m/m3) = '); readln (p3);
write (' Wind speed (m/s) = '); readln (p4);
write (' Air temperature (Kelvin) = '); readln (p5);
write (' Air pressure (mb) = '); readln (p6);
y := ion_needlesink (p1, p2, p3, p4, p5, p6);
writeln ('Ion_needlesink = ', y:0:6, ' s-1');
end;
procedure test_particle_needlesink;
var p1, p2, p3, p4, p5, p6, p7, y : real;
begin
writeln ('Sink of particles on cylindical needles e.g. in forest');
writeln ('particle_needlesink (dparticle, gcm3, dneedle, L, wind, T, p)');
writeln ('Enter values of arguments:');
write (' Diameter of particles (nm) = '); readln (p1);
write (' Density of particles (g cm-3) = '); readln (p2);
write (' Diameter of needles (mm) = '); readln (p3);
write (' Total length of needles in 1 m^3 (m/m3) = '); readln (p4);
write (' Wind speed (m/s) = '); readln (p5);
write (' Air temperature (Kelvin) = '); readln (p6);
write (' Air pressure (mb) = '); readln (p7);
y := particle_needlesink (p1, p2, p3, p4, p5, p6, p7);
writeln ('Particle_needlesink = ', y:0:6, ' s-1');
end;
procedure test_ion_aerosolsink;
var p1, p2, p3, p4, p5, p6, y : real;
begin
writeln ('Sink of ions on large particles e.g. on background aerosol');
write ('ion_aerosolsink (mobility, polarity,');
writeln (' n_particle, d_particle, q_particle, T)');
writeln ('Enter values of arguments:');
write (' Mobility of ions (cm2V-1s-1) = '); readln (p1);
write (' Polarity of ions (+1 or -1) = '); readln (p2);
write (' Number concentration of particles (cm-3) = '); readln (p3);
write (' Diameter of particles (nm) = '); readln (p4);
write (' Algebraic mean charge number of particles = '); readln (p5);
write (' Air temperature (Kelvin) = '); readln (p6);
y := ion_aerosolsink (p1, p2, p3, p4, p5, p6);
writeln ('Sink of ions on large particles = ', y:0:6, ' s-1');
end;
procedure test_approximation;
const p = 1013; t = 10;
di : array [1..11] of real = (0.6, 0.8, 1, 1.5, 2, 3, 5, 7, 10, 20, 50);
var i : integer;
d, mill, ht95, refe, appr : real;
begin
writeln ('Electrical mobility (cm2V-1s-1) at 1013 mb and 10 C:');
writeln
(' d:nm Millikan+ Tammet95 Reference Approx T95/Ref Mkn+/Ref Appr/Ref');
for i := 1 to 11 do begin
d := di [i];
mill := plain_millikan_mobility (p, t, d + 0.3);
ht95 := 1.6022 * mechanical_mobility (28.96, 0.00171, 0.3036, 44, 0.8,
p, t + 273.15, 2, 1, d);
refe := electrical_mobility_1_air (p, t, d);
appr := Z_approx_air (p, t, d);
writeln (d:5:1, mill:10:5, ht95:10:5, refe:10:5, appr:10:5,
ht95/refe:9:4, mill/refe:9:4, appr/refe:9:4);
end;
writeln;
end;
procedure usertest; {should be written and edited by the user}
var z, d : real;
begin
// writeln ('Should be written and edited by the user');
end;
{==============================================================================}
//* MAIN PROGRAM
Var n : integer;
c : char;
Begin try {test of procedures}
initconsole; // delete this line when outside of Delphi!
repeat
writeln ('Control letters and procedure blocks:');
writeln (' A) particle size and mobility in air');
writeln (' B) attachment, coagulation and charge distribution');
writeln (' C) particle growth and sink');
writeln (' D) particle size and mechanical mobility in gases');
writeln (' U) user written test procedure usertest');
writeln (' X) exit and close the program');
write ('Enter a letter: '); readln (c); c := upcase (c);
if not (c in ['A','B','C','D','T','U']) then halt;
if c in ['A','B','C','D'] then writeln ('Control numbers and procedures:');
case c of
'A': begin
writeln (' 1) electrical_mobility_air');
writeln (' 2) mass_diameter_air');
writeln (' 3) electrical_mobility_1_air');
writeln (' 4) mass_diameter_1_air');
writeln (' 5) Z_approx_air');
writeln (' 6) d_approx_air');
writeln (' 7) plain_Millikan_mobility');
writeln (' 8) plain_Millikan_diameter');
writeln (' 0) go to the main menu');
write ('Enter a number: ');
readln (n); writeln;
if (n > 0) and (n mod 2 = 1) then test_mobility_air;
if (n > 0) and (n mod 2 = 0) then test_diameter_air;
end;
'B': begin
writeln (' 1) attachment_coefficient (ion-to-particle)');
writeln (' 2) coagulation_coefficient {particle-to-particle)');
writeln (' 3) Boltzmann_probability (of particle charges)');
writeln (' 4) Wiedensohler_probability (of particle charges)');
writeln (' 5) charge_probability (advanced approximation)');
writeln (' 6) ion_equilibrium (in free air without nucleation)');
writeln (' 0) go to the main menu');
write ('Enter a number: ');
readln (n); writeln;
case n of
1: test_attachment;
2: test_coagulation;
3..5: test_distribution;
6: test_equilibrium;
end;
end;
'C': begin
writeln (' 1) Knudsen_growthrate (free molecule regime)');
writeln (' 2) growthfactor0 (non-evaporating substance, polarizability)');
writeln (' 3) growthfactor1 (non-evaporating substance, Nadykto & Yu)');
writeln (' 4) growthfactor2 (condensing of evaporating substance)');
writeln (' 5) ion_needlesink (e.g. in forest)');
writeln (' 6) particle_needlesink (e.g. in forest)');
writeln (' 7) ion_aerosolsink (e.g. on background aerosol)');
writeln (' 0) go to the main menu');
write ('Enter a number: ');
readln (n); writeln;
case n of
1: test_Knudsengrowth;
2..4: test_growthfactors;
5: test_ion_needlesink;
6: test_particle_needlesink;
7: test_ion_aerosolsink;
end;
end;
'D': begin
writeln (' 1) mechanical_mobility (corrected version of 1995 procedure)');
writeln (' 2) mass_diameter (improved version of 1995 procedure)');
writeln (' 0) go to the main menu');
write ('Enter a number: ');
readln (n); writeln;
if n = 1 then test_mobility_air;
if n = 2 then test_diameter_air;
end;
'T': test_approximation; {a hidden procedure}
'U': usertest; {edit it yourself!}
end;
writeln;
until false;
except end End.