% This routine calculates tidal characteristics (Genz et al. 2010)
% (U, Umax, ebb and flood duration, Assimetry index, tidal prism, etc.)
% It needs output files from ADCP_ProcED (disch_res.mat and disch_res2adcp.txt)
% and, if you have an error matrix, output files from ADCP_ProcED_error (error_res.mat).
%
% U = averaged cross-section velocity [m/s]
% u_bar = averaged vertical velocity [m/s] - barotropic
%
% day time

clear all;

% % open file
% fid= fopen('disch_res1.txt');
% R = textscan(fid,'%s%*[^\n]');
% fclose(fid);
% 
% % set the last survey time
% hora=str2double(R{1});
% hora_fim=hora(end);
% clear R; clear hora; clear fid;

% open file
fid= fopen('dacv_res.txt');
R = textscan(fid,'%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s');
fclose(fid);

% calculates Umax for ebb and flood tide
for i=1:25
u_bar(:,i)=str2double(R{i});
end
u_bar_max=max(nanmax(u_bar(4:end,:)));
u_bar_min=min(nanmin(u_bar(4:end,:)));
clear R; clear fid;

load disch_res

op1=input('\n ##### Do you have an error matrix? yes => 1 :   ');
if op1 == 1
	load error_res
end

%%%%%%%%% Calculates averaged U, Umax, tidal prism, IAdv, etc. %%%%%%%%%%

n=length(new_hour);
day=datenum(data);
if length(day) > 1
    chour=diff(new_hour);
    indh1=find(chour < 0);
    new_hour(indh1+1:end)=new_hour(indh1+1:end)+24;   
    dayhour=new_hour/24;
    daytime=day(1) + dayhour;
else
    dayhour=new_hour/24;
    daytime=day + dayhour;
end

% identify the moment of velocity signal change (ebb to flood or
% vice-versa) 
for i=1:n-1
mod_signal(i)=sign(new_mean_speed(i))+sign(new_mean_speed(i+1));
end

% find time where U=0
ind= find(mod_signal==0);
for i=1:length(ind)
    t_vel0(i)=interp1(new_mean_speed(ind(i):ind(i)+1),daytime(ind(i):ind(i)+1),0);
end

% calculats tidal cycle durations (total, ebb and flood)
t_total=abs(daytime(end)-daytime(1))*24;
t1=(t_vel0(2)-t_vel0(1))*24;
t2=t_total-t1;

% define dt
d=diff(daytime);
dt=d(1);

% set ebb and flood
a=find(daytime >= t_vel0(1) & daytime <= t_vel0(2));% ligado a t1
b=find(daytime < t_vel0(1) | daytime > t_vel0(2));% ligado a t2

% Velocities - std e max
Ustda=std(new_mean_speed(a));
Umxa=max(new_mean_speed(a));
Umna=min(new_mean_speed(a));
Ustdb=std(new_mean_speed(b));
Umxb=max(new_mean_speed(b));
Umnb=min(new_mean_speed(b));

% Discharge -  max
Qmxa=max(new_discharge(1,a));
Qmna=min(new_discharge(1,a));
Qmxb=max(new_discharge(1,b));
Qmnb=min(new_discharge(1,b));

nm1=length(new_discharge(:,1));

% calculates tidal prism and river discharge
    for i=1:nm1
    Ia(i,:)=trapz(new_discharge(i,a)); 
    end
    % reorganize data to be continue (if necessary)
    nn=length(b);
    Qmare_des=zeros(nm1,nn-1);
    Hmare_des=zeros(nm1,nn-1);
    ind1=find(daytime > t_vel0(2));
    ind2=find(daytime < t_vel0(1));
    n1=length(ind1);
    n2=length(ind2);
    Qmare_des(:,1:n1)=new_discharge(:,ind1);
    Qmare_des(:,n1:end)=new_discharge(:,ind2); % preserva o valor de T=1
    for i=1:nm1
    Ib(i,:)=trapz(Qmare_des(i,:));
    end
    % calculates initial/end triangle area - ebb ou flood tide
