%run_t_tide_loop.m run t_tide analysis & prediction 
%
%See also: tide, t_tide
clear all

station = {'Honolulu','HVH'};
lims    = {[-.5:.1:.6],[-2:.2:2]};
yr0     = [2007 2013];
package = 't_tide';
timer.analysis   = [];
timer.prediction = [];

for ist=1:2
        
%% load & analyse & store: <station>_noos_<yearfrom>_<package>_analysis.xml

    for iyr=1:2
        
    [D(iyr).time,D(iyr).data,D(iyr).meta] = noos_read(['..\measurements\',station{ist},'_noos_',num2str(yr0(iyr)),'_measurements.txt']);
     D(iyr).meta = matroos_noos_header(D(iyr).meta);
     tic;
    [A(iyr),~]=t_tide(D(iyr).data,...
        'interval',min(unique(diff(D(iyr).time)))*24,... % in hours
        'sort'    ,'-amp',...
        'start'   ,D(iyr).time(1),...
        'err'     ,'wboot',...
        'lsq'     ,'normal',...
        'output'  ,[station{ist},'_noos_',num2str(yr0(iyr)),'_',package,'_analysis.asc'],...
        'lat'     ,D(iyr).meta.lat);
     timer.analysis(ist) = toc;     
    end

    for iyr=1:2
       T(iyr) = tide_iho.from_t_tide_tidestruct(A(iyr));
       ihoname = [station{ist},'_noos_',num2str(yr0(iyr)),'_',package,'_analysis.xml'];
       T(iyr).to_iho_xml(ihoname);
    end

%% predict & store: <station>_noos_<yearfrom>_<yearto>_<package>_prediction.txt

    for iyr=1:2
    for i1=1:2
    P(i1,iyr).meta                 = D(i1).meta;
    P(i1,iyr).meta.tanalysis       = D(iyr).time(1);
    P(i1,iyr).meta.datenumanalysis = [];
    P(i1,iyr).time                 = D(i1).time(1):1./24./6:D(i1).time(end); % prediction always 10 minutes
    % NB lat required for nodal corrections to be taken into account
    tic;
    P(i1,iyr).data=t_predic(P(i1,iyr).time ,A(iyr).name,A(iyr).freq,A(iyr).tidecon,A(iyr).lat);
    if iyr==i1;timer.prediction(ist) = toc;end
    end
    end
 

    for iyr=1:2
    for i1=1:2
    noos_write(P(i1,iyr).time, P(i1,iyr).data, ...
         'headerlines',matroos_noos_header(P(i1,iyr).meta),...
         'filename'   ,[station{ist},'_noos_',num2str(yr0(iyr)),...
                                       '_',num2str(yr0(i1)),'_',package,'_prediction.txt']);
    end

%% plot    
    
    close
    lim = lims{ist}([1 end]);
%     subplot(2,1,1)
	plot(D(iyr).time,D(iyr).data,'Displayname','data')
    hold on
    plot(P(iyr,iyr).time, P(iyr,iyr).data,'r','Displayname','prediction')
    datetick('x')
    grid on
    title(station{ist})
    ylim(lim)
    
%     subplot(2,1,2)
%     scatterbin(D(iyr).data,P(iyr,iyr).data,lims{ist},lims{ist})
%     hold on
%     %plot(D(iyr).data,P(iyr,iyr).data,'k.') % takes too long to print
%     plot(lim,lim,'k--')
%     colorbar
%     axis equal
%     axis([lim lim])
%     grid on
    print2a4([station{ist},'_noos_',num2str(yr0(iyr)),'_',package,'_prediction.png'],'v','t',200,'o');
    
    end
end

%% write performance log

fid = fopen([package,'.log'],'w');
fprintf(fid,'station    %s %s\n',station{1},station{2});
fprintf(fid,'analysis   %f %f sec.\n',timer.analysis(1),timer.analysis(2));
fprintf(fid,'prediction %f %f sec.\n',timer.prediction(1),timer.prediction(2));
fclose(fid);