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');
