-- test_lsfit.adb  using stuff from  pde_nl33.adb  for test

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Real_Arrays; -- types real, real_vector, real_matrix
use  Real_Arrays;
with Array3d;
use  Array3d;
with Integer_Arrays; -- types integer_vector
use  Integer_Arrays;
with Lsfit;
use  Lsfit;

procedure Test_Lsfit is
  Nx : Integer := 5;      -- number of data in each dimension
  Ny : Integer := 5;
  Nz : Integer := 5;
  Xg : Real_Vector(1..Nx);    -- x grid, does not need to be uniform
  Yg : Real_Vector(1..Ny);    -- y grid
  Zg : Real_Vector(1..Nz);    -- z grid
  Ug : Real_Matrix3d(1..Nx,1..Ny,1..Nz); -- computed solution at xyz grid,
  Ua : Real_matrix3d(1..nx,1..ny,1..nz); -- known includes boundary
  Hx, Hy, Hz : real;
  xmin : Real := -1.0;
  xmax : Real :=  1.0;
  ymin : Real := -1.0;
  ymax : Real :=  1.0;
  zmin : Real := -1.0;
  zmax : Real :=  1.0;
  X, Y, Z, U : Real := 0.0;

  function S(i:Integer; ii:Integer; iii:Integer) return Integer is
  begin
    return (i-2)*(ny-2)*(nz-2)+(ii-2)*(nz-2)+(iii-2)+1; --make single subscript
  end S;

  function ub(x:Real; y:Real; z:Real) return real is -- boundary
  begin
    return x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*x*y*z+5.0*x*y+z+2.0;
  end ub;

  package Real_IO is new Float_IO(Real);
  use Real_IO;

begin
  Put_Line("test_lsfit.adb");

  hx := (xmax-xmin)/real(nx-1);
  Put_Line("nx="&Integer'Image(nx)&", step="&Real'Image(hx));
  for i in 1..nx loop
    xg(i) := xmin+Real(i-1)*hx;
  end loop;
  hy := (ymax-ymin)/real(ny-1);
  Put_Line("ny="&Integer'Image(ny)&", step="&Real'Image(hy));
  for ii in 1..ny loop
    yg(ii) := ymin+Real(ii-1)*hy;
  end loop;
  hz := (zmax-zmin)/real(nz-1);
  Put_Line("nz="&Integer'Image(ny)&", step="& Real'Image(hz));
  for iii in 1..nz loop
    zg(iii) := zmin+Real(iii-1)*hz;
  end loop;

  for i in 1..nx loop
    for ii in 1..ny loop
      for iii in 1..nz loop
        Ua(i,ii,iii) := ub(xg(i),yg(ii),zg(iii));
        Put("xg=");
        Put(xg(i),3,5,0);
        Put(", yg=");
        Put(yg(ii),3,5,0);
        Put(", zg=");
        Put(zg(iii),3,5,0);
        Put(", Ua("&Integer'Image(i)&","&Integer'Image(ii)&","&
            ","&Integer'Image(iii)&")=");
        Put(Ua(i,ii,iii),3,5,0);
        New_Line;
        if i=1 or i=nx or ii=1 or ii=ny or iii=1 or iii=nz then
          ug(i,ii,iii) := ub(xg(i),yg(ii),zg(iii));
        else
          ug(i,ii,iii) := 0.0;
        end if;
      end loop;
    end loop;
  end loop;
  New_Line;

  Fit_Init(4,3); -- max power, dimension

  Fit_Bpoly(Nx, Ny, Nz, Xg, Yg, Zg, Ug);

  -- test fit
  for i in 2..nx-1 loop
     X := Xg(i);
     for ii in 2..ny-1 loop
      Y := Yg(ii);
      for iii in 2..nz-1 loop
        Z := Zg(iii);
        U := Eval_Poly(X, Y, Z);
        Put("expected=");
        Put(Ua(i,ii,iii),8,3,0);
        Put("  fit=");
        Put(U,8,3,0);
        Put("  error="&Real'Image(Ua(i,ii,iii)-U));
        New_Line;
      end loop;
    end loop;
  end loop;
  New_Line;
  Put_Line("test_lsfit.adb finished");
end Test_Lsfit;