Ia1(:,:)=0.5*abs(t_vel0(1)-daytime(a(1)))*new_discharge(:,a(1))+0.5*abs(t_vel0(2)-daytime(a(end)))*new_discharge(:,a(end));
Ib1(:,:)=0.5*abs(t_vel0(1)-daytime(ind(1)))*new_discharge(:,ind(1))+0.5*abs(t_vel0(2)-daytime(ind(2)+1))*new_discharge(:,ind(2)+1);
Pa(:,:)=(Ia*dt+Ia1)*3600*24; %% prism
Pb(:,:)=(Ib*dt+Ib1)*3600*24; %% prism

%

% if the ebb tide is entire 
if sign(new_discharge(ind(1)))<0 
    fprintf('\nTidal cycle duration [hours] =  %3.2f ',t_total);
    fprintf('\nFlooding time =  %3.2f ',t2);
    fprintf('\nEbbing time =  %3.2f \n',t1);
    % IAdv
    Ua=mean(new_mean_speed(a));%ebb
    Ub=mean(new_mean_speed(b));
    Ad=(t1-t2)/(t1+t2);
    Av=(Ua-abs(Ub))/(abs(Ua)+abs(Ub));
    AIdv=Ad+Av
    prism=Pb
    res_disch=(Pa-abs(Pb))/3600/t_total 
    ind_mx_vz=find(new_discharge(1,:) == Qmxa);  
    ind_mx_en=find(new_discharge(1,:) == Qmnb);
    if nm1 > 1
		error_prism=(Pb(2)/Pb(1)-1)*100;
		error_res_disch=(res_disch(3)/res_disch(1)-1)*100;
    end
    % save the output file
    fout_name='tidal_synthesis.txt';
    ff=fopen(fout_name,'w');
    fprintf(ff,'Begin time = %s \n',datestr(daytime(1),'dd/mm/yy HH:MM'));
    fprintf(ff,'End time = % s \n',datestr(daytime(end),'dd/mm/yy HH:MM'));
    fprintf(ff,'Tidal cycle duration [hours] = %3.2f \n',t_total);
    fprintf(ff,'Low water time = %s \n',datestr(t_vel0(2),'dd/mm/yy HH:MM'));
    fprintf(ff,'High water time = %s \n',datestr(t_vel0(1),'dd/mm/yy HH:MM'));
    fprintf(ff,'Ebbing time (te) [hours] =  %3.2f \n',t1);
    fprintf(ff,'Ebb average U [m/s] =  %3.2f \n',Ua);
    fprintf(ff,'Ebb max U [m/s]=  %3.2f \n',Umxa);
    fprintf(ff,'Ebb max u_bar [m/s] =  %3.2f \n',u_bar_max);
    fprintf(ff,'Ebb Qmax [m3/s] =  %10.2f \n',Qmxa);
    fprintf(ff,'Ebb Qmax time = %s \n',datestr(daytime(ind_mx_vz),'dd/mm/yy HH:MM'));
    fprintf(ff,'Flooding time (tf) [hours] =  %3.2f \n',t2);
    fprintf(ff,'Flood average U [m/s] =  %3.2f \n',Ub);
    fprintf(ff,'Flood max U [m/s] =  %3.2f \n',Umnb);
    fprintf(ff,'Flood max u_bar [m/s] =  %3.2f \n',u_bar_min);
    fprintf(ff,'Flood Qmax [m3/s] =  %10.2f \n',Qmnb);
    fprintf(ff,'Flood Qmax time = %s \n',datestr(daytime(ind_mx_en),'dd/mm/yy HH:MM'));
    fprintf(ff,'tf/te =  %5.3f \n',t2/t1);
    fprintf(ff,'AIdv =  %5.3f \n',AIdv);
    fprintf(ff,'Ad =  %5.3f \n',Ad);
    fprintf(ff,'Av =  %5.3f \n',Av);    
    fprintf(ff,'Residual discharge =  %10.2f %10.2f %10.2f \n',res_disch);
    fprintf(ff,'Prism [m3] =  %12.1f %12.1f %12.1f \n',abs(Pb));
	if nm1 > 1
	fprintf(ff,'Prism Error (percent) =  %4.2f \n',error_prism);
	fprintf(ff,'Residual discharge Error (percent) =  %4.2f \n',abs(error_res_disch));
	end
    fclose(ff);

    % if the flood tide is entire
