interface

type
 myVector2f = record
   x, y: real;
 end;
 myVector2i = record
   x, y: integer;
 end;

 myBCell = record
   btm: real;                     // отметка "дна" в данной ячейке
   pts: integer;                  // количество частиц в данной ячейке
   mov: array [0..7] of integer;  // буфер хранящий количество "переезжающих" частиц
   v: myVector2f;                 // вектор перемещения (для сбора статистики, например)
 end;

 myBHydro = object
   cell: array of myBCell;
   h, w: integer;
   function idx_to_crd(i: integer): myVector2i;
   function crd_to_idx(v: myVector2i): integer;
 public
   h_link, v_link: boolean; // на будущее - флаги определяющие замкнута ли карта по горизонтали и вертикали
   density: integer;  // количество частиц в одной единице объёма. теоретически - чем больше, тем лучше, но чем больше - тем медленнее
   min, max: real;
   procedure init(width, height: integer);    // инициализация автомата

   // взаимодействие с моделью

   procedure bottom(x, y: integer; el: real); // задаём отметку дна в указаной точке
   procedure drop(x, y, n: integer);          // добавляет в заданую ячейку n частиц
   function liquid(x, y: integer): real;      // возвращаем толщину слоя жидкости в указаной точке
   function ground(x, y: integer): real;      // возвращаем отметку дна в указаной точке, если кто-то вдруг забыл

   // собсна, сердце модели

   procedure flow(t: real);                   // рассчёт перемещения за отрезок времени t. если t = 1 - переместятся все частицы (если это возможно)
 end;

{ myBHydro }

...

procedure myBHydro.flow(t: real);
var
 a, elc, eln, sum: real;
 c: integer;
 i, j, k, idx: integer;
 v, n: myVector2i;
 M_DIST: array [0..7] of real;     // поправка на расстояние
 SHIFT_X: array [0..7] of integer;  // для определения координат
 SHIFT_Y: array [0..7] of integer;  // тоже
 FINAL: array [0..7] of real;      // для подсчёта итога по текущей ячейке
begin
 // при t больших единицы будет применена рекурсия, поскольку в рамкох одного цикла значения больше единицы не имеют смысла
 if t > 1
   then a := 1
   else a := t;
 // индексация в масках
 // 0 1 2
 // 7 * 3
 // 6 5 4
 // заполняем маски
 for i := 0 to 7 do
   if (i mod 2) = 0
     then M_DIST[i] := 1 / sqrt(2)
     else M_DIST[i] := 1;

 SHIFT_X[0] := -1; SHIFT_Y[0] := -1;
 SHIFT_X[1] := 0;  SHIFT_Y[1] := -1;
 SHIFT_X[2] := 1;  SHIFT_Y[2] := -1;
 SHIFT_X[3] := 1;  SHIFT_Y[3] := 0;
 SHIFT_X[4] := 1;  SHIFT_Y[4] := 1;
 SHIFT_X[5] := 0;  SHIFT_Y[5] := 1;
 SHIFT_X[6] := -1; SHIFT_Y[6] := 1;
 SHIFT_X[7] := -1; SHIFT_Y[7] := 0;

 // заполняем буфер обмена
 for i := 0 to high(cell) do begin
   // если клетка пуста - пропускаем
   if cell[i].pts < 1 then continue;
   // считаем сколько частиц готовы к перемещению
   c := round(a * cell[i].pts);
   if c < 1 then continue;
   // координаты клетки
   v := idx_to_crd(i);
   // инициализируем маску весов
   for j := 0 to 7 do begin
     n.x := v.x + SHIFT_X[j];
     n.y := v.y + SHIFT_Y[j];
     if (n.x < 0) or (n.x >= w) or (n.y < 0) or (n.y >= h)
       then FINAL[j] := 0
       else FINAL[j] := 1;
   end;
   // заполняем ненулевые позиции разницей наша отметка - соседская отметка
   elc := cell[i].btm + cell[i].pts / density;
   sum := 0;
   for j := 0 to 7 do begin
     if FINAL[j] = 0 then continue;
     idx := crd_to_idx(Vector2i(v.x + SHIFT_X[j], v.y + SHIFT_Y[j]));
     eln := cell[idx].btm + cell[idx].pts / density;
     // если там больше, чем здесь - туда не пойдём ни в коем случае
     if eln > elc then begin
       FINAL[j] := 0;
       continue;
     end;
     // сразу введём поправку на расстояние между ячейками и случайный фактор...
     FINAL[j] := (elc - eln) * M_DIST[j];
     sum := sum + FINAL[j];
   end;
   // приводим маску весов к единичному основанию
   for k := 0 to 7 do begin
     try
       FINAL[k] := FINAL[k] / sum;
     except
       FINAL[k] := i / i;
     end;
   end;
   // вписываем значения в буфер обмена и считаем вектор движения
   cell[i].v.x := 0; cell[i].v.y := 0;
   for j := 0 to 7 do begin
     cell[i].mov[j] := round(c * FINAL[j]);
     cell[i].v.x := cell[i].v.x + FINAL[j] * SHIFT_X[j];
     cell[i].v.y := cell[i].v.y + FINAL[j] * SHIFT_Y[j];
   end;
 end;

 // финальный аккорд :: чистая арифметика :: складываем и вычитаем
 for i := 0 to high(cell) do begin
   v := idx_to_crd(i);
   for j := 0 to 7 do begin
     if cell[i].mov[j] < 1 then continue;
     idx := crd_to_idx(Vector2i(v.x + SHIFT_X[j], v.y + SHIFT_Y[j]));
     cell[idx].pts := cell[idx].pts + cell[i].mov[j];
     cell[i].pts := cell[i].pts - cell[i].mov[j];
   end;
 end;

 {max := 0;
 min := maxint;
 for i := 0 to high(cell) do begin
   if cell[i].pts > max then max := cell[i].pts;
   if cell[i].pts < min then min := cell[i].pts;
 end;
 max := max / density;
 min := min / density;      }

 // а вот и рекурсия. на случай если нужно обсчитать более одного кадра за проход
 if t > a then begin
   t := t - a;
   flow(t);
 end;
