function out = tha(in)
% tha Tidal Harmonic Analysis
%
% Usage:
%   out = tha(in);
%
% Description:
%   Performs a simple tidal harmonic analysis by fitting an equation
%   of the form:
%       xfit = xmean + f*a*cos(speed*t + v+u - phase)         (1)
%    where 
%       a - constituent amplitude
%       f - nodal factor
%       v+u - Equilibrium argument
%       phase - constituent phase
%
%   Complex time series are simplified into two separate linear regressions
%   similar to equation (1). 
%
% Input:
%   in - Structure array with the following fields
%     x - real or complex vector [arbitrary units]
%     t - time in hrs from zero [hours startin from zero]
%     start_time - starting time as [yyyy,mm,dd,HH,MM,SS] 
%     constituents - optional input constituent names
%     rayleigh - optional rayleigh criteria (<=1)
%     secular - optional secular correction. Either 'linear' or 'none'.
%                If set to 'linear', the slope is included as a correction
%                to the time series.
%
% Output:
%  out - Structure array with different fields depending on whether the 
%      input time series in real or complex.
%    amplitudes - Tidal constituent amplitudes (real time series)
%    constituents - Tidal constituent names
%    equilibrium_arguments - Equilibrium arguments (v+u) for each
%                            constituent and year
%    frequencies - Tidal constituent frequencies
%    frequency_units - 'cycles/hour'
%    nodal_factors - Nodal factors for each constituent and year.
%    phases - Tidal constituent phases
%    phase_units = 'degrees'    
%    percent_variance_explained - Percent variance explained of the
%                          measured time series by the fitted time series.
%    residuals - difference between measured and fit time series
%    slope - fitted slope(s)
%    speeds - Tidal constituent speeds in degrees/hour
%    speed_units = 'degrees/hour'
%    time - Julien time in hours from start of year
%    time_units = 'hours'
%    (variables bewlow are only for complex time series only) 
%    major_amplitude - Semi-major axis amplitude
%    minor_amplitude - Semi-minor axis amplitude
%    positive_amplitude - Positive axis amplitude (a+)
%    negative_amplitude - Negative axis amplitude (a-)
%    positive_phase - Positive axis phase (phase+) [deg]
%    negitive_phase - Negative axis phase (phase-) [deg]
%    ellipse_inclination - Ellipse orientation or inclination [deg]
%    eccentricity
%
% Disclaimer:
%  This program is provided as is without any guaranty or warranty.
%  The user assumes all risk involved in the application of this software.
%
% Version: 3.0
% Date: November 5, 2012
% Author: Alex Sanchez, USACE-CHL

if nargin==0
   tha_test
   return
end

%---- Set defaults ----
ray = 1.0;
tol = 0.0001;   
cnamein = []; 
secular = 'none';
tide_data_file = 'tide_data_v1.mat';

%--- Collect input options ----
if isfield(in,'rayleigh')
    ray = in.rayleigh;
end
if isfield(in,'tolerance')
    tol = in.tolerance;    
end
if isfield(in,'constituents')
    if iscell(in.constituents)
        cnamein = char(in.constituents);
    else
        cnamein = in.constituents;
    end 
end
if isfield(in,'secular')
    secular = in.secular;    
end

%--- Load Tidal Data -------------
load(tide_data_file)  %loads in tide.()
speedall = tide.speeds;
freqall = speedall/360; %cycles per hour
cnameall = tide.constituents;
nf = length(tide.speeds);

%--- Internal variables ---------------------
i = sqrt(-1);
stime = in.start_time;
%Calculate hrs with respect to start time
thr = (datenum(stime)-datenum(stime(1),1,1))*24 + in.t(:); 

%--- Remove NaN's -------
isf = isfinite(in.x);
xx = real(in.x(isf));
yy = imag(in.x(isf));
tt = thr(isf);
nt = length(tt);

if(~isempty(cnamein))
    %--- Extract specified names ------------------------------    
    ind = false(nf,1);
    for k=1:nf
        for j=1:length(cnamein)
            if strcmp(cnameall{k},cnamein(j,:))
                ind(k) = true(1);
            end
        end
    end
    ju = find(ind);    
else
    %--- Check Rayleigh Criteria -------------------------
    timespan = thr(end)-thr(1); %in hours 
    keep = ones(nf,1);
    for k=1:nf
        if(keep(k)==1)
            dfall = abs(freqall-freqall(k));
            ind = (dfall*timespan >= ray | dfall==0);
            keep = keep.*ind;
        end
    end
    ju = find(keep==1);
    % ju = 1:6;
end
mf = length(ju);
cname = cnameall(ju,:);
freq = freqall(ju);
speed = speedall(ju)*pi/180; %convert to radians
ind = (tide.years==stime(1));  %yearly data
f = tide.nodal_factors(ind,ju)';
vu = tide.equilibrium_arguments(ind,ju)';
r = range(xx);
bad = 1;
iter = 0;
while (bad==1 && iter<3)
    %--- Construct Matrix -------------------------------------
    if strcmpi(secular,'linear')
        M = [ones(nt,1), tt, cos(tt*speed'), sin(tt*speed')];
        nc = 2;
    else %none
        M = [ones(nt,1), cos(tt*speed'), sin(tt*speed')];
        nc = 1;
    end
    %--- Perform regressions ---------------------------
    [x0,dx,cx,sx,xfit,xres,pvx] = doreg(M,xx,mf,nc);
    a = sqrt(cx.^2 + sx.^2); 
    if ~isreal(in.x) %for complex time series only
        [y0,dy,cy,sy,yfit,yres,pvy] = doreg(M,yy,mu,nc);
        a = a + sqrt(cx.^2 + sx.^2); %used for while loop
    end
    %Check for small amplitudes        
    ind = (a<r*tol);
    if any(ind)
        mf = sum(~ind);
        cname = cname(~ind,:);
        freq = freq(~ind);
        speed = speed(~ind);
        f = f(~ind);
        vu = vu(~ind);
        bad = 1;
    else
        bad = 0;
    end
    iter = iter + 1;
