Unit Rekenen;
{
Automated Balancing of Chemical Equations
=========================================
This software has been designed and it is
CopyLefted by Han de Bruijn:

                          (===)
                         @-O^O-@
                          #/_\#
                           ###

Convert Data structure to Solution vector
-----------------------------------------
}
INTERFACE

Uses Breuken, Algemeen, Numeriek;

function Berekenen(data : blok; hoe : integer; var uit : vektor) : byte;

IMPLEMENTATION

procedure afmetingen(data : blok; var ii,jj : integer);
{
  Matrix size derived from data
}
var
  L,k,i,j : integer;
begin
  L := Length(data);
  if debug then
  for k := 0 to L-1 do
  begin
     Write(data[k].e:5);
     Write(data[k].p:5);
     Write(data[k].L:5);
     Write(data[k].r:5);
     Write(data[k].m:5);
     Write(data[k].g:5);
     Writeln;
  end;
  ii := 0; jj := 0;
  for k := 0 to L-1 do
  begin
    i := data[k].r;
    if i > ii then ii := i;
    j := data[k].m;
    if j > jj then jj := j;
  end;
  ii := ii + 1;
  jj := jj + 1;
end;

function Berekenen(data : blok; hoe : integer; var uit : vektor) : byte;
{
  All calculations
}
var
  C : matrix;
  H : matrix; d : vektor;
  A : Bmatrix; b : Bvektor;
  F : Dmatrix; v : Dvektor;
  i,j,k,L,max : integer;
  ii,jj,r,rest,Det : integer;
  P : double;
  OK : boolean;
  fout : byte;
begin
{ Convert data structure to
  [#atoms,#molecules] matrix }
  afmetingen(data,ii,jj);
  SetLength(C,ii,jj);
  for i := 0 to ii-1 do
  for j := 0 to jj-1 do
    C[i,j] := 0;
  L := Length(data);
  for k := 0 to L-1 do
  begin
    i := data[k].r;
    j := data[k].m;
    C[i,j] := C[i,j] + data[k].g;
  end;
  if debug then
  begin
    Writeln;
    Writeln(ii,' x ',jj,' matrix');
    for i := 0 to ii-1 do
    begin
      for j := 0 to jj-1 do
        Write(C[i,j]:4);
      Writeln;
    end;
    Writeln;
  end;

  ii := ii - 1;
  jj := jj - 1;
  SetLength(H,jj,jj);
  SetLength(d,jj);
  if (hoe = 2) and (jj = ii+1) then
  begin
  { Cramer's rule eventually }
    for i := 0 to jj-1 do
    for j := 0 to jj-1 do
      H[i,j] := C[i,j];
    for i := 0 to jj-1 do
      d[i] := - C[i,jj];
  end else
  begin
  { Least Squares Method }
    for i := 0 to jj-1 do
    for j := i to jj-1 do
    begin
      H[i,j] := 0;
      for k := 0 to ii do
        H[i,j] := H[i,j] + C[k,i]*C[k,j];
      H[j,i] := H[i,j];
    end;
    for i := 0 to jj-1 do
    begin
      d[i] := 0;
      for k := 0 to ii do
         d[i] := d[i] - C[k,i]*C[k,jj];
    end;
    if debug then
      Afdrukken(H,d);
  end;
  SetLength(b,jj);

{ Fractions arithmetic }
  if hoe = 0 then
  begin
    SetLength(A,jj,jj);
    for i := 0 to jj-1 do
    for j := 0 to jj-1 do
      A[i,j] := maak(H[i,j],1);
    for i := 0 to jj-1 do
      b[i] := maak(d[i],1);
    Bsolve(A,b);
  end;

{ Floating point method }
  if hoe = 1 then
  begin
    SetLength(F,jj,jj);
    SetLength(v,jj);
    for i := 0 to jj-1 do
    for j := 0 to jj-1 do
      F[i,j] := H[i,j];
    for i := 0 to jj-1 do
      v[i] := d[i];
    Dsolve(F,v,P);
    for j := 0 to jj-1 do
    begin
      b[j] := maak(Round(v[j]*P),Round(P));
      b[j] := vereenvoudig(b[j]);
    end;
  end;
  
{ Cramer's Rule Solution
  extremely inefficient! }
  if hoe = 2 then
  begin
    DirekteInverse(H,Det);
    Vermenigvuldig(H,d);
  { Afdrukken(H,d); }
    for j := 0 to jj-1 do
    begin
      b[j] := maak(d[j],Det);
      b[j] := vereenvoudig(b[j]);
    end;
  end;

  if debug then
  begin
    for j := 0 to jj-1 do
      Write(' ',afdruk(b[j]));
    Writeln;
  end;
  SetLength(b,jj+1);
  b[jj] := maak(1,1);
{ Establish Solution }
  SetLength(uit,jj+1);
  max := 0;
  for j := 0 to jj do
    if b[j].noemer > max
    then max := b[j].noemer;
  for j := 0 to jj do
    b[j] := maal(b[j],maak(max,1));
  for j := 0 to jj do
    uit[j] := b[j].teller;
  if debug then  
  begin
    for j := 0 to jj do
      Write(' ',uit[j]);
    Writeln;
  end;

{ Balancing Check }
  rest := 0;
  for i := 0 to ii do
  begin
    r := 0;
  { For each Element }
    for j := 0 to jj do
      r := r + C[i,j] * uit[j];
    rest := rest + abs(r);
  end;
  if rest > 0 then OK := false
              else OK := true;
{ Improper Balancing }
  fout := 0;
  if not OK then
  begin
  { Multiple solutions }
    for j := 0 to jj do
      if b[j].noemer = 0 then fout := 1;
  { Zero solution only }
    if fout = 0 then fout := 2;
    for j := 0 to jj do
      uit[j] := 0;
  end;
  Berekenen := fout;
end;

END.