-- test_triquad.adb

with Ada.Text_IO;
use  Ada.Text_IO;
with Real_Arrays;
use Real_Arrays;
with Integer_Arrays;
use Integer_Arrays;
with Triquad;
use Triquad;

procedure test_triquad is
  x1 : Real := 1.0;
  y1 : Real := 1.0;
  x2 : Real := 2.5;
  y2 : Real := 3.0;
  x3 : Real := 4.5;
  y3 : Real := 0.5;
  T  : Real_Vector(0..5);
  T2 : Real_Vector(0..5);
  T3 : Real_Vector(0..5);
  X  : Real_Vector(0..35);
  Y  : Real_Vector(0..35);
  W  : Real_Vector(0..35);
  nn : Integer;
  u, a, Area : Real;

  -- functions and procedures used locally
  function f1(x:Real; y:Real) return Real is
  begin
    return x+2.0*y+3.0;
  end f1;

  function f2(x:Real; y:Real) return Real is
  begin
    return x*x + 2.0*x*y + 3.0*y*y + 4.0*x + 5.0*y + 6.0;
  end f2;

  function f3(x:Real; y:Real) return Real is
  begin
    return x*x*x + 2.0*x*x*y + 3.0*x*y*y + 4.0*y*y*y +
           5.0*x*x + 6.0*x*y + 7.0*y*y + 8.0*x + 9.0*y + 10.0;
  end f3;

  function f4(x:Real; y:Real) return Real is
  begin
    return x*x*x*x + 2.0*x*x*x*y + 3.0*x*x*y*y + 4.0*x*y*y*y + 5.0*y*y*y*y +
        6.0*x*x*x + 7.0*x*x*y + 8.0*x*y*y + 9.0*y*y*y +
        10.0*x*x + 11.0*x*y + 12.0*y*y + 13.0*x + 14.0*y + 15.0;
  end f4;

  function f5(x:Real; y:Real) return Real is
  begin
    return x*x*x*x*x + 2.0*x*x*x*x*y + 3.0*x*x*x*y*y + 4.0*x*x*y*y*y +
        5.0*x*y*y*y*y + 6.0*y*y*y*y*y + 7.0*x*x*x*x +
        8.0*x*x*x*y + 9.0*x*x*y*y + 10.0*x*y*y*y + 11.0*y*y*y*y +
        12.0*x*x*x + 13.0*x*x*y + 14.0*x*y*y + 15.0*y*y*y +
        16.0*x*x + 17.0*x*y + 18.0*y*y + 19.0*x + 20.0*y + 21.0;
  end f5;

  function f6(x:Real; y:Real) return Real is
  begin
    return x*x*x*x*x*x + 2.0*x*x*x*x*x*y + 3.0*x*x*x*x*y*y +
        4.0*x*x*x*y*y*y + 5.0*x*x*y*y*y*y + 6.0*x*y*y*y*y*y + 7.0*y*y*y*y*y*y +
        8.0*x*x*x*x*x + 9.0*x*x*x*x*y + 10.0*x*x*x*y*y + 11.0*x*x*y*y*y +
        12.0*x*y*y*y*y + 13.0*y*y*y*y*y + 14.0*x*x*x*x +
        15.0*x*x*x*y + 16.0*x*x*y*y + 17.0*x*y*y*y + 18.0*y*y*y*y +
        19.0*x*x*x + 20.0*x*x*y + 21.0*x*y*y + 22.0*y*y*y +
        23.0*x*x + 24.0*x*y + 25.0*y*y + 26.0*x + 27.0*y + 28.0;
  end f6;

  -- tri_split  split triangles, each one becomes four triangles
  --            similar to the original
  procedure tri_split(ntri:Integer;
                      t1:in out Integer_Vector;
                      t2:in out Integer_Vector;
                      t3:in out Integer_Vector;
                      ncoord:Integer;
                      xc:in out Real_Vector;
                      yc:in out Real_Vector;
                      rtn:out Integer_Vector) is
    j12, j23, j31 : Integer;
    nvert, nntri : Integer;
    i1, i2, i3 : Integer;
    x1, y1, x2, y2, x3, y3 : Real;
  begin
    nvert:=ncoord;
    nntri:=ntri;
    for i in 0..ntri-1 loop
      i1:=t1(i);
      i2:=t2(i);
      i3:=t3(i);
      x1:=xc(i1);
      y1:=yc(i1);
      x2:=xc(i2);
      y2:=yc(i2);
      x3:=xc(i3);
      y3:=yc(i3);

      j12:=nvert;
      nvert := nvert+1;
      xc(j12):=(x1+x2)/2.0;
      yc(j12):=(y1+y2)/2.0;
      -- prevent duplicate vertex
      for k in 0..nvert-2 loop
        if abs(xc(K)-xc(j12))<1.0e-6 and abs(yc(k)-yc(j12))<1.0e-6 then
          j12 := k;
          nvert := nvert-1;
          exit;
        end if;
      end loop;
      j23:=nvert;
      nvert := nvert+1;
      xc(j23):=(x2+x3)/2.0;
      yc(j23):=(y2+y3)/2.0;
      for k in 0..nvert-2 loop
        if abs(xc(K)-xc(j23))<1.0e-6 and abs(yc(k)-yc(j23))<1.0e-6 then
          j23 := k;
          nvert := nvert-1;
          exit;
        end if;
      end loop;
      j31:=nvert;
      nvert := nvert+1;
      xc(j31):=(x3+x1)/2.0;
      yc(j31):=(y3+y1)/2.0;
      for k in 0..nvert-2 loop
        if abs(xc(K)-xc(j31))<1.0e-6 and abs(yc(k)-yc(j31))<1.0e-6 then
          j31 := k;
          nvert := nvert-1;
          exit;
        end if;
      end loop;

      t1(i):=i1; t2(i):=j12; t3(i):=j31;
      t1(nntri):=j12; t2(nntri):=i2; t3(nntri):=j23;
      nntri := nntri+1;
      t1(nntri):=j31; t2(nntri):=j12; t3(nntri):=j23;
      nntri := nntri+1;
      t1(nntri):=j31; t2(nntri):=j23; t3(nntri):=i3;
      nntri := nntri+1;
    end loop;
    rtn(0) := nntri;
    rtn(1) := nvert;
    Put_Line("Tri_Split returns nntri="&Integer'Image(nntri)&
             ", nvert="&Integer'Image(nvert));
  end Tri_Split;

  -- tri_split_test  see how accuracy varies with splitting
  --                 use center points of new triangles
  --                 equally weighted
  procedure tri_split_test(T:Real_Vector) is
    sz : constant :=16386;
    ntri : Integer;                 -- number of triangles
    t1 : Integer_Vector(0..sz);     -- vertex numbers on triangles
    t2 : Integer_Vector(0..sz);     -- vertex numbers on triangles
    t3 : Integer_Vector(0..sz);     -- vertex numbers on triangles
    Ncoord : Integer;               -- number of coordinates
    xc : Real_Vector(0..sz);      -- coordinates in vertex order
    yc : Real_Vector(0..sz);      -- coordinates in vertex order
    x, y, u, a, area : Real;
    rtn : Integer_Vector(0..1);     -- return ntri, ncoord
  begin
    Put_line(" ");
    Put_line("tri_split_test  center point used for numerical quadrature");

    ntri := 1;
    t1(0) := 0;
    t2(0) := 1;
    t3(0) := 2;
    ncoord := 3;
    xc(0) := T(0);
    yc(0) := T(1);
    xc(1) := T(2);
    yc(1) := T(3);
    xc(2) := T(4);
    yc(2) := T(5);
    a := T(0)*T(3) - T(2)*T(1) + T(2)*T(5) - T(4)*T(3) + T(4)*T(1) - T(0)*T(5);
    area := abs(a)/2.0;

    for nsplit in 1..7 loop
      tri_split(ntri, t1, t2, t3, ncoord, xc, yc, rtn);
      ntri := rtn(0);
      ncoord := rtn(1);
      u := 0.0;
      for i in 0..ntri-1 loop
        x := (xc(t1(i))+xc(t2(i))+xc(t3(i)))/3.0;
        y := (yc(t1(i))+yc(t2(i))+yc(t3(i)))/3.0;
        u := u + f6(x,y);
      end loop;
      u := area*u/Real(ntri);
      Put_line(Integer'Image(ntri)&" triangles, integral over T3 of f6:="&
               Real'Image(u));
    end loop;
  end Tri_Split_Test;

  Package Real_IO is new Float_IO(Real);
  use Real_IO;
begin

  Put_line("test_triquad.out");
  Put_line("x1="&Real'Image(x1)&", y1="&Real'Image(y1)&
           ", x2="&Real'Image(x2)&", y2="&Real'Image(y2));
  Put_line("x3="&Real'Image(x3)&", y3="&Real'Image(y3));

  T3(4) := X1; T2(2) := X1; T(0) := x1;
  T3(5) := Y1; T2(3) := Y1; T(1) := y1;
  T3(0) := X2; T2(4) := X2; T(2) := x2;
  T3(1) := Y2; T2(5) := Y2; T(3) := y2;
  T3(2) := X3; T2(0) := X3; T(4) := x3;
  T3(3) := Y3; T2(1) := Y3; T(5) := y3;

  Put_line("T2(0):="&Real'Image(T2(0))&", T2(1):="&Real'Image(T2(1))&
           ", T2(2):="&Real'Image(T2(2))&", T2(3):="&Real'Image(T2(3)));
  Put_line("T2(4):="&Real'Image(T2(4))&", T2(5):="&Real'Image(T2(5)));
  Put_line("T3(0):="&Real'Image(T3(0))&", T3(1):="&Real'Image(T3(1))&
           ", T3(2):="&Real'Image(T3(2))&", T3(3):="&Real'Image(T3(3)));
  Put_line("T3(4):="&Real'Image(T3(4))&", T3(5):="&Real'Image(T3(5)));
  a := x1*y2 - x2*y1 + x2*y3 - x3*y2 + x3*y1 - x1*y3;
  area := abs(a)/2.0;
  Put_line("area:="&Real'Image(area));

  Triquad_msg; -- prints message, runs tests
  nn := 0;
  Put_line("triquad tint(1, T, nn, X, Y, W);");
  tint(1, T, nn, X, Y, W);
  for I in 0..nn-1 loop
     Put_line("x:="&Real'Image(X(i))&", y:="&Real'Image(Y(i))&
              ", w:="&Real'Image(W(i)));
  end loop;
  New_line;

  Put_line("triquad tint(2, T, nn, X, Y, W);");
  tint(2, T, nn, X, Y, W);
  for i in 0..nn-1 loop
     Put_line("x:="&Real'Image(X(i))&", y:="&Real'Image(Y(i))&
              ", w:="&Real'Image(W(i)));
  end Loop;
  New_Line;

  Put_line("triquad tint(3, T, nn, X, Y, W);");
  tint(3, T, nn, X, Y, W);
  for i in 0..nn-1 loop
     Put_line("x:="&Real'Image(X(i))&", y:="&Real'Image(Y(i))&
              ", w:="&Real'Image(W(i)));
  end Loop;
  New_Line;

  Put_line("triquad tint(4, T, nn, X, Y, W);");
  tint(4, T, nn, X, Y, W);
  for i in 0..nn-1 loop
     Put_line("x:="&Real'Image(X(i))&", y:="&Real'Image(Y(i))&
              ", w:="&Real'Image(W(i)));
  end Loop;
  New_Line;

  Put_line("triquad tint(5, T, nn, X, Y, W);");
  tint(5, T, nn, X, Y, W);
  for i in 0..nn-1 loop
     Put_line("x:="&Real'Image(X(i))&", y:="&Real'Image(Y(i))&
              ", w:="&Real'Image(W(i)));
  end Loop;
  New_Line;

  Put_line("triquad tint(6, T, nn, X, Y, W);");
  tint(6, T, nn, X, Y, W);
  for i in 0..nn-1 loop
     Put_line("x:="&Real'Image(X(i))&", y:="&Real'Image(Y(i))&
              ", w:="&Real'Image(W(i)));
  end Loop;
  New_Line;
  New_Line;

  Put_line("f1:=x+2.0*y+3.0");
  for  n in 1..6 loop
    tint(n, T, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f1(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral of f1:="&Real'Image(u));
  end loop;
  New_Line;

  Put_line("f2:=x*x + 2.0*x*y + 3.0*y*y + 4.0*x + 5.0*y + 6.0");
  for  n in 1..6 loop
    tint(n, T, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f2(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral of f2:="&Real'Image(u));
  end loop;
  New_Line;

  Put_line("f3:=x*x*x + 2.0*x*x*y + 3.0*x*y*y + 4.0*y*y*y +");
  Put_line("5.0*x*x + 6.0*x*y + 7.0*y*y + 8.0*x + 9.0*y + 10.0");
  for  n in 1..6 loop
    tint(n, T, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f3(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral of f3:="&Real'Image(u));
  end loop;
  New_Line;

  Put_line("f4:=x*x*x*x + 2.0*x*x*x*y + 3.0*x*x*y*y + 4.0*x*y*y*y + 5.0*y*y*y*y +");
  Put_line("6.0*x*x*x + 7.0*x*x*y + 8.0*x*y*y + 9.0*y*y*y +");
  Put_line("10.0*x*x + 11.0*x*y + 12.0*y*y + 13.0*x + 14.0*y + 15.0");
  for  n in 1..6 loop
    tint(n, T, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f4(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral of f4:="&Real'Image(u));
  end loop;
  New_Line;

  Put_line("f5:=x*x*x*x*x + 2.0*x*x*x*x*y + 3.0*x*x*x*y*y + 4.0*x*x*y*y*y +");
  Put_line("5.0*x*y*y*y*y + 6.0*y*y*y*y*y + 7.0*x*x*x*x +");
  Put_line("8.0*x*x*x*y + 9.0*x*x*y*y + 10.0*x*y*y*y + 11.0*y*y*y*y +");
  Put_line("12.0*x*x*x + 13.0*x*x*y + 14.0*x*y*y + 15.0*y*y*y +");
  Put_line("16.0*x*x + 17.0*x*y + 18.0*y*y + 19.0*x + 20.0*y + 21.0");
  for  n in 1..6 loop
    tint(n, T, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f5(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral of f5:="&Real'Image(u));
  end loop;
  New_Line;

  Put_line("f6:=x*x*x*x*x*x + 2.0*x*x*x*x*x*y + 3.0*x*x*x*x*y*y +");
  Put_line("4.0*x*x*x*y*y*y + 5.0*x*x*y*y*y*y + 6.0*x*y*y*y*y*y + 7.0*y*y*y*y*y*y +");
  Put_line("8.0*x*x*x*x*x + 9.0*x*x*x*x*y + 10.0*x*x*x*y*y + 11.0*x*x*y*y*y +");
  Put_line("12.0*x*y*y*y*y + 13.0*y*y*y*y*y + 14.0*x*x*x*x +");
  Put_line("15.0*x*x*x*y + 16.0*x*x*y*y + 17.0*x*y*y*y + 18.0*y*y*y*y +");
  Put_line("19.0*x*x*x + 20.0*x*x*y + 21.0*x*y*y + 22.0*y*y*y +");
  Put_line("23.0*x*x + 24.0*x*y + 25.0*y*y + 26.0*x + 27.0*y + 28.0");
  for  n in 1..6 loop
    tint(n, T, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f6(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral of f6:="&Real'Image(u));
  end loop;
  New_Line;

  for  n in 1..6 loop
    tint(n, T2, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f6(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral over T2 of f6:="&Real'Image(U));
  end loop;
  New_Line;

  for  n in 1..6 loop
    tint(n, T3, nn, X, Y, W);
    u := 0.0;
    for j in 0..nn-1 loop
      u := u + W(j)*f6(X(j),Y(j));
    end loop;
    Put_line(Integer'Image(nn)&" point integral over T3 of f6:="&Real'Image(u));
  end loop;
  New_Line;

  tri_split_test(T); -- equal area numerical quadrature

  Put_line(" ");
  Put_line("compare accuracy of this 16384 equal area integration");
  Put_line("coordinates with the accuracy of the 36 Gauss Legendre");
  Put_line("integration coordinates in the previous set.");
  Put_line(" ");
  Put_line("P.S. No, you do not get better accuracy using many ");
  Put_line("triangles with a higher order method used in each triangle.");
  Put_line(" ");

end Test_Triquad;