-- fourier.adb   for unequally spaced points, FFT may not work,
--               here we do a dFT discrete Fourier Transform to
--               get the Fourier series to approximate y = f(x)
--               from a data file with x,y values not equally
--               spaced. x in increasing order.
--
-- y_approx(xj) = a0/2 + sum_i=j,n( ai*cos(i*2*pi*xj/T)+
--                                  bi*sin(i*2*pi*xj/T) )
--            ai = 2/T * integral_0,T( f(x)*cos(i*2*pi*xj/T)dx
--            bi = 2/T * integral_0,T( f(x)*sin(i*2*pi*xj/T)dx
--
-- for numerical computation, x is scaled to be in range (0,1)
-- xs = (x-xmin)/(xmax-xmin)  T=1
-- then y_approx is computed  y_approx((x-xmin)/(xmax-xmin))
--
-- real data has coefficients aliased at and above nyquist
-- thus n max is (npoints-2)/2
-- a "smoother" fit can be obtained by a Fejer approximation.
-- many coefficients may be near zero.


with ada.text_io;
use  ada.text_io;
with ada.numerics;
use  ada.numerics;
with ada.numerics.Long_elementary_Functions;
use  ada.numerics.Long_elementary_Functions;
with real_arrays;
use real_arrays; -- defines type real


procedure fourier is
  debug : boolean := true;
  pi : long_float := 3.1415926535897932384626433832795028841971;
  x, y, a, b: real_vector(0..1000);
  xmin, xmax, ymin, ymax, xs, xs1, xs2 : real;
  err, avgerr, rmserr, maxerr, sum, sumsq : real := 0.0;
  npts, ntrm : integer;
  tpi : real := 2.0*pi;
  dx, ti: real;

  function min(x:real;y:real) return real is
  begin
     if x<y then return x; end if;
     return y;
  end min;

  function max(x:real;y:real) return real is
  begin
     if x>y then return x; end if;
     return y;
  end max;

  function min(x:integer;y:integer) return integer is
  begin
     if x<y then return x; end if;
     return y;
  end min;

  package real_io is new float_io(real);
  use real_io ;

begin
  -- read data
  npts := 0;
  ntrm := 6;

  begin -- 1
    loop
      begin -- 2
        get(x(npts));
        get(y(npts));
        if debug then
          put_line("npts="&integer'image(npts)&", x="&real'image(x(npts))&
                   ", y="&real'image(y(npts)));
        end if;
        npts := npts+1;
      exception
        when data_error =>
          put_Line(" bad inpUT daTa ");
          skip_Line ;
      end; -- begin 2
    end loop ;
  exception
    when end_error =>
      put_Line ( " end of input data " ) ;
  end; -- begin 1

  if ntrm = -1 then
    ntrm := npts;
  end if;

  ntrm := min(ntrm,(npts-2)/2);
  put_line(integer'image(npts)&" points read, using "&
           integer'image(ntrm)&" terms");

  xmin := x(0);
  xmax := x(0);
  ymin := y(0);
  ymax := y(0);
  for j in 1..npts-1 loop
    xmin := min(x(j),xmin);
    xmax := max(x(j),xmax);
    ymin := min(y(j),ymin);
    ymax := max(y(j),ymax);
  end loop;

  put_line("xmin="&real'image(xmin)&", xmax="&real'image(xmax));
  put_Line("ymin="&real'image(ymin)&", ymax="&real'image(ymax));

  -- integration, simple trapezoidal rule
  put_line("coefficients");
  for i in 0..ntrm-1 loop
    a(i) := 0.0;
    b(i) := 0.0;
    ti := real(i);
    for j in 0..npts-2 loop
      xs1 := (x(j)-xmin)/(xmax-xmin);
      xs2 := (x(j+1)-xmin)/(xmax-xmin);
      dx := xs2-xs1;
      a(i) := a(i) + dx*(y(j)*cos(ti*tpi*xs1)+y(j+1)*
                              cos(ti*tpi*xs2));
      b(i) := b(i) + dx*(y(j)*sin(ti*tpi*xs1)+y(j+1)*
                              sin(ti*tpi*xs2));
    end loop;
    put_line("a("&integer'image(i)&")="&real'image(a(i))&
             ", b("&integer'image(i)&")="&real'image(b(i)));
  end loop;
  put_line(" ");
  -- check fit

  for j in 0..npts-1 loop
    sum := 0.0;
    sumsq := 0.0;
    xs := (x(j)-xmin)/(xmax-xmin);
    for i in 0..ntrm-1 loop
      ti := real(i);
      sum := sum + a(i)*cos(ti*tpi*xs) + b(i)*sin(ti*tpi*xs);
    end loop;
    err := y(j)-sum;
    sumsq := sumsq + err*err;
    avgerr := avgerr + abs(err);
    maxerr := max(abs(err),maxerr);
    if debug then
      put_line("y("&integer'image(j)&")="&real'image(y(j))&
               ", fit="&real'image(sum)&", err"&real'image(err));
    end if;
  end loop;
  rmserr := sqrt(sumsq/real(npts));
  avgerr := avgerr/real(npts);
  put_line("avgerr="&real'image(avgerr)&", rmserr="&real'image(Rmserr)&
           ", maxerr="&real'image(maxerr));
  put_line("fourier.adb finished ");
end Fourier;