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;