% This routine calculates velocities standart deviation based on the error
% velocity from ADCP data.
%
% First run the ProcED_II routine in order to obtain 
% tidal_cycle_matrix, dist_i, depth_i, files
%
% use the same ProcED_II input file
%
% The following routine files must be in the same directory as ProcED_II
% BlWhRe_colormap, rot.m, smartmean.m, isoline.m, filter_5.m
%
% for full information see: 
% Genz,F., Cirano, M., Lessa, G.C. 2010.ProcED: a MATLAB package for 
% Processing ADCP Estuarine Data. Revista Brasileira de Geofísica
% 28(2):183-192.
%
% for information on how to run, data format and initial steps, please
% refer to the manual: "ProcED_II_user_Manual.doc"
%
% - version II - upgraded by Ricardo M. Domingues , ricardormd@gmail.com
%
% looping indexes:
% i - depth 
% j - cross section
% k - time


clear all; close all;

dir_path = input('\nEnter the full path of data directory \n ==> ','s');
%Change the current directory to the input path
cd(dir_path);
%---ADCP_ProcED_processing_options file--------------
%           
%   - refer to the manual
naminfo=dir('ProcED_processing_options.txt');
if isempty(naminfo),
  fprintf('ERROR******* Can''t find file ProcED_processing_options');
  return;
end;

% plot type set
style=3;

%%%%%%%%%%%%%%%%%%%
% read input file %
%%%%%%%%%%%%%%%%%%%
fid_Op=fopen('ProcED_processing_options.txt');
fgetl(fid_Op);files_header=fgetl(fid_Op);
fgetl(fid_Op);label=fgetl(fid_Op);
fgetl(fid_Op);flag=fscanf(fid_Op,'%f');
fgetl(fid_Op);ini_margin=fscanf(fid_Op,'%f');
fgetl(fid_Op);dist_inicio=fscanf(fid_Op,'%f');
fgetl(fid_Op);ang_manual=fscanf(fid_Op,'%f');
fgetl(fid_Op);ctds_inf=fscanf(fid_Op,'%f');
fgetl(fid_Op);dist_ctd=fscanf(fid_Op,'%f');
fclose(fid_Op);

if(ctds_inf>0);n_ctd=ctds_inf;ctds_inf=1;end

clear fid_Op naminfo

files=dir([files_header,'*']);
n_perf=length(files);
if(ini_margin==1),inverte=zeros(1,n_perf);else, inverte=ones(1,n_perf);end
for fl=1:n_perf,file_names{fl}=files(fl).name;end
if(length(dist_inicio)==1),dist_inicio(1:n_perf)=dist_inicio;end

clear files fl files_header
   
% Defining the size of the matrix

load tidal_cycle_matrix dist_i depth_i

    n_vert=length(dist_i); % number of columns
    
    c_vert=depth_i(1)-depth_i(2); % vertical resolution

    z=0:0.1:1;                              % adimensional Z depths
    
    masknan=ones(length(depth_i),length(dist_i),n_perf);
    mask=ones(length(depth_i),length(dist_i),n_perf);
    speed_i=ones(length(depth_i),length(dist_i),n_perf);
    V_z=ones(length(z),length(dist_i),n_perf);

