| core_proc(frame1,frame2,filtcth_ll,cth_ll,cth_ul,cth_st)
|
function core_proc(frame1,frame2,filtcth_ll,cth_ll,cth_ul,cth_st)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% 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: %%
%% core_proc.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 core numcores sw2D cores pix_m u_tau x vort2D nu delta_99 rho_air u v
global corethresh maxsw
global corefrm NZ
global lenconvfac areaconvfac cellarea
corefrm = zeros(2,NZ);
null_pdf_vars;
core = [];
core=struct('length',{},'width',{},'area',{},'circ',{},'ke',{},'avgvort',{},'maxsw',{},'maxswx',{},'maxswy',{});
corethresh = struct('value',{},'core',{});
disp(['Identifying cores ... ']);
m = 1;
maxsw = getswmax;
maxsw50 = 880.8177;
if(maxsw <= 0)
disp('ERROR: Zero maximum swirl value. Aborting !!!');
return;
end
for cthresh=cth_ll:cth_st:cth_ul
disp(['Running for cthresh = ', num2str(cthresh), ' ...']);
core = [];
core=struct('length',{},'width',{},'area',{},'circ',{},'ke',{},'avgvort',{},'maxsw',{},'maxswx',{},'maxswy',{});
corecount = 1;
for frame=frame1:frame2
V2 = (squeeze(u(frame,:,:)).^2 + (squeeze(v(frame,:,:))).^2);
circ = (squeeze(vort2D(frame,:,:)) * cellarea) ; % compute circulation for each cell for this frame
% temp = squeeze(sw2D(frame,:,:));
% maxsw = max(max(temp));
% findcores_new(frame, maxsw*cthresh , maxsw, 50, 1200, 4)
% findcores_new(frame, maxsw50*cthresh , maxsw, 50, 1200, 1) % maxsw50 is the upper threshold for z+ = 50
findcores_lltrial(frame, maxsw50*filtcth_ll, maxsw50*cthresh , maxsw, 4) % maxsw50 is the upper threshold for z+ = 50
for idtag=1:numcores
[xc yc] = find(cores == idtag);
core(corecount).length = (max(xc) - min(xc) + 1) * lenconvfac; % normalized length
core(corecount).width = (max(yc) - min(yc) + 1) * lenconvfac; % normalized width
core(corecount).area = size(xc,1) * areaconvfac; % normalized area
core(corecount).maxsw = max(max( squeeze(sw2D(frame,:,:)) .* (cores==idtag) ) );
[core(corecount).maxswx core(corecount).maxswy] = find(squeeze(sw2D(frame,:,:))==core(idtag).maxsw);
% above find call could be a potential problem if there two points with max swirl value
% i.e. we have a flat peak covering more than one cell
core(corecount).avgvort = sum(sum( squeeze(vort2D(frame,:,:)).*(cores == idtag) )); % Net Vorticity for this core
% core(corecount).circ = core(corecount).avgvort * (size(xc,1) * cellarea) ;
core(corecount).circ = sum(sum(circ .* (cores == idtag))) ;
core(corecount).ke = 0.5*rho_air*sum(sum( V2 .* (cores == idtag) ));
corecount = corecount+1;
end
corefrm(m,frame) = numcores;
%calc_area_pdf; %calc_maxsw_pdf;
if(rem(frame,25) == 0) % Display a message on screen for every 25th frame
disp(['Done past frame #',num2str(frame,'%.4d')]);
end
end
corethresh(m).value = cthresh;
corethresh(m).core = core;
m = m + 1;
end
return;
function null_pdf_vars()
global core_apdf core_maxswpdf
core_apdf = [];
core_maxswpdf = [];
return;
function calc_area_pdf()
global core_apdf core
core_apdf = [core_apdf [core(:).area]];
return;
function calc_maxsw_pdf()
global core_maxswpdf core
core_maxswpdf = [core_maxswpdf [core(:).maxsw]];
return;
|
|
