%% User inputs
close all; clear all; clc;

Num2Average = 100; % Average these number of frames.
EndNumFrames = floor(1000/Num2Average); 
% The number of frames you hope to end up with.
folder_path = 'Cam_Date=211010_Time=142613_30percent_Fixed_blackpaint/VectorFields';

for i = 1:EndNumFrames
    beginPoint = (i-1)*Num2Average + 1;
    endPoint = (((i)*Num2Average)-1) + 1;
    img_mean = beginPoint:endPoint;
    disp(i)
    savefilename = fullfile(['142613_30percent_Fixed_blackpaint_' num2str(Num2Average) '_frames_set_' num2str(i, '%02d') '_of_' num2str(EndNumFrames)]);
    save_averaged_sections(savefilename, img_mean, folder_path)

end
% img_mean = 1:1000;
% savefilename = fullfile('142613_30percent_Fixed_blackpaint_1000frames');
% save_averaged_sections(savefilename, img_mean, folder_path)



load chirp.mat;
sound(y);
return;

function save_averaged_sections(savefilename, img_mean, folder_path)
nu_air = 1.544e-5; % m^2/s
%% Get list of all vc7 files inside folder
fileList = dir([folder_path '/*.vc7']);
vector_files = {fileList.name}';
nv = length(vector_files);
nt = length(img_mean);

% load first vector file to get size of the field of view
s_field = readimx([folder_path '/' vector_files{1}]);
vf = create2DVec(s_field.Frames{1});
x = vf.X;
z = vf.Y;
    
[nx,ny] = size(vf.Y);
u_all = zeros(nt,nx,ny);
w_all = u_all;
o_all = u_all;
g_all = u_all;
cc = 1; % counter for placing data into u_all, w_all

for i = img_mean
    %% Extract vector field data from vf structure to MATLAB variables
    s_field = readimx([folder_path '/' vector_files{i}]);
    vf = create2DVec(s_field.Frames{1});
    u = vf.U;
    w = vf.V;

    %% Remove outliers
    uol = isoutlier(u);
    u(uol) = 0;
    w(uol) = 0;

    wol = isoutlier(w);
    u(wol) = 0;
    w(wol) = 0;

    u(u==0) = nan;
    w(u==0) = nan;
    u(w==0) = nan;
    w(w==0) = nan;
    
    %% Calculate the curl of u-w
    % for some reason the meshgrid from DAVIS x, z doesn't work on the function
    % curl, so I had to make one using the MATLAB function meshgrid
    dx = x(2,1)-x(1,1);
    dz = z(1,1)-z(1,2);

    xc = [min(x(:)):dx:max(x(:))];
    zc = [min(z(:)):dz:max(z(:))];

    [xc,zc] = meshgrid(zc,xc);
    [omeg,gamma] = curl(xc,zc,u,w);
    
    % I think Martin just started trying this. It wasn't here last week. 
    
    % omeg = Vorticity_AnyStencil(dx,dz,u,w,1);
    %  Fernando Zigunov (2021). 2D Vorticity with any stencil
    %  (https://www.mathworks.com/matlabcentral/fileexchange/71612-2d-vorticity-with-any-stencil),
    %  MATLAB Central File Exchange. Retrieved September 14, 2021.   
    
    %% Place into u_all, w_all data files
    u_all(cc,:,:) = u;
    w_all(cc,:,:) = w;
    o_all(cc,:,:) = omeg;
    g_all(cc,:,:) = gamma;
    cc = cc+1;
end

%% Plot vector field
fs = 1; % Take every image. You can adjust this if there is something periodic in the data and you want to see it.
u_mean = squeeze(nanmean(u_all(1:fs:end,:,:),1));
w_mean = squeeze(nanmean(w_all(1:fs:end,:,:),1));
o_mean = squeeze(nanmean(o_all(1:fs:end,:,:),1));  % Why do you have two of these??  They are not the same as the ones calculated below.
g_mean = squeeze(nanmean(g_all(1:fs:end,:,:),1));

save(savefilename, 'x', 'z', 'u_mean', 'w_mean');  % You need to enclose the variable names with qoutations.
end