end;

// заполняем буфер обмена
for i := 0 to high(cell) do begin

FINAL: array [0..7] of real;

for k := 0 to 7 do begin
     try
       FINAL[k] := FINAL[k] / sum;
     except
       FINAL[k] := i / i;
     end;
   end;

for k := 8 downto 1 do
...

 dim := 25;
 mdl.init(dim, dim);
 mdl.density := 1000;
 // initializing water level
 for x := 0 to dim - 1 do
   for y := 0 to dim - 1 do
     mdl.drop(x, y, 4500 + random(1000));

// приводим маску весов к единичному основанию
   for k := 0 to 7 do begin
     writeln(k);
     //try
       FINAL[k] := FINAL[k] / sum;
    { except
       FINAL[k] := i / i;
     end; }
   end;
   Writeln;

Microsoft Windows [Version 6.1.7601]
(c) Корпорация Майкрософт (Microsoft Corp.), 2009. Все права защищены.

C:\Delphi7\Projects\wasser>dbg
model size: 5
0
1
2
3
4
5
6
7

0
1
2
3
4
5
6
7

0
1
2
3
4
5
6
7

0
1
2
3
4
5
6
7

0
Exception EInvalidOp in module dbg.exe at 0004D24E.
Invalid floating point operation.

cell #6
pts = 4578
c = 262
elevation =  4.57800000000000E+0000
eln #0 =  5.11400000000000E+0000
  ignoring this direction
eln #1 =  5.47800000000000E+0000
  ignoring this direction
eln #2 =  5.23000000000000E+0000
  ignoring this direction
eln #3 =  5.39600000000000E+0000
  ignoring this direction
eln #4 =  5.32100000000000E+0000
  ignoring this direction
eln #5 =  5.19300000000000E+0000
  ignoring this direction
eln #6 =  4.99500000000000E+0000
  ignoring this direction
eln #7 =  5.23100000000000E+0000
  ignoring this direction
normalizing mask: 0Exception EInvalidOp in module dbg.exe at 0004D61E.
Invalid floating point operation.

NewVal = (Src[x-1,y-1]  + Src[x, y-1] + Src[x+1, y-1] +
         Src[x-1, y]   +             + Src[x+1, y]   +
         Src[x-1, y+1] + Src[x, y+1] + Src[x+1, y+1]) * 0.25 - Dst[x,y];
Dst[x,y] := NewVal - (NewVal * FadeFactor);