function [x,y,u_avg,v_avg,CHC_tot] = piv_averages(prefix,suffix,Nstart,Nfinish,interp)

% This program reads in a series of instantaneous PIV vector fields from

% Insight and averages them. The user has the option of excluding

% interpolated vectors, which have CHC > 1. (interp = 0 means do not

% interpolate, while interp = 1 means interpolate).

% Create file name for each image

c_exit=3.007; %speed of sound

x0=1024; %origin of the jet in pixels

y0=53.8; %origin of the jetin pixels

D=923.71; %diameter of jet exit in pixels

v_shift=0;

for i = Nstart:Nfinish

Nstring = int2str(i); % Convert iteration number to a character string

if i < 10

filename_inst = strcat(prefix,'0000',Nstring,suffix);

elseif i < 100

filename_inst = strcat(prefix,'000',Nstring,suffix);

elseif i < 1000

filename_inst = strcat(prefix,'00',Nstring,suffix);

elseif i < 10000

filename_inst = strcat(prefix,'0',Nstring,suffix);

else

filename_inst = strcat(prefix,Nstring,suffix);

end

% Read file name

A_inst = csvread(filename_inst,1,0);

x = A_inst(:,1); % x-position (mm)

y = A_inst(:,2); % y-position (mm)

u = A_inst(:,3); % x-velocity (m/s)

v = A_inst(:,4); % y-velocity (m/s)

chc = A_inst(:,5); % number of good vectors at this location

N = size(x,1); % Length of entire vector array

% Initialize output variables if this is the first file

if i == Nstart

u_tot = zeros(N,1);

v_tot = zeros(N,1);

CHC_tot = zeros(N,1);

end

for j = 1:N

if interp == 0

if chc(j,1) == 1

u_tot(j,1) = u_tot(j,1) + u(j,1);

v_tot(j,1) = v_tot(j,1) + v(j,1);

CHC_tot(j,1) = CHC_tot(j,1) + 1;

end

elseif interp == 1

if chc(j,1) > 0

u_tot(j,1) = u_tot(j,1) + u(j,1);

v_tot(j,1) = v_tot(j,1) + v(j,1);

CHC_tot(j,1) = CHC_tot(j,1) + 1;

end

end

end

end

for j = 1:N

u_avg(j,1) = u_tot(j,1)/CHC_tot(j,1);

v_avg(j,1) = v_tot(j,1)/CHC_tot(j,1);

end

% Set origin to jet exit centerline

x_c = x - x0;

y_c = y - y0;

% Shift by convective velocity

v = v - v_shift;

% Nondimensionalize variables

x_non = x_c/D; % Nondimensionalize using jet diameter

y_non = y_c/D; % Nondimensionalize using jet diameter

u_non = u_avg/c_exit; % Nondimensionalize using sonic speed

v_non = v_avg/c_exit; % Nondimensionalize using sonic speed

%%%%%%%FOR H/D=2%%%%%%%%%

% 1) Choose a threshold

yThreshold = 400; %accept all vectors that start at or above y = yThreshold;

% 2) identify all quiver arrows that meet the criteria

acceptIdx = y >= yThreshold;

% 3) plot the quiver arrows, but only the ones accepted

figure(4)

quiver(x(acceptIdx), y(acceptIdx), u_avg(acceptIdx), v_avg(acceptIdx), 5)

% add reference line at threshold if you'd like

line(x,yThreshold);

% Plot nondimensional vector field

figure(2)

quiver(x_non(acceptIdx),y_non(acceptIdx),u_non(acceptIdx),v_non(acceptIdx),4)

axis([-1 1 0 2])

xticks([-1 0 1])

yticks([0 1 2])

set(gca,'YtickLabel',2:-1:0)

set(gca,'XAxisLocation','top','YAxisLocation','left');

xlabel('x/D')

ylabel('y/d')

title('Nondimensional velocity field')

% Plot averaged vector field

figure(1)

quiver(x(acceptIdx),y(acceptIdx),u_avg(acceptIdx),v_avg(acceptIdx),4)

axis([0 2000 0 2000])

set(gca,'YtickLabel',2000:-200:0)

set(gca,'XAxisLocation','top','YAxisLocation','left');

xlabel('z (pixels)')

ylabel('x (pixels)')

title('Dimensional velocity field')

% Plot averaged vector field

figure(3)

quiver(x(acceptIdx)/48.11,(y(acceptIdx)/48.11),u_avg(acceptIdx),v_avg(acceptIdx),4)

axis([0 42 0 42])

yticks([0 5 10 15 20 25 30 35 40])

set(gca,'XAxisLocation','top','YAxisLocation','left');

% set(gca,'ydir','reverse')

xlabel('z (mm)')

ylabel('x (mm)')

title('Dimensional velocity field')

% Output averaged data

output_avg = [x y u_avg v_avg CHC_tot];

dlmwrite('average_vels.dat','xyUVC');

dlmwrite('average_vels.dat',output_avg,'-append');