function []=calcsw(frame1, frame2)
% USAGE:
%
% swirl(frame1)
% swirl(frame1,frame2)
%
% Input Arguments:
% frame1 first frame in the range of frames for which swirl needs to be computed
% frame2 first frame in the range of frames for which swirl needs to be computed
%
% Output: Computed swirl data
%
% REMARK: If only one argument is provided then swirl is computed for that frame only. In addition it also
% generate a contour plot of the 2D swirl.
% If a range is provided then swirl is computed for all the frames in that range.
% The argument must be such that frame1 < frame2 and
%
% EXAMPLE:
% calcsw(100,200)
%
% This function computes the 2D swirl for all the frames between 100 upto 200.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% Copyright (C) 2005-2007 Anurag Singh %%
%% %%
%% This program/code snippet/file (hence forth refered as "library") %%
%% is free software; you can redistribute it and/or %%
%% modify it under the terms of the GNU Lesser General Public %%
%% License as published by the Free Software Foundation; either %%
%% version 2.1 of the License, or (at your option) any later version. %%
%% %%
%% This library is distributed in the hope that it will be useful, %%
%% but WITHOUT ANY WARRANTY; without even the implied warranty of %%
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU %%
%% Lesser General Public License for more details. %%
%% %%
%% You should have received a copy of the GNU Lesser General Public %%
%% License along with this library; if not, write to the Free Software %%
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA %%
%% %%
%% FILENAME: %%
%% calcsw.m %%
%% %%
%% AUTHOR: %%
%% %%
%% Anurag Singh, %%
%% MS, 2007 %%
%% Aerospace Engineering & Mechanics %%
%% University of Minnesota - Twin Cities. %%
%% Minneapolis, MN 55455 (USA) %%
%% %%
%% (currently working at: Exa Corporation, Burlington, MA 01803) %%
%% %%
%% CONTACT/EMAIL: %%
%% %%
%% anurag@aem.umn.edu %%
%% anurag9@gmail.com %%
%% %%
%% SOURCE CONTROL INFORMATION: %%
%% None (since I was planning on putting it under source control since it has %%
%% reached the critical file system size. Would be a good thing to put it under %%
%% source control while making changes. %%
%% %%
%% DESCRIPTION: %%
%% %%
%% Need to add %%
%% %%
%% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global u v x y NX sw2D vecimfreq u_norm piv_ustd u_tau nu vort2D dudx dudy dvdx dvdy stdsw2d stdvort2d showplot swcores cores
disp ([' SHOWPLOT Status = ',num2str(showplot),' ( 1 - display swirl plot, 0 - no plot display)']);
switch nargin
case 0
disp (' ERROR: Insufficient arguments.');
help calcsw;
case nargin*(nargin==1 || nargin==2)
if(nargin==1)
frame2 = frame1;
end
for frame=frame1:frame2
dx = x(2) - x(1);
dy = dx;
dudx = zeros(NX);
dudy = zeros(NX);
dvdx = zeros(NX);
dvdy = zeros(NX);
U = squeeze(u(frame,:,:));
V = squeeze(v(frame,:,:));
% Compute derivatives
% Compute all the inside points
% Compute the derivatives inside the boundaries
dudx(2:NX-1,:) = (U(3:NX,:) - U(1:NX-2,:))/(2*dx);
dudy(:,2:NX-1) = (U(:,3:NX) - U(:,1:NX-2))/(2*dy);
dvdx(2:NX-1,:) = (V(3:NX,:) - V(1:NX-2,:))/(2*dx);
dvdy(:,2:NX-1) = (V(:,3:NX) - V(:,1:NX-2))/(2*dy);
% Compute the edges
% left & bottom edges of the image, use 3 point forward difference
dudx(1,:) = ( 4*U(2,:) - 3*U(1,:) - U(3,:) )/(2*dx);
dudy(:,1) = ( 4*U(:,2) - 3*U(:,1) - U(:,3) )/(2*dy);
dvdx(1,:) = ( 4*V(2,:) - 3*V(1,:) - V(3,:) )/(2*dx);
dvdy(:,1) = ( 4*V(:,2) - 3*V(:,1) - V(:,3) )/(2*dy);
% Top & right edges of the image, use 3 point backward difference
dudx(NX,:) = (-4*U(NX-1,:) + 3*U(NX,:) + U(NX-2,:) )/(2*dx);
dudy(:,NX) = (-4*U(:,NX-1) + 3*U(:,NX) + U(:,NX-2) )/(2*dy);
dvdx(NX,:) = (-4*V(NX-1,:) + 3*V(NX,:) + V(NX-2,:) )/(2*dx);
dvdy(:,NX) = (-4*V(:,NX-1) + 3*V(:,NX) + V(:,NX-2) )/(2*dy);
% Compute the descriminant matrix at all the points and hence the 2D swirl
for i=1:NX
for j=1:NX
D = [dudx(i,j) dudy(i,j); dvdx(i,j) dvdy(i,j)];
P = -trace(D);
Q = det(D);
descr = P^2 - 4*Q;
if( descr < 0)
sw2D(frame,i,j) = sqrt(abs(descr))/2;
else
sw2D(frame,i,j) = 0;
end
vort2D(frame,i,j) = dvdx(i,j) - dudy(i,j);
end
end
end
stdsw2d = mean(mean(std(sw2D(frame1:frame2,:,:))));
stdvort2d = mean(mean(std(vort2D(frame1:frame2,:,:))))
% trial = abs(squeeze(vort2D(frame,:,:)));
% trial = trial/(u_tau^2/nu); % non dimensionalize
% trial3 = trial.*(trial>0.04); % cut-off threshold
% caxlim = [-3*stdvort2d 3*stdvort2d];
trial = abs(squeeze(sw2D(frame,:,:)));
trial = trial .* cores; % use cores identified by findcores_new previously
% maxsw = max(max(trial));
% trial3 = trial.*(trial> 880.8177*0.15 ); % cut-off threshold
% trial3 = trial.*(trial> 580*0.15 ); % cut-off threshold
caxlim = [100 7*stdsw2d];
trial2 = squeeze(vort2D(frame,:,:));
trial2 = trial2/(u_tau^2/nu); % non dimensionalize
%trial3 = trial2.*(abs(trial2)>0.04); % cut-off threshold
trial3 = trial .* trial2 ./ abs(trial2) ;
t1 = (trial3 > 0);
t2 = -1*(trial3 < 0);
trial3 = t1 + t2;
trial = squeeze(sw2D(frame,:,:));
trial2 = squeeze(vort2D(frame,:,:));
trial3 = trial .* (cores > 0) .* sign(trial2); % use cores identified by findcores_new previously
if(nargin==1 && showplot == 1)
colormap(hsvmap(30,5,[1 0 0],[0 0 1]));
pcolor(x*u_tau/nu, y*u_tau/nu, squeeze(u_norm(frame,:,:))); % the squeeze is to squeeze down a 3d matrix (where one of the dimensions is one) to a 2d matrix
caxis([-3*piv_ustd 3*piv_ustd]) % i am setting the colorbar based on the rms of the signal
%contour(x*u_tau/nu,y*u_tau/nu,trial3, 21);
%meshc(x*u_tau/nu,y*u_tau/nu,trial3);
%caxis(caxlim) % i am setting the colorbar based on the rms of the signal
%surfc(x*u_tau/nu,y*u_tau/nu,trial3);
cb = colorbar('vert'); % display the colorbar
shading interp; % interpolated shading
lighting phong; % nice spline fit to the shading
title(['Image frame # ',num2str(frame),'at time ', num2str(frame/vecimfreq, '%.5f'), ' s'])
xlabel('x^{+} ->');
ylabel('y^{+} ->');
zlabel('Swirl Values ->');
daspect([1 1 1]); % set aspect ratio to one (equal x and y scale)
hold on;
% contour(x*u_tau/nu,y*u_tau/nu,squeeze(sw2D(frame,:,:)).*(squeeze(sw2D(frame,:,:))>2*stdsw2d),12);
% contour(x*u_tau/nu,y*u_tau/nu,(squeeze(vort2D(frame,:,:))).*((squeeze(vort2D(frame,:,:)))>1.2*stdvort2d), 12 );
% contour(x*u_tau/nu,y*u_tau/nu,-1*(squeeze(vort2D(frame,:,:))).*(-1*(squeeze(vort2D(frame,:,:)))>1.2*stdvort2d), 12 );
% contour(x*u_tau/nu,y*u_tau/nu,abs(squeeze(vort2D(frame,:,:))).*(abs(squeeze(vort2D(frame,:,:)))>1.5*stdvort2d), 12 );
% contour(x*u_tau/nu,y*u_tau/nu,squeeze(sw2D(frame,:,:)).*(squeeze(sw2D(frame,:,:))),12);
quiver(x*u_tau/nu,y*u_tau/nu,squeeze(u_norm(frame,:,:)), squeeze(v(frame,:,:)),'k');
[c h]=contour(x*u_tau/nu,y*u_tau/nu,trial3, 2);
set(h,'linewidth',3);
hold off;
end
otherwise
disp (' ERROR: Insufficient arguments.');
help calcsw;
end