else 
    fprintf('\nTidal cycle duration [hours] =  %3.2f ',t_total);
    fprintf('\nFlooding time =  %3.2f ',t1);
    fprintf('\nEbbing time =  %3.2f \n',t2);
    % IAdv
    Ua=mean(new_mean_speed(a));%flood
    Ub=mean(new_mean_speed(b));
    Ad=(t2-t1)/(t1+t2);
    Av=(Ub-abs(Ua))/(abs(Ua)+abs(Ub));
    AIdv=Ad+Av
    prism=Pa
    res_disch=(Pb-abs(Pa))/3600/t_total %
    ind_mx_vz=find(new_discharge(1,:) == Qmxb); 
    ind_mx_en=find(new_discharge(1,:) == Qmna);
	if nm1 > 1
    error_prism=(Pa(2)/Pa(1)-1)*100;
    error_res_disch=(res_disch(3)/res_disch(1)-1)*100;
    end
	
    % save the output file
    fout_name='tidal_synthesis.txt';
    ff=fopen(fout_name,'w');
    fprintf(ff,'Begin time = %s \n',datestr(daytime(1),'dd/mm/yy HH:MM'));
    fprintf(ff,'End time = %s \n',datestr(daytime(end),'dd/mm/yy HH:MM'));
    fprintf(ff,'Tidal cycle duration [hours] = %3.2f \n',t_total);
    fprintf(ff,'Low water time = %s \n',datestr(t_vel0(1),'dd/mm/yy HH:MM'));
    fprintf(ff,'High water time = %s \n',datestr(t_vel0(2),'dd/mm/yy HH:MM'));
    fprintf(ff,'Ebbing time (te) [hours] =  %3.2f \n',t2);
    fprintf(ff,'Ebb average U [m/s] =  %3.2f \n',Ub);
    fprintf(ff,'Ebb max U [m/s]=  %3.2f \n',Umxb);
    fprintf(ff,'Ebb max u_bar [m/s] =  %3.2f \n',u_bar_max);
    fprintf(ff,'Ebb Qmax [m3/s] =  %10.2f \n',Qmxb);
    fprintf(ff,'Ebb Qmax time = %s \n',datestr(daytime(ind_mx_vz),'dd/mm/yy HH:MM'));
    fprintf(ff,'Flooding time (tf) [hours] =  %3.2f \n',t1);
    fprintf(ff,'Flood average U [m/s] =  %3.2f \n',Ua);
    fprintf(ff,'Flood max U [m/s] =  %3.2f \n',Umna);
    fprintf(ff,'Flood max u_bar [m/s] =  %3.2f \n',u_bar_min);
    fprintf(ff,'Flood Qmax [m3/s] =  %10.2f \n',Qmna);
    fprintf(ff,'Flood Qmax time = %s \n',datestr(daytime(ind_mx_en),'dd/mm/yy HH:MM'));
    fprintf(ff,'tf/te =  %5.3f \n',t1/t2);
    fprintf(ff,'AIdv =  %5.3f \n',AIdv);
    fprintf(ff,'Ad =  %5.3f \n',Ad);
    fprintf(ff,'Av =  %5.3f \n',Av);    
	fprintf(ff,'Residual discharge =  %10.2f %10.2f %10.2f \n',res_disch);
    fprintf(ff,'Prism [m3] =  %12.1f %12.1f %12.1f \n',abs(Pa));
    if nm1 > 1
	fprintf(ff,'Prism Error (percent) =  %4.2f \n',error_prism);
	fprintf(ff,'Residual discharge Error (percent) =  %4.2f \n',abs(error_res_disch));
	end
    fclose(ff);

    
end

figure
plot(daytime,new_discharge)
xlabel('Time');
ylabel('Discharge [m3/s]');


