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.