Calculate Real Time Pendulum

const
  RungeKutta = True; // true=RKM4 false=RKM1
  Dim = 4;           // Dimension des Problems = 4

type
  TData = array[0..Dim-1] of Extended;
  
private
    t, x1, y1, x2, y2: Extended;
    u: TData;
    procedure Reset;
    procedure Start;
    procedure Stop;
    procedure Draw;
    procedure CalcTimeStep;
    
const
 u0 : TData = (-pi/4,-pi/2,0,0); // Anfangsbedingung (Winkel & Winkelgeschw.)
 fps = 60;
 N = 256;      // Zwischenschritte pro Frame (Genauigkeit)
 dt = 1/fps/N;
 g = 9.81;     // Gravitation
 l = 1.0;      // Länge in m
 Scale = 70.0; // Pixel pro m
 
//

function f(u: TData; t: Extended): TData;
var
 c, s, a, b: Extended;
begin
 c := cos(u[0]-u[1]);
 s := sin(u[0]-u[1]);
 a := sqr(u[3])*s-2*g/l*cos(u[0]);
 b := -sqr(u[2])*s-g/l*cos(u[1]);

 Result[0] := u[2];
 Result[1] := u[3];
 Result[2] := (a-c*b)/(sqr(c)-2);
 Result[3] := (2*b-c*a)/(sqr(c)-2);
end;

function Mult(factor: Extended; u: TData): TData;
var
 I: Integer;
begin
  for I := 0 to length(u)-1 do
   Result[I] := factor*u[I];
end;

function Add(u, v: TData): TData;
var
 I: Integer;
begin
 for I := 0 to length(u)-1 do
  Result[I] := u[I]+v[I];
end;

function RKM1(u: TData; t: Extended; dt: Extended): TData; // Euler
var
 du: TData;
begin
  du := Mult(dt, f(u, t));
  Result := Add(u, du);
end;

function RKM4(u: TData; t: Extended; dt: Extended): TData; // Runge Kutta 4
var
 v, du: TData;
begin
  du := f(u, t);
  Result := du;
  v := Add(u, Mult(dt/2, du));
  du := f(v, t+dt/2);
  Result := Add(Result, Mult(2, du));
  v := Add(u, Mult(dt/2, du));
  du := f(v, t+dt/2);
  Result := Add(Result, Mult(2, du));
  v := Add(u, Mult(dt, du));
  du := f(v, t+dt);
  Result := Add(Result, du);
  Result := Add(u, Mult(dt/6 ,Result));
end;

procedure TForm1.CalcTimeStep;
var
 I: Integer;
begin
  for I := 1 to N do
  begin
   t := t + dt;
   if RungeKutta then
    u := RKM4(u, t, dt) else
    u := RKM1(u, t, dt);
  end;
end;

// Steuerung & GUI

procedure TForm1.Draw;
 procedure PointLine(Canvas: TCanvas; x, y: Integer);
 begin
   with Canvas do
   begin
    Brush.Color := clBlack;
    LineTo(x,y);
    Rectangle(x-2,y-2,x+3,y+3);
   end;
 end;
var
 x0, y0: Integer;
begin
  with Image1, Image1.Canvas do
  begin
   x0 := Width div 2;
   y0 := Height div 2;
   x1 := l*cos(u[0]);
   y1 := l*sin(u[0]);
   x2 := x1+l*cos(u[1]);
   y2 := y1+l*sin(u[1]);
   Brush.Color := clWhite;
   FillRect(ClipRect);
   MoveTo(x0,y0);
   PointLine(Canvas, x0, y0);
   PointLine(Canvas, Round(x0+x1*Scale), Round(y0+y1*Scale));
   PointLine(Canvas, Round(x0+x2*Scale), Round(y0+y2*Scale));
   Brush.Color := clWhite;
   TextOut(2,2,Format('t = %.1f s',[t]));
   // Anzahl Umdrehungen:
   with StatusBar2 do
   begin
    Panels[0].Text := Format('U. Pendel 1: %d',[Floor(u[0]/pi/2+1/4)]);
    Panels[1].Text := Format('U. Pendel 2: %d',[Floor(u[1]/pi/2+1/4)]);
   end;
  end;
end;

procedure TForm1.Reset;
var
 I: Integer;
begin
  t := 0;
  for I := 0 to Length(u0)-1 do
   u[I] := u0[I];
  Draw;
end;

procedure TForm1.Start;
begin
  Timer1.Interval := Round(dt*N*1000);
  Timer1.Enabled := True;
end;

procedure TForm1.Stop;
begin
  Timer1.Enabled := False;
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
  ControlStyle := ControlStyle + [csOpaque];
  Reset;
  with StatusBar1 do
  begin
   Panels[0].Text := Format('1 m = %.0f Pixel',[Scale]);
   Panels[1].Text := Format('L = %.0f m',[l]);
   Panels[2].Text := Format('g = %.2f m/s²',[g]);
  end;
end;

procedure TForm1.Button1Click(Sender: TObject);
begin
  Start;
end;

procedure TForm1.Button2Click(Sender: TObject);
begin
  Stop;
end;

procedure TForm1.Button3Click(Sender: TObject);
begin
  Reset;
end;

procedure TForm1.Timer1Timer(Sender: TObject);
begin
  CalcTimeStep;
  Draw;
end;