end

%Process regresion results into output variables
if isreal(in.x)
    amp = a./f; %amplitude with nodal correction
    theta = atan2(-sx,cx);
    phase = mod(vu - theta*180/pi,360); %phase in degrees with correction
else %Ellipse parameters
    c1 = cx + sy;
    c2 = cy - sx;
    ap = 0.5*sqrt(c1.^2 + c2.^2); %a+ coefficients
    epsp = atan2(c2,c1); %positive phase angle, radians
    c1 = cx - sy;
    c2 = cy + sx;
    am = 0.5*sqrt(c1.^2 + c2.^2); %a- coefficients
    epsm = atan2(c2,c1); %negative phase angle, radians
    ap = ap./f; %positive (cyclonic) amplitude with nodal correction
    am = am./f; %negative (anti-cyclonic) amplitude with nodal correction
    fmaj = ap + am; %major axis amplitude
    fmin = ap - am; %minor axis amplitude
    gp = mod(vu - epsp*180/pi, 360); %positive (cyclonic) phase in degrees with nodal correction
    gm = mod(vu + epsm*180/pi, 360); %negative (anti-cyclonic) phase in degrees with nodal correction
    finc = mod((epsp + epsm)*180/pi, 360)/2; % ellipse inclination in degrees
    phase = mod(gp + finc, 360); %phase in degrees
    rot = (ap.^2 - am.^2)./(ap.^2 + am.^2); %Indicates rotation direction
    eccentricity = (ap-am)/(ap+am);
end

% %--- Predicted ------------------------------------------
% xc = x0 + xx*0; %Initialize
% for k=1:mf    
%     xc = xc + f(k)*amp(k)*cos(speed(k)*tt + (vu(k) - phase(k))*pi/180);
% end
% out.predicted = xc;

%--- Collect Output Variables -------------------------
out.time = tt; %hours
out.constituents = cname;
out.frequencies = freq;
out.speeds = speed;
out.nodal_factors = f;
out.equilibrium_arguments = vu;
out.time_units = 'hours';
out.frequency_units = 'cycles/hours';
out.speed_units = 'degrees/hour';
out.equilibrium_argument_units = 'degrees';
out.phase_units = 'degrees';
if isreal(in.x)    
    out.x  = xx;
    out.mean = x0;
    out.slope = dx;
    out.fit = xfit;
    out.residuals = xres;
    out.percent_variance_explained = pvx;
    out.amplitudes = amp;
    out.phases = phase;
else
    out.rotation = rot;
    out.fit = xfit +  i*yfit;
    out.mean = x0 + i*y0;
    out.slope = dx + i*dy;
    out.fit = xfit + i*fit;
    out.residuals = xres + i*yres;
    out.percent_variance_explained = pvx + i*pvy;
    out.major_amplitude = fmaj;
    out.minor_amplitude = fmin;
    out.positive_amplitude = ap;
    out.negative_amplitude = am;
    out.positive_phase = gp;
    out.negitive_phase = gm;
    out.ellipse_inclination = finc;
    out.phase = phase;   
    out.eccentricity = eccentricity;
end

return

%=========================================================
function [z0,dz,cz,sz,zfit,zres,pvz] = doreg(M,zz,mu,nc)
% Perform Regressions 
% written by Alex Sanchez, USACE-CHL
A = M\zz;  %Solve matrix
% A = robustfit(M(:,2:end),zz,'ols',4.685); %note that robustfit prepends M with a column
z0 = A(1); %Mean
if nc==2
    dz = A(2); %Trend
else
    dz = 0; %No trhend
end
n = 2*mu + nc; %Total number of coefficients
cz = A(nc+1:mu+nc); %cosine coefficients
sz = A(mu+nc+1:n);  %sine coefficients
zfit = M*A;  %Compute fit
zres = zz - zfit; %Residuals
pvz = 100*(1-cov(zres)/cov(zz)); %Percent variance explained
return

%================================================================
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);

%Simple Tidal Harmonic Analysis
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
in.t = (tday-tday(1))*24; %hrs with respect to zero
in.rayleigh = 0.5; 
in.secular = 'none';
in.tolerance = 0.01;
in.constituents = {'M2','S2','N2','K1','O1'}; %use to select constituents
in.start_time = datevec(tday(1)); %Starting time vector [yyyy,mm,dd,HH,MM,SS]
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
subplot(3,1,1)
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')

subplot(3,1,2)
b = [ampm,out.amplitudes];
bar(b,'grouped');
set(gca,'xticklabel',cnamem)
ylabel('Amplitude [m]')
legend('NOAA','Computed',0)

subplot(3,1,3)
b = [phasem,out.phases];
bar(b,'grouped');
set(gca,'xticklabel',cnamem,'ylim',[0 360])
ylabel('Phases [deg]')
legend('NOAA','Computed',0)

return




