%This program is designed to process uc, vc, and wc into a presentable
%format:

factor = 0.877; %um per pixel
[yd,xd,zd] = size(uc); %dimensions of data

%replace data outliers with NaN's: assume average is within 2 of actual
%(low density of outliers)

%uc:
out = 0;
for h=1:yd
    average = nanmean(uc(:,h));
    outliers = find(uc(:,h) > average+2 | uc(:,h) < average-2);
    uc(outliers,h) = NaN;
    vc(outliers,h) = NaN;
    wc(outliers,h) = NaN;
    out = out+length(outliers)*3;
end

xn = x.*factor;
yn = y.*factor;
ucn = uc.*factor;
vcn = vc.*factor;
wcn = wc.*factor;

NNaN = 0;
for i =1:zd
%Loop through layers and count NaN values:
%Take rms average, std deviation, and plot:

%UC:
    B = isnan(ucn(:,:,i));
    n = find(B==1);
    U_RMS(i) = sqrt(nanmean(ucn(:).^2));
    U_STD(i) = nanstd(ucn(:));
    NNaN = NNaN + length(n);
   
    figure;
    pcolor(xn(:,:,i),yn(:,:,i),ucn(:,:,i));
    title(['\fontsize{20}X Displacement UC in Layer #',num2str(i),',um']),xlabel('\fontsize{20}X Position (um)'),ylabel('\fontsize{20}Y Position (um)');
    colorbar;
    %shading interp

%VC:
    B = isnan(vcn(:,:,i));
    n = find(B==1);
    V_RMS(i) = sqrt(nanmean(vcn(:).^2));
    V_STD(i) = nanstd(vcn(:));
    NNaN = NNaN + length(n);
    
    figure;
    pcolor(xn(:,:,i),yn(:,:,i),vcn(:,:,i));
    title(['\fontsize{20}Y Displacement VC in Layer #',num2str(i),',um']),xlabel('\fontsize{20}X Position (um)'),ylabel('\fontsize{20}Y Position (um)');
    colorbar;
    
%WC:
    B = isnan(wcn(:,:,i));
    n = find(B==1);
    W_RMS(i) = sqrt(nanmean(wcn(:).^2));
    W_STD(i) = nanstd(wcn(:));
    NNaN = NNaN + length(n);
    
    figure;
    pcolor(xn(:,:,i),yn(:,:,i),wcn(:,:,i));
    title(['\fontsize{20}Z Displacement WC in Layer #',num2str(i),',um']),xlabel('\fontsize{20}X Position (um)'),ylabel('\fontsize{20}Y Position (um)');
    colorbar;
end

%Calculate Local Strains (Lagrangian):
for i=2:xd-1
    
    for j=2:yd-1
        
        for k=1:zd
            
            %Partials of x
            if(xd>2)
                u_x = (uc(j,i+1,k)-uc(j,i-1,k))/(x(j,i+1,k)-x(j,i-1,k));
                v_x = (vc(j,i+1,k)-vc(j,i-1,k))/(x(j,i+1,k)-x(j,i-1,k));
                w_x = (wc(j,i+1,k)-wc(j,i-1,k))/(x(j,i+1,k)-x(j,i-1,k));
            end
            
            %Partials of y
             if(yd>2)

                u_y = (uc(j+1,i,k)-uc(j-1,i,k))/(y(j+1,i,k)-y(j-1,i,k));
                v_y = (vc(j+1,i,k)-vc(j-1,i,k))/(y(j+1,i,k)-y(j-1,i,k));
                w_y = (wc(j+1,i,k)-wc(j-1,i,k))/(y(j+1,i,k)-y(j-1,i,k));
             end
            
            %Partials of z (only w)
             if(zd>2 && k >1 && k<zd)

                w_z = (wc(j,i,k+1)-wc(j,i,k-1))/(z(j,i,k+1)-z(j,i,k-1));

             end
            e_xx(j,i,k) = u_x + 0.5*(u_x^2 + v_x^2 + w_x^2);
            e_yy(j,i,k) = v_y + 0.5*(u_y^2 + v_y^2 + w_y^2);
            e_xy(j,i,k) = 0.5*(u_y + v_x) + 0.5*(u_x*u_y + v_x*v_y + w_x*w_y);
            if((exist('w_z', 'var') ~= 0)) 
                e_zz(j,i,k) = w_z; 
            end
        end
    end
end

%trim zeros from starting at 2:
e_xx = e_xx(2:(yd-1),2:(xd-1),:);
e_yy = e_yy(2:(yd-1),2:(xd-1),:);
e_xy = e_xy(2:(yd-1),2:(xd-1),:);

%Export Results:
fprintf('\nNumber of Uncorrelated Subsets: %d\n',NNaN-out);
fprintf('\nNumber of Outliers/Badly Correlated Subsets: %d\n',out);
fprintf('\n(Note that these are replaced by NaN as well to keep from skewing averages\n');
fprintf('\nThey are defined as more than 2 pixels from the general trend of correlation)\n');

%total averages:

U_RMST = mean(U_RMS);
V_RMST = mean(V_RMS);
W_RMST = mean(W_RMS);
U_STDT = mean(U_STD);
V_STDT = mean(V_STD);
W_STDT = mean(W_STD);
B = isnan(e_xx);
e_xxRMS = mean2(e_xx(B==0));
e_xxSTD = std2(e_xx(B==0));
B = isnan(e_yy);
e_yyRMS = mean2(e_yy(B==0));
e_yySTD = std2(e_yy(B==0));
B = isnan(e_xy);
e_xyRMS = mean2(e_xy(B==0));
e_xySTD = std2(e_xy(B==0));
%B = isnan(e_zz);
%e_zzRMS = sqrt(mean2(e_zz(B==0).^2));
%e_zzSTD = std2(e_zz(B==0));

M = [U_RMS;V_RMS;W_RMS;U_STD;V_STD;W_STD];
N = [U_RMST,V_RMST,W_RMST;U_STDT,V_STDT,W_STDT;e_xxRMS,e_yyRMS,e_xyRMS;e_xxSTD,e_yySTD,e_xySTD];
E = [e_xx;e_yy;e_xy];

csvwrite('Averages.csv',M)
csvwrite('Totals.csv',N)
csvwrite('LocalStrains.csv',E)

% figure;
% pcolor(x(2:16,2:16,1),y(2:16,2:16,1),e_xx(:,:,1));
% title(['\fontsize{20}Strain in exx, Layer #',num2str(1)]),xlabel('\fontsize{20}X Position (um)'),ylabel('\fontsize{20}Y Position (um)');
% colorbar;
% 
% figure;
% pcolor(x(2:4,2:4,1),y(2:4,2:4,1),e_yy(:,:,1));
% title(['\fontsize{20}Strain in eyy, Layer #',num2str(1)]),xlabel('\fontsize{20}X Position (um)'),ylabel('\fontsize{20}Y Position (um)');
% colorbar;