unit Unit1;

{$mode objfpc}{$H+}

interface

uses

  Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs, ExtCtrls;

type

  { TForm1 }

  TForm1 = class(TForm)

    Timer1: TTimer;

    procedure FormCreate(Sender: TObject);

    procedure FormPaint(Sender: TObject);

    procedure Timer1Timer(Sender: TObject);

  private

    { private declarations }

  public

    { public declarations }

  end;

  TBall = record

    x, y, vx, vy, ax, ay, r, m: double;

  end;

const

  ballsCount = 5;

  g = 100;

var

  Form1: TForm1;

  balls: array[0..ballsCount] of TBall;

  ringX: double = 500;

  ringY: double = 500;

implementation

uses Math;

{$R *.lfm}

function RingNormalize(a, ringSize: double): double;

begin

  if ringSize = 0 then Result := a else

    Result := a - ringSize*Floor(a/ringSize);

end;

function RingDelta(a, b, ringSize: double): double;

begin

  if ringSize = 0 then Result := a - b else begin

    Result := RingNormalize(a - b, ringSize);

    if Result > ringSize/2 then Result := Result - ringSize;

  end;

end;

{ TForm1 }

procedure TForm1.FormCreate(Sender: TObject);

var

  i: integer;

  r: double;

begin

  ringX := ClientWidth;

  ringY := ClientHeight;

  Randomize;

  for i := 0 to ballsCount-1 do begin

    r := 10 + 30*Random;

    balls[i].x := ringX * Random;

    balls[i].y := ringY * Random;

    balls[i].vx := 0;

    balls[i].vy := 0;

    balls[i].r := r;

    balls[i].m := r*r*r;

  end;

end;

procedure TForm1.FormPaint(Sender: TObject);

var

  i, j, k: integer;

  x, y, r: double;

begin

  for i := 0 to ballsCount-1 do begin

    r := balls[i].r;

    for j := -1 to 1 do for k := -1 to 1 do begin

      x := balls[i].x + ringX * j;

      y := balls[i].y + ringY * k;

      Canvas.Ellipse(round(x-r), round(y-r), round(x+r), round(y+r));

    end;

  end;

end;

procedure TForm1.Timer1Timer(Sender: TObject);

var

  t, i, j: integer;

  dx, dy, d, nx, ny, ax, ay, dt, dvx, dvy, dv

  , v0, v1: double;

begin

  dt := Timer1.Interval/1000/1000;

  for t := 1 to 1000 do begin

  for i := 0 to ballsCount-1 do begin

    balls[i].ax := 0;

    balls[i].ay := 0;

  end;

  for i := 0 to ballsCount-1 do begin

    for j := i+1 to ballsCount-1 do begin

      dx := RingDelta(balls[j].x, balls[i].x, ringX);

      dy := RingDelta(balls[j].y, balls[i].y, ringY);

      d := sqrt(dx*dx + dy*dy);

      nx := dx / d;

      ny := dy / d;

      if d > balls[i].r + balls[j].r then begin

        ax := g*balls[i].m*balls[j].m*nx/(d*d);

        ay := g*balls[i].m*balls[j].m*ny/(d*d);

        balls[i].ax := balls[i].ax + ax;

        balls[i].ay := balls[i].ay + ay;

        balls[j].ax := balls[j].ax - ax;

        balls[j].ay := balls[j].ay - ay;

      end else begin

        dvx := balls[j].vx - balls[i].vx;

        dvy := balls[j].vy - balls[i].vy;

        dv := dvx*nx + dvy*ny;

        if dv<0 then begin

          v0 :=  2*dv*balls[j].m/(balls[i].m + balls[j].m);

          v1 := -2*dv*balls[i].m/(balls[i].m + balls[j].m);

          balls[i].vx := balls[i].vx + v0*nx;

          balls[i].vy := balls[i].vy + v0*ny;

          balls[j].vx := balls[j].vx + v1*nx;

          balls[j].vy := balls[j].vy + v1*ny;

        end;

      end;

    end;

  end;

  for i := 0 to ballsCount-1 do begin

    balls[i].vx := balls[i].vx + balls[i].ax*dt/balls[i].m;

    balls[i].vy := balls[i].vy + balls[i].ay*dt/balls[i].m;

    balls[i].x := RingNormalize(balls[i].x + balls[i].vx*dt, ringX);

    balls[i].y := RingNormalize(balls[i].y + balls[i].vy*dt, ringY);

  end;

  end;

  Refresh;

end;

end.