-- fit_sin.adb demonstrate that a lsfit to a specified power,
-- does not give the same coefficients as a truncated approximation,
-- and is a more accurate fit.
--
-- sin(X+Y) = X + Y - (X^3 + 3X^2*Y + 3X*Y^2 + Y^3)/6
--
-- sin(X+Y) = X + Y - (X^3 + 3X^2*Y + 3X*Y^2 + Y^3)/6 +
--            (X^5 + 5X^4*Y + 10X^3*Y^2 + 10X^2*Y^2 + 5X*Y^4 + Y^5)/120
--
-- sin(X+Y+Z) = X + Y + Z -
--              (X^3 + Y^3 + Z^3 + 6X*Y*Z +
--               3X^2*Y + 3X^2*Z + 3X*Y^2 + 3X*Z^2 + 3*Y^2Z + 3*Y*Z^2)/6
--
with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Lsfit;

procedure Fit_Sin is
  subtype Real is Long_Float;
  Hx, Hy, Hz : Real := 0.1;
  xmin : Real := 0.0;
  xmax : Real := 1.0;
  ymin : Real := 0.0;
  ymax : Real := 1.0;
  zmin : Real := 0.0;
  zmax : Real := 1.0;
  X, Y, Z, W3, W5, We : Real := 0.0;
  Err3, Err5, Maxerr : Real := 0.0;
  Last : Boolean := False;
  N2 : Integer := 10;
  N3 : Integer := 7;
  Imax, Jmax, Kmax : Integer := 0;

