% function tha_test
%Testing subroutine for tha.m

%% Read measured data
% fid = fopen('Aug_9419750.txt','r');
% C = textscan(fid,'%4d/%2d/%2d %*s %d:%d %f','headerlines',15);
% fclose(fid);
fid = fopen('TimeSeries_9419750.txt','r');
C = textscan(fid,'%d %4d%2d%2d %d:%d %f %f %f','headerlines',8);
fclose(fid);

%% Simple Tidal Harmonic Analysis
% Prepare input
%for Aug_9419750.txt
% tday = datenum(double([C{1},C{2},C{3},C{4},C{5},C{1}*0])); %Time in Julien days
% in.x = C{6}; %Input water level time series
%for TimeSeries_9419750.txt
tday = datenum(double([C{2},C{3},C{4},C{5},C{6},C{1}*0])); %Time in Julien days
in.x = C{8}; %Input water level time series
in.t = (tday-tday(1))*24; %hrs with respect to zero
in.rayleigh = 0.5; 
% in.secular = 'none';
in.secular = 'linear';
in.tolerance = 0.01;
in.constituents = {'M2','S2','N2','K1','O1','P1'}; %use to select constituents
in.start_time = datevec(tday(1)); %Starting time vector [yyyy,mm,dd,HH,MM,SS]
%Run tha
out = tha(in);

%% Read NOAA constituents
fid = fopen('Constituents_9419750.txt','r');
C = textscan(fid,'%*s %s %f %f %f','headerlines',18);
fclose(fid);
cnamem = C{1}; ampm = C{2}; phasem = C{3}; speedm = C{4};
nfm = length(speedm);
mf = length(out.constituents);
ind = false(nfm,1);
for k=1:nfm
    for j=1:mf
        if strcmp(cnamem{k},out.constituents{j})
            ind(k) = true(1);
        end
    end
end
cnamem = cnamem(ind);
ampm = ampm(ind);
phasem = phasem(ind);

disp('Name     Measured       Calculated')
disp('        Amp   Phase    Amp   Phase')
for i=1:mf
    str = sprintf('%5s  %6.3f %6.2f  %6.3f %6.2f',cnamem{i},...
        ampm(i),phasem(i),out.amplitudes(i),out.phases(i));
    disp(str)
end

%% Predict tide using output constituents
xc = out.mean + in.x*0; %Initialize
for i=1:mf    
    xc = xc + out.nodal_factors(i)*out.amplitudes(i)*cos(out.speeds(i)*out.time ...
        + (out.equilibrium_arguments(i) - out.phases(i))*pi/180);
end

%% Plotting
clf
axes('position',[0.08, 0.61, 0.87, 0.335])
inc = 4;
plot(out.time(1:inc:end)/24,in.x(1:inc:end),'k.',...
    out.time/24,out.fit,'-b',...
    out.time(1:inc:end)/24,out.fit(1:inc:end)-in.x(1:inc:end),'r')
legend('Measured','Fitted','Residual',0)
% plot(out.time,in.x,'k.-',out.time,out.fit,'.-b',...
%     out.time,xc,'g.-',out.time,out.fit-in.x,'r')
% legend('Measured','Fitted','Predicted','Residual',0)
xlabel('Julien time [days]')
ylabel('Water Level [m]')
title('Crescent City water levels, August 2012')

axes('position',[0.08, 0.32, 0.87, 0.22])
b = [ampm,out.amplitudes];
bar(b,'grouped');
set(gca,'xticklabel',cnamem)
ylabel('Amplitude [m]')
legend('NOAA','Computed',0)

axes('position',[0.08, 0.055, 0.87, 0.22])
b = [phasem,out.phases];
bar(b,'grouped');
set(gca,'xticklabel',cnamem,'ylim',[0 360])
ylabel('Phases [deg]')
legend('NOAA','Computed',0)

return