function [OKall,OKlegend] = t_tide_test
%T_TIDE_TEST  test for t_tide, incl. non-canonical variable dt option
%
% Testst 324 combinations of :
% * 2 solutions methods
% * 3 synthesis options
% * 3 start time options
% * 2 types of time series (odd and even length)
% * 3 types of dt (constant, nearly constant (eps) & non constant)
% * 3 parameters: amplitudes, phases, fitted/synthetic timeseries
%
%See also: t_tide

err        = 'cboot'; %'linear'% 'cboot' does not allow any dt variations
rnd        = 1e-9; % dt variations are neglected if smaller than eps('single')

tstart     = datenum(1968,1,10); % [datenum]
nt         = 12;                 % [#] %% if 12 for 24hr series, canonical t_tide 1.3 crashes
tend       = tstart+1.333;       % [datenum] % test for short period not aligned with S2, to test centraltime
T          = 12/24;              % [day] % set to S2
dtfac      = 2;                  % [-] day dt by dt/dtfac

OK         = [];
doplot     = 0;

% test 2 x 3 x 3 = 18 combinations
% NOTE when k=1, i.e.start_time=[]: centraltime is unknown, only
% internally, so phases cannot be checked.

values.lsq        = {'normal','direct'};
values.synthesis  = {-1,0,1};
values.start_time = {[],0,tstart}; % canonical t_predic has test for jdmid>1
for i=1:length(values.lsq)
for j=1:length(values.synthesis)
for k=1:length(values.start_time)

lsq        = values.lsq{i};
synthesis  = values.synthesis{j};
start_time = values.start_time{k};
OKlegend{i,j,k} = ['lsq=',lsq,'    synthesis=',num2str(synthesis),'    start\_time=',num2str(start_time)];

%% odd
   to  = linspace(tstart,tend,nt+(rem(nt,2)+1));
   to  = to + (rand(size(to))-.5).*rnd;
   dto = (to(2) - to(1)).*24;
   to(1)=tstart;
   to(to<tstart) = tstart;to(to>tend) = tend;to(end)=tend;to = sort(to);
   zo  = sin(2*pi.*to./T);
   disp('-----odd----------t----------------------')
   [Do1,Eo1] = t_tide(zo,'int',diff(to).*24,'output','none','error',err,'Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   disp('-----odd----------dt---------------------')
   [Do2,Eo2] = t_tide(zo,'int',dto         ,'output','none','error',err,'Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);

%% even

   te  = linspace(tstart,tend,nt+(rem(nt,2)));
   te  = te + (rand(size(te))-.5).*rnd;
   dte = (te(2) - te(1)).*24;
   te(1)=tstart;te(te<tstart) = tstart; te(te>tend) = tend;te(end)=tend;te = sort(te);
   ze  = sin(2*pi.*te./T);
   disp('-----even---------t----------------------')
   [De1,Ee1] = t_tide(ze,'int',diff(te).*24,'output','none','error',err,'Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   disp('-----even---------dt---------------------')
   [De2,Ee2] = t_tide(ze,'int',dte         ,'output','none','error',err,'Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   
%% test forced crash when dt variation exceeds tolerance 1e-3
%  except when error measure is 'linear'

   disp('-----odd----------dt rand----------------')
   %diff(to)
   to3 = to + (rand(size(to))-.5).*dto./dtfac;
   to3(1)=tstart;to3(to3<tstart) = tstart; to3(to3>tend) = tend;to3(end)=tend;to3 = sort(to3);
   zo3 = sin(2*pi.*to3./T);
   %diff(to3)
   try
   [Do3,Eo3] = t_tide(zo3,'int',diff(to3).*24,'output','none','error',err  ,'Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   OK(end+1) = 0;
   catch
   [Do3,Eo3] = t_tide(zo3,'int',diff(to3).*24,'output','none','error','lin','Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   OK(end+1) = 1;
   end
   
   try
   [Do4,Eo4] = t_tide(zo3,'int',diff(to3).*24,'output','none','error','lin','Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   OK(end+1) = 1;
   end

   disp('-----even---------dt rand----------------')
   %diff(te)
   te3 = te + (rand(size(te))-.5).*dte./dtfac;
   
   te3(1)=tstart;te3(te3<tstart) = tstart; te3(te3>tend) = tend;te3(end)=tend;te3 = sort(te3);
   ze3 = sin(2*pi.*te3./T);
   %diff(te3)
   try
   [De3,Ee3] = t_tide(ze3,'int',diff(te3).*24,'output','none','error',err  ,'Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   OK(end+1) = 0;
   catch
   [De3,Ee3] = t_tide(ze3,'int',diff(te3).*24,'output','none','error','lin','Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   OK(end+1) = 1;
   end

   try
   [De4,Ee4] = t_tide(ze3,'int',diff(te3).*24,'output','none','error','lin','Rayleight','S2','lsq',lsq,'start',start_time,'syn',synthesis);
   OK(end+1) = 1;
   end

%% do numeric machine tests

   % fits or synthetic reconstructions
   eps_series = eps('single');
   OK(end+1) = max(abs(Eo1-zo )) < eps_series;
   OK(end+1) = max(abs(Eo2-zo )) < eps_series;
   OK(end+1) = max(abs(Ee1-ze )) < eps_series;
   OK(end+1) = max(abs(Ee2-ze )) < eps_series;

   OK(end+1) = max(abs(Eo3-zo3)) < eps_series;
   OK(end+1) = max(abs(Eo4-zo3)) < eps_series;
   OK(end+1) = max(abs(Ee3-ze3)) < eps_series;
   OK(end+1) = max(abs(Ee4-ze3)) < eps_series;
   
   % amplitudes
   eps_amp = eps('single');
   OK(end+1) = max(abs(Do1.tidecon(1) - 1)) < eps_amp;
   OK(end+1) = max(abs(Do2.tidecon(1) - 1)) < eps_amp;
   OK(end+1) = max(abs(De1.tidecon(1) - 1)) < eps_amp;
   OK(end+1) = max(abs(De2.tidecon(1) - 1)) < eps_amp;

   OK(end+1) = max(abs(Do3.tidecon(1) - 1)) < eps_amp;
   OK(end+1) = max(abs(Do4.tidecon(1) - 1)) < eps_amp;
   OK(end+1) = max(abs(De3.tidecon(1) - 1)) < eps_amp;
   OK(end+1) = max(abs(De4.tidecon(1) - 1)) < eps_amp;

   % amplitudes
   % all phases are only 90 when start time is supplied,
   % otherwise start times are ragged
   eps_phase = 1e2.*eps('single');
   if isempty(start_time)
       OK = [OK 2 2 2 2  2 2 2 2]; % undefined, depends on internal centraltime
   else
   OK(end+1) = max(abs(Do1.tidecon(3) - 90)) < eps_phase;
   OK(end+1) = max(abs(Do2.tidecon(3) - 90)) < eps_phase;
   OK(end+1) = max(abs(De1.tidecon(3) - 90)) < eps_phase;
   OK(end+1) = max(abs(De2.tidecon(3) - 90)) < eps_phase;

   OK(end+1) = max(abs(Do3.tidecon(3) - 90)) < eps_phase;
   OK(end+1) = max(abs(Do4.tidecon(3) - 90)) < eps_phase;
   OK(end+1) = max(abs(De3.tidecon(3) - 90)) < eps_phase;
   OK(end+1) = max(abs(De4.tidecon(3) - 90)) < eps_phase;
   end
   
   disp('Amplitude   phase')
   disp(num2str(Do1.tidecon([1 3])))
   disp(num2str(Do2.tidecon([1 3])))
   disp(num2str(De1.tidecon([1 3])))
   disp(num2str(De2.tidecon([1 3])))
   disp(num2str(Do3.tidecon([1 3])))
   disp(num2str(Do4.tidecon([1 3])))
   disp(num2str(De3.tidecon([1 3])))
   disp(num2str(De4.tidecon([1 3])))
   
   
%% test plot

   if doplot
   TMP = figure('name',mfilename);maximize(TMP);
   plot(to ,zo ,'ro--','DisplayName','nobs=odd,  raw','markersize',4,'markerfacecolor','r')
   hold on
   plot(to ,Eo1,'ro-','DisplayName','nobs=odd,  t','markersize',8)
   plot(to ,Eo2,'ro-','DisplayName','nobs=odd,  dt~constant','markersize',14)
   plot(to3,zo3,'rv:','DisplayName','nobs=odd,  dt variable raw','markersize',7,'markerfacecolor','r')
   plot(to3,Eo4,'rv:','DisplayName','nobs=odd,  dt variable','markersize',14)
   
   plot(te ,ze ,'bo--','DisplayName','nobs=even, raw','markersize',4,'markerfacecolor','b')
   plot(te ,Ee1,'bo-','DisplayName','nobs=even, t','markersize',8)
   plot(te ,Ee2,'bo-','DisplayName','nobs=even, dt~constant','markersize',14)
   plot(te3,ze3,'bv:','DisplayName','nobs=even, dt variable raw','markersize',7,'markerfacecolor','b')
   plot(te3,Ee4,'bv:','DisplayName','nobs=even, dt variable','markersize',14)
   
   title(OKlegend{i,j,k})
   legend show
   grid on
   pausedisp
   try,close(TMP),end
   end

   OKall{i,j,k} = OK;OK=[];
   
end % lsq
end % synthesis
end % start_time

%num2str((to ([1 end])))
%num2str((te ([1 end])))
%num2str((to3([1 end])))
%num2str((te3([1 end])))

disp('    timesrs          phases ')
disp('try         amplitde        ')
disp('----########--------########')
disp('                            ')
disp('    - means     # means that')
disp('    nt=odd      nt=even     ')
disp('--##----####----####----####')
disp('                            ')
disp('    - means     # means that')
disp('    dt small    DT LARGE    ')
disp('    ~ constant  O(dt)       ')
disp('####--##--##--##--##--##--##')
for i=1:size(OKall,1)
for j=1:size(OKall,2)
for k=1:size(OKall,3)
    disp([num2str(OKall{i,j,k},'%1d'),'   ',OKlegend{i,j,k}])
end
end
end

% 1111111111111111111122222222   lsq=normal    synthesis=-1    start\_time=
% 1111111111111111111111111111   lsq=normal    synthesis=-1    start\_time=0
% 1111111111111111111111111111   lsq=normal    synthesis=-1    start\_time=718807
% 1111111111111111111122222222   lsq=normal    synthesis=0    start\_time=
% 1111111111111111111111111111   lsq=normal    synthesis=0    start\_time=0
% 1111111111111111111111111111   lsq=normal    synthesis=0    start\_time=718807
% 1111111111111111111122222222   lsq=normal    synthesis=1    start\_time=
% 1111111111111111111111111111   lsq=normal    synthesis=1    start\_time=0
% 1111111111111111111111111111   lsq=normal    synthesis=1    start\_time=718807
% 1111111111111111111122222222   lsq=direct    synthesis=-1    start\_time=
% 1111111111111111111111111111   lsq=direct    synthesis=-1    start\_time=0
% 1111111111111111111111111111   lsq=direct    synthesis=-1    start\_time=718807
% 1111111111111111111122222222   lsq=direct    synthesis=0    start\_time=
% 1111111111111111111111111111   lsq=direct    synthesis=0    start\_time=0
% 1111111111111111111111111111   lsq=direct    synthesis=0    start\_time=718807
% 1111111111111111111122222222   lsq=direct    synthesis=1    start\_time=
% 1111111111111111111111111111   lsq=direct    synthesis=1    start\_time=0
% 1111111111111111111111111111   lsq=direct    synthesis=1    start\_time=718807