function OK = harmanal_test(varargin)
%HARMANAL_TEST   test for harmanal
%
%See also: harmanal_test

D.T       = [ .5  1  2]; % day
D.phi     = [10  20 30]; % deg
D.A       = [ 1   2  3]; % any unit
D.w       = 2*pi./D.T;   % rad/day
D.A_noise = [2];         % any unit

OPT.eps = min(D.A)./1e12;

%% no noise, with linear transform
   slope     = 2;

   T.t     = 0:(1/24):max(D.T);
   T.h     = 0.*T.t+4;
   for j=1:length(D.T)
   T.h     = T.h + D.A(j).*cos(D.w(j).*T.t - deg2rad(D.phi(j)));
   end
   FIT = harmanal(T.t,T.h,'omega',D.w,'screenoutput',0,'transformFun',@(x) slope*x,'transformFunInv',@(x) x./slope);
   
   FIG(1) = figure;
   plot  (T.t,T.h,'.-')
   hold on
   xlabel('day')
   xlabel('signal')
   plot  (T.t,FIT.hfit,'go')
   ylim([0 12])
   title('linear transform')
   
   OK = 1;
   for i=1:length(D.A)
      if abs(D.A(i).*slope-FIT.hamplitudes(i)) > OPT.eps
      OK = 0;return
      end
      if abs(D.phi(i)-rad2deg(FIT.hphases(i))) > OPT.eps
      OK = 0;return
      end
   end

%% no noise, with poly transform
   poly    = 2.0; % higher values more weight
   poly    = 1/2; % lower  values more weight

   T.t     = 0:(1/24):max(D.T);
   T.h     = 0.*T.t+4;
   for j=1:length(D.T)
   T.h     = T.h + D.A(j).*cos(D.w(j).*T.t - deg2rad(D.phi(j)));
   end
   FIT = harmanal(T.t,T.h,'omega',D.w,'screenoutput',0,'transformFun',@(x) x.^poly,'transformFunInv',@(x) x.^(1/poly));
   
   FIG(2) = figure;
   plot  (T.t,T.h,'.-')
   hold on
   xlabel('day')
   xlabel('signal')
   plot  (T.t,FIT.hfit,'go')
   ylim([0 12])
   title('poly transform')
   
%% no noise, with log transform
%  lower  values more weight

   T.t     = 0:(1/24):max(D.T);
   T.h     = 0.*T.t+4;
   for j=1:length(D.T)
   T.h     = T.h + D.A(j).*cos(D.w(j).*T.t - deg2rad(D.phi(j)));
   end
   try % log10
   FIT = harmanal(T.t,T.h,'omega',D.w,'screenoutput',0,'transformFun',@(x) log10(x),'transformFunInv',@(x) 10.^(x));
   catch
   OK = 0;
   end

   try % log
   FIT = harmanal([nan nan],[nan nan],'omega',nan,'screenoutput',0,'transformFun',@(x)   log(x),'transformFunInv',@(x) 10.^(x));
   OK = 0;
   catch
   end
   
   FIG(3) = figure;
   plot  (T.t,T.h,'.-')
   hold on
   xlabel('day')
   xlabel('signal')
   plot  (T.t,FIT.hfit,'go')
   ylim([0 12])
   title('log transform')

%% with noise

   T.t     = 0:(1/24):max(D.T);
   T.h     = 0.*T.t+4;
   for j=1:length(D.T)
   T.h     = T.h + D.A(j).*cos(D.w(j).*T.t - deg2rad(D.phi(j)));
   end
   D.noise = D.A_noise*(-1+2*rand(size(T.t)));
   T.h     = T.h + D.noise;
   FIT  = harmanal(T.t,T.h,'omega',D.w,'screenoutput',0,'timeseries',D.noise);
   FIT2 = harmanal(T.t,T.h,'omega',D.w,'screenoutput',0);
   
   FIG(4) = figure;
   plot  (T.t,T.h,'r.-','DisplayName','clean data + noise')
   hold on
   plot  (T.t,T.h - D.noise,'b.-','DisplayName','clean data')
   xlabel('day')
   xlabel('signal')
   plot  (T.t,FIT.hfit ,'go' ,'DisplayName','fit of (clean data + noise) when feeding noise timeseries')
   plot  (T.t,FIT2.hfit,'g.-','DisplayName','fit of (clean data + noise) to get rid of noise')
   ylim([0 12])
   legend show
   title('noise timeseries fitted too')

   for i=1:length(D.A)
      if abs(D.A(i)-FIT.hamplitudes(i)) > OPT.eps & ...
            (FIT.c-1) > OPT.eps % noise is fitted to its exact amplitude
      OK = 0;break
      end
      if abs(D.phi(i)-rad2deg(FIT.hphases(i))) > OPT.eps
      OK = 0;break
      end
   end

%% close
   
   if nargin>0
   pausedisp
   end

   try
      close(FIG)
   end