begin
  Put_Line("fit_sin.adb  demonstrate that a lsfit to a specified power,");
  Put_Line("does not give the same coefficients as a truncated approximation");
  Put_Line("and is a more accurate fit.");

  Put_Line("approximation 2D power 3, sin(X+Y) =");
  Put_Line(" X + Y - (X^3 + 3X^2*Y + 3X*Y^2 + Y^3)/6");
  Put_Line("  1.000 X");
  Put_Line("  1.000 Y");
  Put_Line(" -0.166 X^3");
  Put_Line(" -0.500 X^2*Y");
  Put_Line(" -0.500 X*Y^2");
  Put_Line(" -0.166 Y^3");
  New_Line;
  Put_Line("lsfit terms");
  Lsfit.Fit_Init(3,2); -- power, dimension
  Last := False;
  Maxerr := 0.0;
  for I in 0..N2 loop
    for J in 0..N2 loop
      X := Xmin+Hx*Real(I);
      Y := Ymin+Hy*Real(J);
      We := Sin(X+Y);
      W3 := X+Y-(X*X*X+3.0*X*X*Y+3.0*X*Y*Y+Y*Y*Y)/6.0;
      Err3 := abs(We-W3);
      if Err3 > Maxerr then Maxerr := Err3; end if;
      if Err3>0.1  then
        Put_Line("X="&Real'Image(X)&", Y="&Real'Image(Y));
        Put_Line("      Err3="&Real'Image(Err3));
      end if;
      if I=N2 and J=N2 then Last := True; end if;
      if Last then Put_Line("lsfit coefficients"); end if;
      Lsfit.Fit2_Data(X, Y, We, Last);
    end loop;
  end loop;
  Put_Line("approximation 2D power 3 maxerr="&Real'Image(Maxerr));
  New_Line;
  -- check lsfit 3,2
  Maxerr := 0.0;
  for I in 0..N2 loop
    for J in 0..N2 loop
      X := Xmin+Hx*Real(I);
      Y := Ymin+Hy*Real(J);
      We := Sin(X+Y);
      W3 := Lsfit.Eval2_Poly(X, Y);
      Err3 := abs(We-W3);
      if Err3 > Maxerr then Maxerr := Err3; end if;
      if Err3>0.1  then
        Put_Line("X="&Real'Image(X)&", Y="&Real'Image(Y));
        Put_Line(" lsfit 2,3  Err3="&Real'Image(Err3));
      end if;
    end loop;
  end loop;
  Put_Line("lsfit 3,2 maxerr="&Real'Image(Maxerr));
  New_Line;

  Put_Line("approximation 2D power 5, sin(X+Y) =");
  Put_Line(" X + Y - (X^3 + 3X^2*Y + 3X*Y^2 + Y^3)/6 +");
  Put_Line(" (X^5 + 5X^4*Y + 10X^3*Y^2 + 10X^2*Y^2 + 5X*Y^4 + Y^5)/120");
  New_Line;
  Put_Line("lsfit terms");
  Lsfit.Fit_Init(5,2);
  Last := False;
  for I in 0..N2 loop
    for J in 0..N2 loop
      X := Xmin+Hx*Real(I);
      Y := Ymin+Hy*Real(J);
      We := Sin(X+Y);
      W5 := X+Y-(X*X*X+3.0*X*X*Y+3.0*X*Y*Y+Y*Y*Y)/6.0+
                (X*X*X*X*X+5.0*X*X*X*X*Y+10.0*X*X*X*Y*Y+
                 10.0*X*X*Y*Y*Y+5.0*X*Y*Y*Y*Y+Y*Y*Y*Y*Y)/120.0;
      Err5 := abs(We-W5);
      if Err5 > Maxerr then Maxerr := Err5; end if;
      if Err5>0.01  then
        Put_Line("X="&Real'Image(X)&", Y="&Real'Image(Y));
        Put_Line("      Err5="&Real'Image(Err5));
      end if;
      if I=N2 and J=N2 then Last := True; end if;
      if Last then Put_Line("lsfit coefficients"); end if;
      Lsfit.Fit2_Data(X, Y, We, Last);
    end loop;
  end loop;
  Put_Line("approximation 2D power 5 maxerr="&Real'Image(Maxerr));
  New_Line;
  -- check lsfit 5,2
  Maxerr := 0.0;
  for I in 0..N2 loop
    for J in 0..N2 loop
      X := Xmin+Hx*Real(I);
      Y := Ymin+Hy*Real(J);
      We := Sin(X+Y);
      W5 := Lsfit.Eval2_Poly(X, Y);
      Err5 := abs(We-W5);
      if Err5 > Maxerr then Maxerr := Err5; end if;
      if Err5>0.001  then
        Put_Line("X="&Real'Image(X)&", Y="&Real'Image(Y));
        Put_Line(" lsfit 5,2  Err5="&Real'Image(Err5));
      end if;
    end loop;
  end loop;
  Put_Line("lsfit 5,2 maxerr="&Real'Image(Maxerr));
  New_Line;

  Put_Line("approximation 3D power 3, sin(X+Y+Z) =");
  Put_Line(" X + Y + Z -");
  Put_Line(" (X^3 + Y^3 + Z^3 + 6X*Y*Z +");
  Put_Line("  3X^2*Y + 3X^2*Z + 3X*Y^2 + 3X*Z^2 + 3*Y^2Z + 3*Y*Z^2)/6");
  New_Line;
  Put_Line("lsfit terms");
  Lsfit.Fit_Init(3,3);
  Last := False;
  Maxerr := 0.0;
  for I in 0..N3 loop
    for J in 0..N3 loop
      for k in 0..N3 loop
        X := Xmin+Hx*Real(I);
        Y := Ymin+Hy*Real(J);
        Z := Zmin+Hz*Real(K);
        We := Sin(X+Y+Z);
        W3 := X+Y+Z-(X*X*X+Y*Y*Y+Z*Z*Z+6.0*X*Y*Z+
                     3.0*X*X*Y+3.0*X*X*Z+3.0*X*Y*Y+
                     3.0*X*Z*Z+3.0*Y*Y*Z+3.0*Y*Z*Z)/6.0;
        Err3 := abs(We-W3);
        if Err3 > Maxerr then Maxerr := Err3; end if;
        if Err3>0.2  then
          Put_Line("X="&Real'Image(X)&", Y="&Real'Image(Y)&
                   ", Z="&Real'Image(Z));
          Put_Line("                      Err3="&Real'Image(Err3));
        end if;
        if I=N3 and J=N3 and K=N3 then Last := True; end if;
        if Last then Put_Line("lsfit coefficients"); end if;
        Lsfit.Fit3_Data(X, Y, Z, We, Last);
      end loop;
    end loop;
  end loop;
  Put_Line("approximation 3D power 3 maxerr="&Real'Image(Maxerr));
  New_Line;
  -- check lsfit 3,3
  Maxerr := 0.0;
  for I in 0..N3 loop
    for J in 0..N3 loop
      for k in 0..N3 loop
        X := Xmin+Hx*Real(I);
        Y := Ymin+Hy*Real(J);
        Z := Zmin+Hz*Real(K);
        We := Sin(X+Y+Z);
        W3 := Lsfit.Eval3_Poly(X, Y, Z);
        Err3 := abs(We-W3);
        if Err3 > Maxerr then Maxerr := Err3; end if;
        if Err3>0.3  then
          Put_Line("X="&Real'Image(X)&", Y="&Real'Image(Y));
          Put_Line(" lsfit 3,3  Err3="&Real'Image(Err3));
        end if;
      end loop;
    end loop;
  end loop;
  Put_Line("lsfit 3,3 maxerr="&Real'Image(Maxerr));
  New_Line;
  Put_Line("fit_sin.adb ending");
end Fit_Sin;