File : n_balls.adb

with Reps, Reps.Ops, Messages, Lin, Lin.IO.Vinci;
use Reps, Reps.Ops, Messages, Lin, Lin.IO.Vinci;
with Lass;

procedure N_Balls is  -- computing volume of N independent balls, also writes .ine file

  NumEqu: constant Integer := 12;
  NumVar: constant Integer := 3;

  NBalls: Integer;
  V: Rep;
  A: RMatPtr := new RMat(1 .. NumEqu, 0 .. NumVar);

begin
  Verbose := 2;
  SetZero(A.all);

  A( 1,0) := Half;  A( 1,1) :=-Half;  A( 1,2) := S3/Two;    A( 1,3) := Zero;
  A( 2,0) := Half;  A( 2,1) := Half;  A( 2,2) := S3/Two;    A( 2,3) := Zero;
  A( 3,0) := Half;  A( 3,1) := Zero;  A( 3,2) := S3/Three;  A( 3,3) := S6/Three;
  A( 4,0) := Half;  A( 4,1) := Zero;  A( 4,2) := S3/Three;  A( 4,3) :=-S6/Three;
  A( 5,0) := Half;  A( 5,1) :=-One;   A( 5,2) := Zero;      A( 5,3) := Zero;
  A( 6,0) := Half;  A( 6,1) := One;   A( 6,2) := Zero;      A( 6,3) := Zero;
  A( 7,0) := Half;  A( 7,1) :=-Half;  A( 7,2) :=-S3/Six;    A( 7,3) := S6/Three;
  A( 8,0) := Half;  A( 8,1) :=-Half;  A( 8,2) :=-S3/Six;    A( 8,3) :=-S6/Three;
  A( 9,0) := Half;  A( 9,1) := Half;  A( 9,2) :=-S3/Six;    A( 9,3) := S6/Three;
  A(10,0) := Half;  A(10,1) := Half;  A(10,2) :=-S3/Six;    A(10,3) :=-S6/Three;
  A(11,0) := Half;  A(11,1) :=-Half;  A(11,2) :=-S3/Two;    A(11,3) := Zero;
  A(12,0) := Half;  A(12,1) := Half;  A(12,2) :=-S3/Two;    A(12,3) := Zero;

  NBalls := PromptInteger("Number of balls: ");

  declare
    NumEqu2: constant Integer := NBalls*NumEqu;
    NumVar2: constant Integer := NBalls*NumVar;

    package Lasserre is new Lass (NumEqu => NumEqu2, NumVar => NumVar2);
    use Lasserre;

    B: RMatPtr := new RMat(1 .. NumEqu2, 0 .. NumVar2);

  begin
    MinDoTrans :=    5;    -- min dimension for translations
    MinVolSave :=    4;    -- min dimension for considering saved volume
    MaxVolSave :=   99;    -- max dimension for using saved volume
    MinVolConv :=    4;    -- min dimension for converting saved volume

    SetZero(B.all);
    for N in 0 .. NBalls-1 loop
      for I in 1 .. NumEqu loop
        B(N*NumEqu+I,0) := A(I,0);
        for J in 1 .. NumVar loop
          B(N*NumEqu+I,N*NumVar+J) := A(I,J);
        end loop;
      end loop;
    end loop;
    Free(A);

    WriteVinci("n_balls.ine",B);

    Vol(B,V);
    Show("Volume of balls is ",V,17,4);
  end;

end N_Balls;