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;