Why won't my contour work?

6 ビュー (過去 30 日間)
Monique Embury
Monique Embury 2020 年 3 月 4 日
コメント済み: Monique Embury 2020 年 3 月 5 日
So I have been able to generatae a contour plot for a full data set. For one of my experiemnts there are bad vectors that should be removed from said data set. I have modified my code as follows but I cannot get the contour plot to work. It says "Z must be at least a 2x2 matrix."
Can someone help me please?
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=2.776; %speed of sound
x0=1040; %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)+0; % x-position (mm)
y = A_inst(:,2)-198.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
% define orifice opening in pixels
orificepl=linspace(0, x0-(D/2));
orificepr=linspace(x0+(D/2), 2016);
yorifice=2016*ones(size(orificepl));
% define orifice opening in mm
orificeml=linspace(0, (x0-(D/2))/48.11);
orificemr=linspace((x0+(D/2))/48.11, 2016/48.11);
% define orifice opening non-dimensional
orificendl=(orificepl-x0)/D;
orificendr=(orificepr-x0)/D;
% 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 = 300; %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(5)
% quiver(x(acceptIdx), y(acceptIdx), u_avg(acceptIdx), v_avg(acceptIdx), 5)
% % add reference line at threshold if you'd like
% line(x(acceptIdx),yThreshold,'linewidth', 4, 'color', 'k');
aa=(yThreshold/2-67.8)*ones(size(acceptIdx));
aa_non=((yThreshold/2)-y0-67.8)/D*ones(size(acceptIdx));
x_non_d=downsample(x_non(acceptIdx),15);
y_non_d=downsample(y_non(acceptIdx),15);
u_non_d=downsample(u_non(acceptIdx),15);
v_non_d=downsample(v_non(acceptIdx),15);
%for contour
N = size(x(acceptIdx),1); % Length of entire vector array
xsize = 1;
x1 = x(acceptIdx(1));
x2 = x(acceptIdx(2));
while (x2 > x1) % Calculate number of x-locations using while loop
xsize = xsize + 1;
x1 = x2;
x2 = x(acceptIdx(xsize+1));
end
ysize = N/xsize; % Calculate number of y-locations
% Parse x- and y- positions to avoid repetition
for i = 1:xsize
x_re(i,1) = x_non(acceptIdx(i,1));
end
for i = 1:ysize
y_re(i,1) = y_non(acceptIdx(1+(i-1)*xsize,1));
end
mag = sqrt((u_non(acceptIdx)).^2 + (v_non(acceptIdx)).^2);
% Reorganize data for contour plots
mag_re=[reshape(mag,xsize,ysize)]';
u_re = [reshape(u_non(acceptIdx),xsize,ysize)]';
v_re = [reshape(v_non(acceptIdx),xsize,ysize)]';
%Contour Quiver
figure(7)
contourf(x_re,y_re,mag_re);
cb = colorbar();
colormap(jet)
ylabel(cb,'|U_P_i_s_t_o_n/U|')
hold on
qp=quiver(x_non_d,y_non_d,u_non_d,v_non_d,3,'k')
line(orificendl,yorifice/1008, 'linewidth', 4, 'color', 'k')
line(orificendr,yorifice/1008, 'linewidth', 4, 'color', 'k')
line(x_non,aa_non, 'linewidth', 3, 'color', 'k') %only for H/d=2
hold off
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('z/D')
ylabel('x/D')
% title('Nondimensional velocity field')
set(findall(gcf,'-property','FontSize'),'FontSize',16)
end
A link to my data set can be found here since it is too large to upload.
Thank you in advance!
  4 件のコメント
Walter Roberson
Walter Roberson 2020 年 3 月 4 日
Well that's nice, but I cannot test it without knowing the input values and having the input files. You have provided one variable in the .mat file, and I, as an outside observer, cannot assume the value of any of the variables.
By the way, you can make your code shorter. Replace
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
with
Nstring = sprintf('%s%05d%s', prefix, i, suffix);
Monique Embury
Monique Embury 2020 年 3 月 5 日
Thank you for that tip! I have attached a link for a sample set of files I use. It will not let me upload them here due to their format.

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeVector Fields についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by