% read ADCP data for the time k

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1st K looping
for k=1:n_perf

    fprintf('\n ******************************************************************** \n\n');    
    file = [dir_path,file_names{k}];
    T = dlmread(file);T = T';%T = load(file);
    
    time=T(1:7,:);
    dx=T(9,:);
    dy=T(8,:);   
    m=(size(T,1)-10)/4;
    
    %Calculating DT for each section
    
    hour(k,1)= time(4,1) + time(5,1)/60;
    hour_i(k,1)=hour(k,1);
    dia = time(3,1);ano=time(1,1);mes=time(2,1);
        
    if(k==1)
        dia_ini = dia;
        data{1}= sprintf('%02g/%02g/%02g',mes,dia,ano);
    elseif(dia>dia_ini)
        data{2}= sprintf('%02g/%02g/%02g',mes,dia,ano);
    end             
    
    if(k>1 && hour(k,1)<hour(k-1,1))
        hour_i(k,1)=hour(k-1,1)+hour(k,1);hour_v = hour(k-1,1);
    elseif(k>1 && hour(k,1)<hour_i(k-1,1)),hour_i(k,1)=hour_i(k-1,1)+ ...
            (hour(k,1)-hour(k-1,1));
    end
    
    DT_sec(k)= (time(3,end)-time(3,1))*24 + [(time(4,end) + time(5,end)/60)...
        - (time(4,1) + time(5,1)/60)];
    hour_end(k)= hour_i(k)+ DT_sec(k);
    
    %Calculating the river angle based on the section angle
    a=mean(diff(dy));b=mean(diff(dx)); ang_sec = (atand(abs(b/a)));

    if(a > 0 && b > 0),ang_sec=ang_sec;end
    if(a < 0 && b > 0),ang_sec=ang_sec + 90;end
    if(a < 0 && b < 0),ang_sec=ang_sec + 180;end
    if(a > 0 && b < 0),ang_sec=ang_sec + 270;end
              
    if(ini_margin==1)
        ang_rio(k)=ang_sec - 90;
    else
        ang_rio(k)=ang_sec + 90;
        if(ang_rio(k)>360),ang_rio(k)=ang_rio(k)-360;end
    end
    %Defining which sections must be inverted
    
    if(k>1 && abs(ang_rio(k)-ang_rio(1))>140)
        if(ang_rio(k)<180)
            ang_rio(k)=ang_rio(k)+180;
        else
            ang_rio(k)=ang_rio(k)-180;
        end
        if(ini_margin==1),inverte(k)=1;else,inverte(k)=0;end
    end                     
    
    % sum distance to the bouy in the cross profile beginning plus ensemble
    % avg distance
    dist_tmp=sqrt(dx.^2+dy.^2) + dist_inicio(k);
    dist = dist_tmp + mean(diff(dist_tmp));
    u=T(11+2*m:11+3*m-1,:);
    v=T(11+m:11+2*m-1,:);
    [m,n]=size(v);Cont=find(v==flag);  
    
    clear u v
    
    depth=-1.*T(11:11+m-1,1);
    depth_int= mean(diff(depth));
    depth_max=-1*T(10,:)';
    ncol1=length(dist);
    nline1=length(depth);
    %%%%%% load error velocity %%%%%%%
    speed=T(11+3*m:11+4*m-1,:);
    speed(Cont)=nan;
    clear T m n
% check if there is ensemble without bottom track and use interpolation to
% fill the gap
    ind0=isnan(depth_max);
    ind1=find(ind0==1);
    ind2=find(ind0==0);
    if ind1 > 0
    depth_max(ind1)=interp1(dist(ind2),depth_max(ind2),dist(ind1),'linear','extrap');
    end
    clear file; clear adcp; clear dx; clear dy; clear u; clear v;
    clear v_cross; clear v_along; clear ind0; clear ind1; clear ind2;

% check if it is necessary to reverse the cross profile direction

    if inverte(k)==1
        speed_tmp(:,1:ncol1)=speed(:,ncol1:-1:1);
        speed=speed_tmp;
        depth_max_tmp(1:ncol1)=depth_max(ncol1:-1:1);
        depth_max= depth_max_tmp;
        for j=1:ncol1
         dist_tmp_inv(j)=mean(diff(dist_tmp))+dist_tmp(ncol1)...
         -dist_tmp(ncol1-j+1);
        end
        dist=dist_tmp_inv;
        clear dist_tmp_inv; clear speed_tmp; clear dist_tmp; clear depth_max_tmp;
    end

% matrix redimension
    for j=1:ncol1
        ind_tmp=isnan(speed(:,j));
        ind0=find(ind_tmp==0);
        if length(ind0)>=1 % check if column has data 
            ind1(j)=1;
        else
            ind1(j)=0;
        end
    end
    ind2= find(ind1==1); % find god data columns
    ncol2= length(ind2); % how many
    speed_z=zeros(length(depth_i),ncol2);
    clear ind2;
    
