unit Unit1;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, Forms, Controls, Graphics, Dialogs, ExtCtrls, StdCtrls,
  ComCtrls, vecMath, Math;

type

  { TForm1 }

  TForm1 = class(TForm)
    Button1: TButton;
    Button2: TButton;
    Button3: TButton;
    Button4: TButton;
    CheckBox1: TCheckBox;
    Image1: TImage;
    Image2: TImage;
    Image3: TImage;
    Label1: TLabel;
    ProgressBar1: TProgressBar;
    RadioGroup1: TRadioGroup;
    ZeichenTimer: TTimer;
    procedure FormActivate(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure ZeichenTimerTimer(Sender: TObject);
    procedure CheckBox1Click(Sender: TObject);
    procedure CheckBox2Click(Sender: TObject);
    procedure Image3MouseDown(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure Image3MouseUp(Sender: TObject; Button: TMouseButton;
      Shift: TShiftState; X, Y: Integer);
    procedure Image3MouseMove(Sender: TObject; Shift: TShiftState; X,
      Y: Integer);
    procedure FormResize(Sender: TObject);
    procedure RadioGroup1Click(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure Button3Click(Sender: TObject);
    procedure Button4Click(Sender: TObject);
  private

  public
    Dim,Koerp:   Integer;
    Kanten:      array of TPoint;
    Ecken:       array of TMatrix;
    ProjektionL,
    ProjektionR,
    Rotation,
    Beschl,Rot:  TMatrix;
    Zoom:        extended;

    procedure Drehen;
    procedure BeschlBelegen;
    procedure RotationInitialisieren;
    procedure ProjektionenBelegen;
    procedure KoerperErzeugen;
    procedure Zeichnen;
    procedure KantenZeichnen(Pte: array of TPoint; PZ: array of extended; drawZ: boolean; drawTo: TCanvas);
    procedure drawColoredLine(wo: TCanvas; X1,Y1: Integer; Z1: extended; X2,Y2: Integer; Z2: Extended);
    function farbVerlauf(x: Extended): TColor;
    procedure SetZoom;
  end;

var
  Form1: TForm1;

implementation

{$R *.lfm}

uses
  unit2;

procedure TForm1.FormActivate(Sender: TObject);
begin
 Dim:=3;
 Koerp:=-1;
 RadioGroup1.OnClick(Sender);
end;

procedure TForm1.FormDestroy(Sender: TObject);
var i: integer;
begin
 for i:=0 to length(Ecken)-1 do
  Ecken[i].destroy;
 ProjektionL.destroy;
 ProjektionR.destroy;
 Rotation.destroy;
 Rot.destroy;
 Beschl.destroy;
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
 setlength(Ecken,0);
 ProjektionL:=TMatrix.create;
 ProjektionR:=TMatrix.create;
 Rotation:=TMatrix.create;
 Rot:=TMatrix.create;
 Beschl:=TMatrix.create;
end;

procedure TForm1.ZeichenTimerTimer(Sender: TObject);
var i: integer;
begin
 Zeichnen;
 for i:=0 to ZeichenTimer.Interval div 10 do
  begin
   Drehen;
   BeschlBelegen;
  end;
end;

procedure TForm1.ProjektionenBelegen;
var z,s: integer;
    tmp: TMatrix;
const alpha = 2*pi/360 * 10;
begin
 ProjektionL.destroy;
 ProjektionL:=TMatrix.create;
 ProjektionL.xdim:=4;
 ProjektionL.ydim:=3;
 for z:=0 to 2 do
  for s:=0 to 3 do
   ProjektionL.a[z,s]:=byte(z=s);

 ProjektionR.destroy;
 ProjektionR:=TMatrix.create;
 ProjektionR.xdim:=4;
 ProjektionR.ydim:=4;
 ProjektionR.assign0;
 ProjektionR.a[0,0]:=cos(alpha);
 ProjektionR.a[0,2]:=sin(alpha);
 ProjektionR.a[3,0]:=-sin(alpha);
 ProjektionR.a[3,2]:=cos(alpha);
 ProjektionR.a[1,1]:=1;
 ProjektionR.a[2,3]:=1;

 ProjektionR._times(ProjektionL);

 tmp:=TMatrix.create;
 tmp.xdim:=dim;
 tmp.ydim:=4;
 for s:=0 to dim-1 do
  for z:=0 to 3 do
   tmp.a[z,s]:=byte(z=s);

 ProjektionR.times_(tmp);
 ProjektionL.times_(tmp);
end;

procedure TForm1.KoerperErzeugen;
function spiegleKante(Knr,Enr: integer; N: TMatrix): integer;
var Es:     array[boolean] of integer;
    i:      integer;
    M,O:    TMatrix;
    bol,nk: boolean;
    tmp:    extended;
begin
  M:=TMatrix.create;
  O:=TMatrix.create;
  Es[false]:=Kanten[Knr].x;
  Es[true]:=Kanten[Knr].y;
  for bol:=false to true do begin
    M.assign(Ecken[Enr]);
    M.minus(Ecken[Es[bol]]);
    M.transp;
    M.times_(N);
    tmp:=M.a[0,0];
    M.assign(N);
    M.times(2*tmp);
    M.plus(Ecken[Es[bol]]);
    tmp:=1;
    for i:=0 to length(Ecken)-1 do begin
      O.assign(M);
      O.minus(Ecken[i]);
      tmp:=min(tmp,O.len);
    end;
    if tmp>1e-6 then begin
      Es[bol]:=length(Ecken);
      setlength(Ecken,length(Ecken)+1);
      Ecken[Es[bol]]:=TMatrix.create;
      Ecken[Es[bol]].assign(M);
    end
    else begin
      for i:=0 to length(Ecken)-1 do begin
        O.assign(M);
        O.minus(Ecken[i]);
        if O.len <=1.1e-6 then begin
          Es[bol]:=i;
          break;
        end;
      end;
    end;
  end;
  result:=-1;
  nk:=true;
  for i:=0 to length(kanten)-1 do
    nk:=nk and ((Kanten[i].x<>Es[false]) or (Kanten[i].y<>Es[true]))
           and ((Kanten[i].y<>Es[false]) or (Kanten[i].x<>Es[true]));
  if nk then begin
    setlength(Kanten,length(Kanten)+1);
    Kanten[length(Kanten)-1].x:=Es[false];
    Kanten[length(Kanten)-1].y:=Es[true];
    result:=length(Kanten)-1;
  end
  else
  for i:=0 to length(kanten)-1 do begin
    if (Kanten[i].x=Es[false]) and (Kanten[i].y=Es[true]) then result:=i;
    if (Kanten[i].y=Es[false]) and (Kanten[i].x=Es[true]) then result:=i;
  end;
  M.destroy;
  O.destroy;
end;

var ral:         array of array[boolean] of TMatrix;

procedure rotierePunkt(Quelle: TMatrix; Achse: integer; umNull: boolean; var Ziel: TMatrix);
var tang,otang: TMatrix;
    tmp:        extended;
begin
  Ziel.assign(Quelle);
  if ral[Achse,true].len<1e-6 then exit;
  if not umNull then Ziel.minus(ral[Achse,false]);
  tang:=TMatrix.create;
  otang:=TMatrix.create;
  tang.assign(Ziel);
  Ziel.transp;
  Ziel.times_(ral[Achse,true]);
  tmp:=Ziel.a[0,0];
  Ziel.assign(ral[Achse,true]);
  Ziel.times(tmp/sqr(ral[Achse,true].len));
  tang.minus(Ziel);
  if tang.len<1e-6 then begin
    if not umNull then Ziel.plus(ral[Achse,false]);
    tang.destroy;
    otang.destroy;
    exit;
  end;
  otang.assign(ral[Achse,true]);
  otang.cross(tang);
  otang.times(tang.len/otang.len*sin(ral[Achse,true].len));
  tang.times(cos(ral[Achse,true].len));
  Ziel.plus(tang);
  Ziel.plus(otang);
  if not umNull then Ziel.plus(ral[Achse,false]);
  tang.destroy;
  otang.destroy;
end;

procedure rotiereKante(Kante,Achse: integer);
var M,N:     TMatrix;
    Es:      array[boolean] of integer;
    bol,gef: boolean;
    i:       integer;
begin
  M:=TMatrix.create;
  N:=TMatrix.create;
  Es[false]:=Kanten[Kante].x;
  Es[true]:=Kanten[Kante].y;
  for bol:=false to true do begin
    rotierePunkt(Ecken[Es[bol]],Achse,false,M);
    gef:=false;
    for i:=0 to length(Ecken)-1 do begin
      N.assign(M);
      N.minus(Ecken[i]);
      if N.len<1e-6 then begin
        gef:=true;
        Es[bol]:=i;
        break;
      end;
    end;
    if not gef then begin
      setlength(Ecken,length(Ecken)+1);
      Ecken[length(Ecken)-1]:=TMatrix.create;
      Ecken[length(Ecken)-1].assign(M);
      Es[bol]:=length(Ecken)-1;
    end;
  end;
  gef:=false;
  for i:=0 to length(Kanten)-1 do
    gef:=gef or ((Kanten[i].x=Es[false]) and (Kanten[i].y=Es[true])) or
                ((Kanten[i].y=Es[false]) and (Kanten[i].x=Es[true]));
  if not gef then begin
    setlength(Kanten,length(Kanten)+1);
    Kanten[length(Kanten)-1].x:=Es[false];
    Kanten[length(Kanten)-1].y:=Es[true];
  end;

  M.destroy;
  N.destroy;
end;

var skl:             array of TPoint;
    i,j,k:           integer;
    M,N:             TMatrix;
    nsa,bol:         boolean;
    b,d,f,t,phi,tmp: extended;
begin
 case Koerp of
  0:
   begin
    if dim<>3 then exit;
//    b:=3/22*(7-sqrt(5));               // kurze Seite      (lange Seite=1)
    d:=3*sqrt((5-sqrt(5))/20);         // kurze Diagonale
    f:=1/11*sqrt((470+156*sqrt(5))/5); // lange Diagonale

    Setlength(Ecken,16);
    for i:=0 to length(Ecken)-1 do begin
      Ecken[i]:=TMatrix.create;
      Ecken[i].xdim:=1;
      Ecken[i].ydim:=dim;
    end;
    Setlength(Kanten,35);
    for i:=0 to length(Kanten)-1 do
      Kanten[i].x:=-1;
    Ecken[0].a[0,0]:=0;
    Ecken[0].a[1,0]:=0;
    Ecken[0].a[2,0]:=sqrt(1-sqr(d*cos(pi/5)/sin(2*pi/5)));  // die berechnete Höhe
    for i:=0 to 4 do begin
      Ecken[i+1].a[0,0]:=cos(2*pi/5*i)*d*cos(pi/5)/sin(2*pi/5);
      Ecken[i+1].a[1,0]:=sin(2*pi/5*i)*d*cos(pi/5)/sin(2*pi/5);
      Ecken[i+1].a[2,0]:=0;
    end;
    M:=TMatrix.create;
    M.xdim:=1;
    M.ydim:=dim;
    N:=TMatrix.create;
    N.xdim:=1;
    N.ydim:=dim;
    for i:=0 to 4 do begin
      Ecken[i+6].assign(Ecken[i+1]);
      Ecken[i+6].plus(Ecken[(i+1) mod 5 +1]);
      Ecken[i+6].times(-1/2);
      Ecken[i+6].plus(Ecken[0]);
      Ecken[i+6].times(-f/Ecken[i+6].len);
      Ecken[i+6].plus(Ecken[0]);
      Ecken[i+11].assign(Ecken[0]);
      M.assign(Ecken[i+6]);
      Ecken[i+11].times(0.5);
      M.times(0.5);
      Ecken[i+11].plus(M);
      M.assign(Ecken[i+1]);
      M.minus(Ecken[0]);
      N.assign(Ecken[(i+1) mod 5+1]);
      N.minus(Ecken[0]);
      M.cross(N);
      M.times(3/M.len);
      Ecken[i+11].plus(M);
    end;
    for i:=0 to 4 do begin
      Kanten[i].x:=0;
      Kanten[i].y:=i+1;
      Kanten[i+5].x:=i+1;
      Kanten[i+5].y:=i+6;
      Kanten[i+10].x:=i+6;
      Kanten[i+10].y:=(i+1) mod 5 +1;
      Kanten[i+15].x:=i+11;
      Kanten[i+15].y:=0;
      Kanten[i+20].x:=i+11;
      Kanten[i+20].y:=i+1;
      Kanten[i+25].x:=i+11;
      Kanten[i+25].y:=i+6;
      Kanten[i+30].x:=i+11;
      Kanten[i+30].y:=(i+1) mod 5 +1;
    end;

    setlength(skl,5);
    for i:=0 to 4 do begin
      skl[i].x:=i+5;
      skl[i].y:=(i+4) mod 5+10;
    end;

    repeat
      nsa:=false;
      i:=length(skl)-1;
      while i>=0 do begin
        M.assign(Ecken[Kanten[skl[i].x].x]);
        M.minus(Ecken[Kanten[skl[i].x].y]);
        N.assign(Ecken[Kanten[skl[i].y].x]);
        N.minus(Ecken[Kanten[skl[i].y].y]);
        N.cross(M);
        N.times(1/N.len);
        j:=length(skl)-1;
        while j>=0 do begin
          setlength(skl,length(skl)+1);
          skl[length(skl)-1].x:=spiegleKante(skl[j].x,Kanten[skl[i].x].x,N);
          skl[length(skl)-1].y:=spiegleKante(skl[j].y,Kanten[skl[i].x].x,N);
          bol:=false;
          for k:=0 to length(skl)-2 do
            bol:=bol or ((skl[length(skl)-1].x=skl[k].x) and
                         (skl[length(skl)-1].y=skl[k].y)) or
                        ((skl[length(skl)-1].y=skl[k].x) and
                         (skl[length(skl)-1].x=skl[k].y));
          if bol then setlength(skl,length(skl)-1)
          else nsa:=true;
          dec(j);
        end;
        dec(i);
      end;
    until not nsa;

    for i:=0 to length(skl)-1 do begin
      M.assign(Ecken[Kanten[skl[i].x].x]);
      M.minus(Ecken[Kanten[skl[i].x].y]);
      N.assign(Ecken[Kanten[skl[i].y].x]);
      N.minus(Ecken[Kanten[skl[i].y].y]);
      N.cross(M);
      N.times(1/N.len);
      j:=length(Kanten)-1;
      while j>=0 do begin
        spiegleKante(j,Kanten[skl[i].x].x,N);
        dec(j);
      end;
    end;
    M.destroy;
    N.destroy;
   end;
  1:
   begin
    if dim<>3 then exit;
    phi:=(sqrt(5)+1)/2;
    t:=1/12*(power(44+12*phi*(9+sqrt(81*phi-15)),1/3)+power(44+12*phi*(9-sqrt(81*phi-15)),1/3)-4);
    b:=2*(1-2*sqr(t))/(1+2*t);         // kurze Seite      (lange Seite=1)
    d:=2*(1-2*sqr(t));                 // Diagonale
    f:=sqrt(b*(b+d));                  // (ganz) kurze Diagonale

    Setlength(Ecken,21);
    for i:=0 to length(Ecken)-1 do begin
      Ecken[i]:=TMatrix.create;
      Ecken[i].xdim:=1;
      Ecken[i].ydim:=dim;
    end;
    Setlength(Kanten,45);
    for i:=0 to length(Kanten)-1 do
      Kanten[i].x:=-1;
    Ecken[0].a[0,0]:=0;
    Ecken[0].a[1,0]:=0;
    Ecken[0].a[2,0]:=sqrt(1-sqr(d*cos(pi/5)/sin(2*pi/5)));  // die berechnete Höhe
    for i:=0 to 4 do begin
      Ecken[i+1].a[0,0]:=cos(2*pi/5*i)*d*cos(pi/5)/sin(2*pi/5);
      Ecken[i+1].a[1,0]:=sin(2*pi/5*i)*d*cos(pi/5)/sin(2*pi/5);
      Ecken[i+1].a[2,0]:=0;
    end;
    M:=TMatrix.create;
    M.xdim:=1;
    M.ydim:=dim;
    N:=TMatrix.create;
    N.xdim:=1;
    N.ydim:=dim;
    for i:=0 to 4 do begin
      Ecken[i+6].assign(Ecken[i+1]);
      Ecken[i+6].plus(Ecken[(i+1) mod 5 +1]);
      Ecken[i+6].times(1/2);
      Ecken[i+6].minus(Ecken[0]);
      Ecken[i+6].times(sqrt(sqr(b)-sqr((d-b)/2))/Ecken[i+6].len);
      N.assign(Ecken[(i+1) mod 5 +1]);
      N.minus(Ecken[i+1]);
      N.times((d-b)/2/N.len);
      Ecken[i+6].plus(N);
      Ecken[i+6].plus(Ecken[i+1]);
      N.times(b/N.len);
      Ecken[i+11].assign(Ecken[i+6]);
      Ecken[i+11].plus(N);
      M.assign(Ecken[i+6]);
      M.plus(Ecken[i+11]);
      M.times(0.3);
      Ecken[i+16].assign(Ecken[0]);
      Ecken[i+16].times(0.7);
      Ecken[i+16].plus(M);
      M.assign(Ecken[i+1]);
      M.minus(Ecken[0]);
      N.assign(Ecken[(i+1) mod 5 +1]);
      N.minus(Ecken[0]);
      M.cross(N);
      M.times(3/M.len);
      Ecken[i+16].Plus(M);
    end;
    for i:=0 to 4 do begin
      Kanten[i].x:=0;
      Kanten[i].y:=i+1;
      Kanten[i+5].x:=i+1;
      Kanten[i+5].y:=i+6;
      Kanten[i+10].x:=i+6;
      Kanten[i+10].y:=i+11;
      Kanten[i+15].x:=i+11;
      Kanten[i+15].y:=(i+1) mod 5 +1;
      Kanten[i+20].x:=0;
      Kanten[i+20].y:=i+16;
      Kanten[i+25].x:=i+1;
      Kanten[i+25].y:=i+16;
      Kanten[i+30].x:=i+6;
      Kanten[i+30].y:=i+16;
      Kanten[i+35].x:=i+11;
      Kanten[i+35].y:=i+16;
      Kanten[i+40].x:=(i+1) mod 5 +6;
      Kanten[i+40].y:=i+16;
    end;

    setlength(ral,2);
    for i:=0 to 3 do
      ral[i div 2,odd(i)]:=TMatrix.create;
    ral[0,false].assign(Ecken[0]);             // die 5er-Symmetrieachse
    ral[0,true].assign(Ecken[0]);
    ral[0,true].times(2*pi/5/ral[0,true].len);

    ral[1,false].assign(Ecken[6]);             // die 3er-Symmetrieachse
    tmp:=arccos((sqrt(3)/6*f)/sqrt(sqr(b)-sqr(f/2))); // Kippwinkel der Rotationsachse gegenüber Flächennormale
    ral[1,true].assign(Ecken[1]);
    ral[1,true].plus(Ecken[11]);
    ral[1,true].times(-1/2);
    ral[1,true].plus(Ecken[6]);
    ral[1,true].times(1/ral[1,true].len);
    N.assign(Ecken[1]);
    N.minus(Ecken[11]);
    N.cross(ral[1,true]);
    N.times(1/N.len);
    ral[1,true].times(sin(tmp));
    N.times(cos(tmp));
    ral[1,true].plus(N);
    ral[1,true].times(2*pi/3/ral[1,true].len);

    repeat
      nsa:=false;
      i:=length(ral)-1;
      while i>=0 do begin
        j:=length(ral)-1;
        while j>=0 do begin
          setlength(ral,length(ral)+1);
          for bol:=false to true do begin
            ral[length(ral)-1,bol]:=TMatrix.create;
            rotierePunkt(ral[j,bol],i,bol,ral[length(ral)-1,bol]);
          end;
          bol:=false;
          for k:=0 to length(ral)-2 do begin
            M.assign(ral[k,false]);
            M.minus(ral[length(ral)-1,false]);
            if M.len<1e-6 then begin
              M.assign(ral[k,true]);
              M.minus(ral[length(ral)-1,true]);
              if M.len<1e-6 then begin
                bol:=true;
                break;
              end;
            end;
          end;
          if bol then begin
            ral[length(ral)-1,false].destroy;
            ral[length(ral)-1,true].destroy;
            setlength(ral,length(ral)-1);
          end
          else nsa:=true;
          dec(j);
        end;
        dec(i);
      end;
    until not nsa;

    for i:=0 to length(ral)-1 do begin
      j:=length(Kanten)-1;
      while j>=0 do begin
        rotiereKante(j,i);
        dec(j);
      end;
    end;

    for i:=0 to 2*length(ral)-1 do
      ral[i div 2,odd(i)].destroy;
    setlength(ral,0);


    M.destroy;
    N.destroy;
   end;
 end{of Case};
 M:=TMatrix.create;
 M.assign0;
 for i:=0 to length(Ecken)-1 do
   M.plus(Ecken[i]);
 M.times(-1/length(Ecken));
 for i:=0 to length(Ecken)-1 do
   Ecken[i].plus(M);
 b:=0;
 for i:=0 to length(Ecken)-1 do
   b:=max(b,Ecken[i].len);
 for i:=0 to length(Ecken)-1 do
   Ecken[i].times(1/2/b);
 M.Destroy;

 for i:=0 to length(Kanten)-1 do
   if Kanten[i].x>Kanten[i].y then begin
     j:=Kanten[i].x;
     Kanten[i].x:=Kanten[i].y;
     Kanten[i].y:=j;
   end;

 for i:=0 to length(Kanten)-2 do begin
   j:=i;
   for k:=i+1 to length(Kanten)-1 do
     if  (Kanten[k].x<Kanten[j].x) or
        ((Kanten[k].x=Kanten[j].x) and
         (Kanten[k].y<Kanten[j].y)) then
       j:=k;
   if j>i then begin
     k:=Kanten[j].x;
     Kanten[j].x:=Kanten[i].x;
     Kanten[i].x:=k;
     k:=Kanten[j].y;
     Kanten[j].y:=Kanten[i].y;
     Kanten[i].y:=k;
   end;
 end;

 SetZoom;
end;

procedure TForm1.CheckBox1Click(Sender: TObject);
begin
 ZeichenTimer.Enabled:=Checkbox1.Checked;
 Image3.Visible:=not Checkbox1.Checked;
end;

procedure TForm1.CheckBox2Click(Sender: TObject);
begin
 ZeichenTimer.Interval:=30;
end;

procedure TForm1.KantenZeichnen(Pte: array of TPoint; PZ: array of extended; drawZ: boolean; drawTo: TCanvas);
var i: integer;
begin
 if drawZ then
  begin
   for i:=0 to length(Kanten)-1 do
    drawColoredLine(drawto,Pte[Kanten[i].X].X,Pte[Kanten[i].X].Y,PZ[Kanten[i].X],
                           Pte[Kanten[i].Y].X,Pte[Kanten[i].Y].Y,PZ[Kanten[i].Y]);
  end
   else
  begin
   for i:=0 to length(Kanten)-1 do
    begin
     drawto.MoveTo(Pte[Kanten[i].X].X,Pte[Kanten[i].X].Y);
     drawto.LineTo(Pte[Kanten[i].Y].X,Pte[Kanten[i].Y].Y);
    end;
  end;
end;

procedure TForm1.drawColoredLine(wo: TCanvas; X1,Y1: Integer; Z1: extended; X2,Y2: Integer; Z2: Extended);
var i:         integer;
    a,b,az,bz: extended;
begin
 if abs(X2-X1)>abs(Y1-Y2) then
  begin
   a:=(Y1-Y2)/(X1-X2);
   b:=Y1-X1*a;
   az:=(Z1-Z2)/(X1-X2);
   bz:=Z1-X1*az;
   if X1<X2 then
    begin
     for I:=X1 to X2 do
      wo.Pixels[I,round(I*a + b)]:=farbVerlauf(I*az + bz);
    end
     else
    begin
     for I:=X2 to X1 do
      wo.Pixels[I,round(I*a + b)]:=farbVerlauf(I*az + bz);
    end;
  end
   else
  begin
   if Y1=Y2 then exit;
   a:=(X1-X2)/(Y1-Y2);
   b:=X1-Y1*a;
   az:=(Z1-Z2)/(Y1-Y2);
   bz:=Z1-Y1*az;
   if Y1<Y2 then
    begin
     for I:=Y1 to Y2 do
      wo.Pixels[round(I*a + b),I]:=farbVerlauf(I*az + bz);
    end
     else
    begin
     for I:=Y2 to Y1 do
      wo.Pixels[round(I*a + b),I]:=farbVerlauf(I*az + bz);
    end
  end;
end;

function TForm1.farbVerlauf(x: Extended): TColor;
begin
 x:=5*x;
 if x<1 then
  begin
   result:=min($FF,max($0,round($FF*x))) shl 16;
   exit;
  end;
 x:=x-1;
 if x<1 then
  begin
   result:=(min($FF,max($0,round($FF*(1-x)))) shl 16) or
           (min($FF,max($0,round($FF*x))) shl 8);
   exit;
  end;
 x:=x-1;
 if x<1 then
  begin
   result:=$00FF00 or
           min($FF,max($0,round($FF*x)));
   exit;
  end;
 x:=x-1;
 if x<1 then
  begin
   result:=$0000FF or
           (min($FF,max($0,round($FF*(1-x)))) shl 8);
   exit;
  end;
 x:=x-1;
 result:=$0000FF or
         (min($FF,max($0,round($FF*x))) shl 16);
end;

procedure TForm1.RotationInitialisieren;
var z,s: integer;
begin
 Rot.xdim:=dim;
 Rot.ydim:=dim;
 repeat
  for z:=0 to dim-1 do
   for s:=0 to dim-1 do
    Rot.a[z,s]:=0.01*random;
  Rotation.assignId(dim);
  Rotation.plus(Rot);
 until Rotation.det>0;
end;

procedure TForm1.BeschlBelegen;
var z,s: integer;
    M:   TMatrix;
begin
 M:=TMatrix.create;
 M.xdim:=dim;
 M.ydim:=dim;
 repeat
  for z:=0 to dim-1 do
   for s:=0 to dim-1 do
    M.a[z,s]:=0.02*random;
  Beschl.assignId(dim);
  Beschl.plus(M);
 until Beschl.det>0;
 Beschl.orthogonalize(true);
 M.destroy;
end;

procedure TForm1.Drehen;
begin
 Rotation.assignId(dim);
 Rotation.plus(Rot);
 Rotation.orthogonalize(true);
 Rot.times_(Beschl);
 ProjektionL.times_(Rotation);
 ProjektionR.times_(Rotation);
end;

procedure TForm1.Image3MouseDown(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
 Image3.Tag:=floor(Y*dim*(dim-1)/2/Image3.Height)+
             X*((dim*(dim-1)) div 2);
end;

procedure TForm1.Image3MouseUp(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
 Image3.Tag:=-1;
end;

procedure TForm1.Image3MouseMove(Sender: TObject; Shift: TShiftState; X,
  Y: Integer);
var M:     TMatrix;
    a1,a2: integer;
const alpha = 2*pi/360 * 1;
begin
 if not(ssLeft in Shift) then exit;
 Y:=Image3.Tag mod ((dim*(dim-1)) div 2);
 X:=X - (Image3.Tag div ((dim*(dim-1)) div 2));
 Image3.Tag:=Image3.Tag + X*((dim*(dim-1)) div 2);
 M:=TMatrix.create;
 M.assignId(dim);
 a1:=0;
 while Y>=0 do
  begin
   inc(a1);
   Y:=Y-dim+a1;
  end;
 a2:=Y+dim;
 dec(a1);
 M.a[a1,a1]:=cos(alpha*X);
 M.a[a1,a2]:=sin(alpha*X);
 M.a[a2,a1]:=-sin(alpha*X);
 M.a[a2,a2]:=cos(alpha*X);
 ProjektionR.times_(M);
 ProjektionL.times_(M);
 M.destroy;
 Zeichnen;
end;

procedure TForm1.Zeichnen;
var Pkt:   TMatrix;
    Pkte:  array[boolean] of array of TPoint;
    PZ:    array of extended;
    i:     integer;
    b:     boolean;
    mi,ma: extended;
begin
 if length(Ecken)=0 then
  exit;
 setlength(Pkte[false],length(Ecken));
 setlength(Pkte[true],length(Ecken));
 setlength(PZ,length(Ecken));
 for i:=0 to length(Pkte[false])-1 do
  for b:=false to true do
   begin
    Pkt:=Ecken[i].duplicate;
    if B then Pkt._times(ProjektionR)
         else Pkt._times(ProjektionL);
    PZ[i]:=Pkt.a[2,0];
    Pkte[b,i]:=Pkt.asPoint(Zoom);
    Pkte[b,i].x:=Pkte[b,i].x + Image1.Width div 2;
    Pkte[b,i].y:=Pkte[b,i].y + Image1.Height div 2;
    Pkt.destroy;
   end;
 mi:=PZ[0];
 ma:=PZ[0];
 for i:=0 to length(PZ)-1 do
  begin
   mi:=min(mi,PZ[i]);
   ma:=max(ma,PZ[i]);
  end;
 ma:=max(ma,mi+0.00001);
 for i:=0 to length(PZ)-1 do
  PZ[i]:=(PZ[i]-mi)/(ma-mi);
 Image1.Canvas.Rectangle(-10,-10,Image1.Width+10,Image1.Height+10);
 KantenZeichnen(Pkte[true],PZ,false,Image1.Canvas);
 Image2.Canvas.Rectangle(-10,-10,Image2.Width+10,Image2.Height+10);
 KantenZeichnen(Pkte[false],PZ,false,Image2.Canvas);
end;

procedure TForm1.FormResize(Sender: TObject);
var L,T: integer;
begin
 L:=Image1.Left;
 T:=Image1.Top;
 Image1.Free;
 Image2.Free;
 Image1:=TImage.Create(Form1);
 Image2:=TImage.Create(Form1);
 Image1.Parent:=Form1;
 Image2.Parent:=Form1;
 Image1.Left:=L;
 Image1.Top:=T;
 Image2.Top:=T;
 Image1.Width:=min(Form1.ClientHeight-Image1.Top,(Form1.ClientWidth-Image1.Left) div 2);
 Image1.Height:=Image1.Width;
 Image2.Width:=Image1.Width;
 Image2.Height:=Image2.Width;
 Image2.Left:=Image1.Left+Image1.Width;
 Progressbar1.Width:=Form1.Clientwidth-Progressbar1.Left;
 SetZoom;
 Zeichnen;
end;

procedure TForm1.SetZoom;
begin
 case Koerp of
  0: Zoom:=min(Image1.Width,Image1.Height)*0.9;
  1: Zoom:=min(Image1.Width,Image1.Height)*0.9;
 end{of Case};
end;

procedure TForm1.RadioGroup1Click(Sender: TObject);
var i: integer;
begin
 if Koerp = Radiogroup1.ItemIndex then exit;
 Koerp:=RadioGroup1.ItemIndex;

 for i:=0 to length(Ecken)-1 do
  Ecken[i].destroy;
 ProjektionenBelegen;
 RotationInitialisieren;
 BeschlBelegen;
 KoerperErzeugen;
 Image3.Canvas.Brush.Color:=Form1.Color;
 Image3.Canvas.Rectangle(-10,-10,Image3.Width+10,Image3.Height+10);
 For i:=1 to (dim*(dim-1)) div 2 do
  begin
   Image3.Canvas.MoveTo(10,round((I-0.5)*(Image3.Height*2/(dim*(dim-1)))));
   Image3.Canvas.LineTo(Image3.Width-10,round((I-0.5)*(Image3.Height*2/(dim*(dim-1)))));
  end;
// Drehen;
// BeschlBelegen;
 Zeichnen;
end;

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

procedure TForm1.Button2Click(Sender: TObject);
var M:         TMatrix;
    I,J,K,L,N: Longint;
    B:         Boolean;
    Memo1:     TMemo;
    s:         String;
    x:         extended;
    fcs:       array of Longint;
begin
  M:=TMatrix.create;
  M.xdim:=length(Kanten)-3;
  for I:=0 to length(Kanten)-1 do begin
    B:=Kanten[I].x=Kanten[I].y;
    for J:=0 to I-1 do
      B:=B or ((Kanten[I].X=Kanten[J].X) and (Kanten[I].Y=Kanten[J].Y))
           or ((Kanten[I].X=Kanten[J].Y) and (Kanten[I].Y=Kanten[J].X));
    if B then M.xdim:=M.xdim-1;
  end;
  M.ydim:=3*length(Ecken)-9;
  if M.xdim<>M.ydim then begin
    if MessageDlg(
      'Die Anzahl der Randbedingungen ('
       +inttostr(M.xdim)+
      ') unterscheidet sich von der Anzahl der Freiheitsgrade ('
      +inttostr(M.ydim)+
      ')! Genauere Erklärung gewünscht?',mterror,[mbYes,mbNo],0) = mrYes then begin
      Memo1:=TMemo.Create(Form1);
      Memo1.Parent:=Form1;
      Memo1.Left:=Image1.Left;
      Memo1.Top:=Image1.Top;
      Memo1.Width:=Form1.ClientWidth-Memo1.Left;
      Memo1.Height:=Form1.ClientHeight-Memo1.Top;
      Memo1.Lines.Clear;
      x:=0;
      for I:=0 to length(Ecken)-1 do
        x:=max(x,Ecken[I].len);
      s:='';
      for I:=0 to length(Kanten)-1 do
        s:=s+' '+inttostr(Kanten[I].x)+'-'+inttostr(Kanten[I].y);
      delete(s,1,1);
      Memo1.Lines.Add(s);
      Memo1.Lines.Add('');
      s:='';
      for I:=0 to length(Ecken)-1 do
        if Ecken[I].len>=0.98*x then
          s:=s+' '+inttostr(I);
      delete(s,1,1);
      Memo1.Lines.Add(s);
      Memo1.Lines.Add('');
      s:='';
      for I:=0 to length(Kanten)-1 do
        if (Ecken[Kanten[I].x].len>=0.98*x) or
           (Ecken[Kanten[I].y].len>=0.98*x) then
          s:=s+' '+inttostr(Kanten[I].x)+'-'+inttostr(Kanten[I].y);
      delete(s,1,1);
      Memo1.Lines.Add(s);
      Memo1.Lines.Add('');
      s:='';
      for i:=0 to length(Ecken)-1 do
        if Ecken[I].len>=0.98*x then begin
          s:=s+' '+inttostr(i)+':';
          for j:=0 to length(Kanten)-1 do
            if (Kanten[j].x=i) or (Kanten[j].y=i) then
              s:=s+inttostr(Kanten[j].x+Kanten[j].y-i)+',';
          delete(s,length(s),1);
        end;
      delete(s,1,1);
      Memo1.Lines.Add(s);
      Application.ProcessMessages;
      MessageDlg('Memo schließen.',mtInformation,[mbOk],0);
      Memo1.Free;
    end;
    exit;
  end;
  M.assign0;

  case Radiogroup1.ItemIndex of
    0: begin
      // Ecken 1,6,11
      setlength(fcs,3);
      fcs[0]:=1;
      fcs[1]:=6;
      fcs[2]:=11;
    end;
    1: begin
      // Ecken 0,1,16
      setlength(fcs,3);
      fcs[0]:=0;
      fcs[1]:=1;
      fcs[2]:=16;
    end;
  end{of case};

  j:=0;
  for i:=0 to length(Kanten)-1 do begin
    K:=Kanten[i].x;
    L:=Kanten[i].y;
    for N:=length(fcs)-1 downto 0 do begin
      if fcs[N]=K then K:=-1;
      if fcs[N]=L then L:=-1;
      if K>=fcs[N] then dec(K);
      if L>=fcs[N] then dec(L);
    end;
    if K<>-1 then
      for N:=0 to 2 do
        M.a[j,3*K+N]:=
          M.a[j,3*K+N]
           + Ecken[Kanten[i].x].a[N,0]
           - Ecken[Kanten[i].y].a[N,0];
    if L<>-1 then
      for N:=0 to 2 do
        M.a[j,3*L+N]:=
          M.a[j,3*L+N]
           + Ecken[Kanten[i].y].a[N,0]
           - Ecken[Kanten[i].x].a[N,0];
    if (K<>-1) or (L<>-1) then inc(j);
  end;

  M.Progressbar:=Progressbar1;
  Messagedlg('det M = '+floattostr(M.det(dmGauss)),mtInformation,[mbOk],0);
  M.Free;
end;

procedure TForm1.Button3Click(Sender: TObject);
var I,J,K,L: Longint;
    x:       extended;
begin
  x:=0;
  for I:=0 to length(Ecken)-1 do
    x:=max(x,Ecken[I].len);
  I:=length(Ecken)-1;
  while I>=0 do begin
    if Ecken[I].len>=0.98*x then begin
      for J:=I+1 to length(Ecken)-1 do
        Ecken[J-1].assign(Ecken[J]);
      Ecken[length(Ecken)-1].Free;
      setlength(Ecken,length(Ecken)-1);
      J:=length(Kanten)-1;
      L:=0;
      while J>=0 do begin
        if (Kanten[J].x=I) or (Kanten[J].y=I) then begin
          inc(L);
          for K:=J+1 to length(Kanten)-1 do
            Kanten[K-1]:=Kanten[K];
          setlength(Kanten,length(Kanten)-1);
        end
        else begin
          if Kanten[J].x>I then dec(Kanten[J].x);
          if Kanten[J].y>I then dec(Kanten[J].y);
        end;
        dec(J);
      end;
      if L<>4 then Form1.Caption:=Form1.Caption+' '+inttostr(L);
    end;
    dec(I);
  end;
  x:=0;
  for I:=0 to length(Ecken)-1 do
    x:=max(x,Ecken[I].len);
  for i:=0 to length(Ecken)-1 do
    Ecken[i].times(1/2/x);
  Zeichnen;
end;

procedure TForm1.Button4Click(Sender: TObject);
const mindim = 4;
      maxdim = 9;
var M:       TMatrix;
    dm:      TDetMethode;
    I,J,K,L: Longint;
    ts,errs: array[mindim..maxdim,TDetMethode] of extended;
    t1,t2:   extended;
    s:       string;
begin
  M:=TMatrix.create;
  t2:=0;
  for I:=mindim to maxdim do begin
    for dm:=low(dm) to high(dm) do begin
      ts[I,dm]:=0;
      errs[I,dm]:=0;
    end;
    for J:=1 to 50 do begin
      M.assignId(I);
      for K:=0 to I-1 do
        for L:=0 to I-1 do
          if random*50<=J then
            M.a[K,L]:=M.a[K,L]+random-0.5;
      for dm:=low(dm) to high(dm) do begin
        ts[I,dm]:=ts[I,dm]-now;
        t1:=M.det(dm);
        ts[I,dm]:=ts[I,dm]+now;
        if dm=low(dm) then t2:=t1;
        errs[I,dm]:=sqr(t1-t2);
      end;
    end;
  end;
  M.Free;
  s:='';
  for i:=mindim to maxdim do begin
    s:=s+inttostr(i)+': ';
    for dm:=low(dm) to high(dm) do
      s:=s+floattostr(errs[i,dm])+' '+floattostr(ts[i,dm]*24*60*60)+';';
    delete(s,length(s),1);
    s:=s+#13;
  end;
  delete(s,length(s),1);
  Messagedlg(s,mtinformation,[mbOk],0);
end;

end.