% find bad ensemble columns
    n_bad= find(ind1==0);
    if length(n_bad) < 1; n_bad=0; end
    bad_seq= diff(n_bad);
    c=0;
    for kk=1:length(bad_seq)
        if bad_seq(kk) == 1
        c=c+1;
        end
    end

fprintf('\nCross profile n. : %d   \n',k);
fprintf('Bad ensemble id :  %d   \n',n_bad);
fprintf('N. of sequencial bad ensembles : %d \n\n',c);
clear ind0; clear ind1; clear n_bad; clear bad_seq; clear kk; clear c;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set matrix speed to the god column data in speed_z   %
% vertical INTERPOLATION - top and bottom              %
% Eliminate cells in the Side Lobe Layer Pulse or      %
% one cell close to bottom                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   Dtop=-depth(1); %depth to the middle of the first cell
   Dsurface= Dtop - (-depth_int)/2; % depth of de Top Discharge Layer

   counter=0; 
   for j=1:ncol1
        ind_tmp=isnan(speed(:,j));
        ind0=find(ind_tmp==0);
        ind1=find(ind_tmp==1);
        speed(ind1,j)=0; % set NaN to ZERO
        if length(ind0)>=1 
           counter=counter+1;
           maxdepth_col(counter)=depth_max(j);
           dist_col(counter)=dist(j);
           Dlg_max=-depth(ind0(end)); % WinRiver II ASCII files delete cells after Dlg_max
                                       % jul/2010
          
           % set bottom layer cells
           ind_pos=find(-depth > Dlg_max & -depth <= -maxdepth_col(counter)); 
        
           if length(ind_pos) >= 1 && ind_pos(1) > 1
               n_pos=length(ind_pos);
               depth_alt(1:ind_pos(1)-1)=depth(1:ind_pos(1)-1); % set cells before Dlg
               speed_alt(1:ind_pos(1)-1,counter)=speed(1:ind_pos(1)-1,j);
               depth_alt(ind_pos(1))=maxdepth_col(counter); % cell after Dlg => bottom depth
               speed_alt(ind_pos(1),counter)=0;                   % speed = 0 at bottom;
               nn1=length(depth)-ind_pos(end);
                for i=1:nn1                              % cell after bottom, modified depth´s
                    depth_alt(ind_pos(1)+i)=depth(ind_pos(end)+i);
                    speed_alt(ind_pos(1)+i,counter)=0;         % speed = 0 after bottom;
                end
                
                for i=1:n_pos-1                    % set last cells depth and speed (null)
                     depth_alt(nline1-n_pos+1+i)=depth(end)+i*depth_int;
                     speed_alt(nline1-n_pos+1+i,counter)=0;
                end
               
           else
               depth_alt=depth;
               speed_alt(:,counter)=speed(:,j);
           end
           
           % new depth interpolation and to the top extrapolation
           speed_z(:,counter)=interp1(depth_alt,speed_alt(:,counter),depth_i','linear');
           speed_top(:,counter)=interp1(depth_alt,speed_alt(:,counter),depth_i','nearest','extrap');
           dif_depth_top=-depth(1)/c_vert;
           ncell_top=floor(dif_depth_top);
           for dd=1:ncell_top+1
               speed_z(dd,counter)=speed_top(dd,counter); % usa a extrap nearest no topo
           end
           
           clear ind_pos; clear depth_alt; clear speed_alt; 
           clear n_pos; clear dd;
        end
        clear ind_tmp; clear ind0; clear ind1; clear dif_depth_top; clear ncell_top;
        clear Dlg_max; clear Dsl; clear nn1;  clear speed_pos;
                
   end

   clear counter; clear maxdepth_tmp; clear depth; clear Dlg_max; clear speed_top;
   clear speed_bott; clear depth_max; clear dist;

%%%%%%%%%%%%%%%%%%%%%%%%
% interpolate the grid %
%%%%%%%%%%%%%%%%%%%%%%%%

speed_i(:,:,k) = interp2(dist_col,depth_i',speed_z,dist_i,depth_i','linear',0);
nline_i= length(depth_i);
ncol_i=length(speed_i(1,:,k));
clear speed_z; clear speed;

% max depth interpolation in the new matrix
maxdepth_i(k,:)=interp1(dist_col,maxdepth_col,dist_i,'linear','extrap');

izero=find(maxdepth_i(k,:)>0); % verifica se houve extrapolacao incorreta
if izero > 0; maxdepth_i(k,izero)=0; end

clear maxdepth_col; clear izero; clear dist_col;

% calculate velocity in Z adimensional depths

for j=1:ncol_i
    ind0=find(depth_i >= maxdepth_i(k,j));
    if sum(ind0 > 1)
    z_tmp=depth_i(ind0)/maxdepth_i(k,j);
    V_z(:,j,k)= interp1(z_tmp',speed_i(ind0,j,k),z','linear','extrap');
    clear ind0; clear z_tmp;
    else
    V_z(:,j,k)= NaN;
    end
end

% define the cross profile mask

    for j=1:ncol_i
    ind0=find(depth_i < maxdepth_i(k,j)); % all cells below max depth are set to NaN
    masknan(ind0,j,k)=NaN;
    clear ind0;
    end

mask(:,:,k)=speed_i(:,:,k).*masknan(:,:,k);


% plot cross section velocity field

%     figure(k);
%     contourf(dist_i,depth_i,mask(:,:,k),'k');
%     caxis([-0.2 0.2]);
%     colorbar;
%     hold on;
%     plot(dist_i,maxdepth_i(k,:),'-k','linewidth',2)
%     hold on;
%    
%     title('ERROR VELOCITY field');
%     xlabel('Distance (m)');
%     ylabel('Depth (m)');

    
    %%#############################################%%
    % estimate ERROR - adcp measured data
    %%#############################################%%

    tmp_matrix=mask(:,:,k);

    ind0=isnan(tmp_matrix);
    ind1=find(ind0==0);
    speed_error(k)=std(tmp_matrix(ind1));
    
    clear ind0; clear ind1;  
    clear tmp_matrix;
       
end %%%%%%%%%%%%%%%%% 1st k looping %%%%%%%%%%%%%%%%%%%%%

min_speed_error=min(speed_error);
max_speed_error=max(speed_error);
   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET AVERAGED CROSS PROFILE VEL FOR THE SAME TIME, DISCHARGE AND VEL     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diff_hour=ones(n_perf-1,1);
for k=2:n_perf
    if((hour_i(k)-hour_end(k-1))<DT_sec)
        diff_hour(k-1)=0;
    end
end

hour=hour_i;
ind_hour=find(diff_hour ~= 0);
last_hour=max(ind_hour)+1;
if ind_hour >=1
    hour_i= hour(ind_hour);
    hour_i(length(ind_hour)+1)=hour(last_hour);
else
    hour_i=hour(1);
end

n_time=length(hour_i); 
n_z=length(z);

% Calculate averaged cross profile vel for the same time (if exists)
    
for k=1:n_time %%%%%%%%%%%%%%%%%%%%% 2nd K Looping %%%%%%%%%%%%%%%%%%%
    
    ind1=find(hour==hour_i(k));

    if length(ind1)>1
    speed_profile(:,:,k)=smartmean(mask(:,:,ind1),3);
    V_z_hour(:,:,k)=smartmean(V_z(:,:,ind1),3);
    maxdepth_hour(k,:)=smartmean(maxdepth_i(ind1,:),1);
    else
    speed_profile(:,:,k)=mask(:,:,ind1);
    V_z_hour(:,:,k)=V_z(:,:,ind1);
    maxdepth_hour(k,:)=maxdepth_i(ind1,:);
    end
    
end %%%%%%%%%%%%%%%%%%%%% 2nd k looping %%%%%%%%%%%%%%%%%%%

  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE TIDAL CYCLE END AND SAME TIME STEP VEL FIELDS %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if n_time > 1
 
    %%%% load discharge informations %%%%%%%%
    load disch_res;

    % new time array

    new_n_time=length(new_hour); %%%%%%%
    
    % check time step
    chour=diff(new_hour);
    indh1=find(chour < 0);
    if length(indh1 > 0)
        new_hour(indh1+1:end)=new_hour(indh1+1:end)+24;   
    end
    
    % interpolate U in Zadm and maxdepth to the new time
    % calculate residual cross section field
     
     for j=1:ncol_i
         V_tmp(1:n_z,1:n_time)= reshape(V_z_hour(:,j,:),[length(V_z_hour(:,1,1)) length(V_z_hour(1,1,:))]);
         new_maxdepth(1:new_n_time,j) = interp1(hour_i,maxdepth_hour(:,j),new_hour','linear','extrap');
         ind1=isnan(V_tmp(1,:));
         ind2=find(ind1==0);
         for i=1:n_z
         new_V_z(i,j,1:new_n_time) = interp1(hour_i(ind2),V_tmp(i,ind2),new_hour','linear','extrap');
         net_speed(i,j)= (0.5*new_V_z(i,j,1)+ sum(new_V_z(i,j,2:end-1))+ 0.5*new_V_z(i,j,end))/24; % campo residual
         end
     end
 
  
    % set bottom cell vel to = 0 in Z adm profiles
 
     ind0=isnan(new_V_z(11,:,:));
     ind1=find(ind0==0);
     new_V_z(11,ind1)=0;
 
   
    %%#############################################%%
    % Estimate error - interpolated data
    %%#############################################%%
    
    for k2=1:new_n_time

    tmp_matrix=new_V_z(:,:,k2);
    ind0=isnan(tmp_matrix);
    ind1=find(ind0==0);
    new_speed_error(k2)=std(tmp_matrix(ind1));
    
    clear ind1; clear ind0; clear tmp_matrix;
        

    end % FIM CICLO K2

    clear speed_temp; clear speed_tmp2; clear m_tmp_1; clear m_tmp_2; clear diff_hour; clear ind_hour; clear last_hour;

end % 

min_new_speed_error=min(new_speed_error);
max_new_speed_error=max(new_speed_error);

new_error_perc=new_speed_error./abs(new_mean_speed);

new_discharge(2,:)=new_discharge(1,:).*(1+new_error_perc);
new_discharge(3,:)=new_discharge(1,:).*(1-new_error_perc);

% restore original time step
if length(indh1 > 0)
    new_hour(indh1+1:end)=new_hour(indh1+1:end)-24;   
end

% save informations

save error_res new_hour new_speed_error new_discharge


fout_name=[dir_path,'tidal_cycle_res_error.txt'];
ff=fopen(fout_name,'w');
var=[hour';speed_error];
fprintf(ff,'%s\n','Time   Error[m/s]');
fprintf(ff,'%3.1f %8.5f \n',var);
var1=[min_speed_error;max_speed_error];
fprintf(ff,'%s\n','Error min[m/s] Error max[m/s]');
fprintf(ff,'%8.3f %8.3f\n',var1);
fclose(ff);
clear fout_name; clear var;

fout_name=[dir_path,'disch_res_interp_error.txt'];
ff=fopen(fout_name,'w');
var=[new_hour;new_mean_speed;new_speed_error;new_discharge];
fprintf(ff,'%s\n','Time   U[m/s] Error[m/s] Disch[m3/s] Disch+Error Disch-Erro');
fprintf(ff,'%3.1f %6.4f %8.5f %8.3f %8.3f %8.3f\n',var);
var1=[min_new_speed_error;max_new_speed_error];
fprintf(ff,'%s\n','Error min[m/s] Error max[m/s]');
fprintf(ff,'%8.3f %8.3f\n',var1);
fclose(ff);
clear fout_name; clear var;clear var1;
%clear all
