% visual.m

figure(2)
u=u.*f;
v=v.*f;
quiver(x,y,u,v,10)
hold on
for i=1:rows,
   for j=1:cols,
      if (f(i,j)==0),
         plot(x(i,j),y(i,j),'o')
      end
   end
end
hold off
-----
% fcalc.m

rho=1.23;	% air density (kg/m^3)
L=.136;		% cylinder length (m)

cols=100;
rows=100;
dT=100;		% time interval (microseconds)

xLen=.03;	% FOV along x (m)
yLen=.03;	% FOV along y (m)
dx=xLen/cols;
dy=yLen/rows;

xpx=30;		% micrometers/pixel along x
ypx=30;		% micrometers/pixel along y

x=data(:,1);
x=reshape(x,cols,rows)';
x=x*xpx;		% pos along x in micrometers

y=data(:,2);
y=reshape(y,cols,rows)';
y=y*ypx;		% pos along x in micrometers

u=data(:,3);	% "vel" along x in pixels
u=reshape(u,cols,rows)';
u=u*xpx/dT;	% vel along x in m/s

v=data(:,4);	% "vel" along y in pixels
v=reshape(v,cols,rows)';
v=v*ypx/dT;	% vel along y in m/s

f=.5*(data(:,5)+1);	% flags indicating valid points
f=reshape(f,cols,rows)';

%%%%%%%%%%%%%%%%%%%%
% Top calculations %
%%%%%%%%%%%%%%%%%%%%
tmpu=u(1,:)';
tmpv=v(1,:)';
tmpf=f(1,:)';

[topu,topv]=vecintrp(tmpu,tmpv,tmpf,cols);

Ftop=dx*[topu,topv]'*topv;

%%%%%%%%%%%%%%%%%%%%%%%
% Bottom calculations %
%%%%%%%%%%%%%%%%%%%%%%%
tmpu=u(end,:)';
tmpv=v(end,:)';
tmpf=f(end,:)';

[botu,botv]=vecintrp(tmpu,tmpv,tmpf,cols);

Fbot=-dx*[botu,botv]'*botv;

%%%%%%%%%%%%%%%%%%%%%
% Left calculations %
%%%%%%%%%%%%%%%%%%%%%
tmpu=u(:,1);
tmpv=v(:,1);
tmpf=f(:,1);

[lftu,lftv]=vecintrp(tmpu,tmpv,tmpf,rows);

Flft=-dy*[lftu,lftv]'*lftu;

%%%%%%%%%%%%%%%%%%%%%%
% Right calculations %
%%%%%%%%%%%%%%%%%%%%%%
tmpu=u(:,end);
tmpv=v(:,end);
tmpf=f(:,end);

[rgtu,rgtv]=vecintrp(tmpu,tmpv,tmpf,rows);

Frgt=dy*[rgtu,rgtv]'*rgtu;

F=[Ftop,Fbot,Flft,Frgt,Ftop+Fbot+Flft+Frgt]*rho*L;
flux=[dx*sum(topv),-dx*sum(botv),-dy*sum(lftu),dy*sum(rgtu)];
flux=[flux,sum(flux)];
fx=F(1,:);
fy=F(2,:);


-----
% fcylb.m

load testcylb00000.vec -ascii
load testcylb00001.vec -ascii
load testcylb00002.vec -ascii
load testcylb00003.vec -ascii
load testcylb00004.vec -ascii
load testcylb00005.vec -ascii
load testcylb00006.vec -ascii
load testcylb00007.vec -ascii
load testcylb00008.vec -ascii
load testcylb00009.vec -ascii

data=testcylb00000;
fcalc;
Fx=fx;
Fy=fy;
Flux=flux;

data=testcylb00001;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00002;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00003;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00004;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00005;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00006;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00007;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00008;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

data=testcylb00009;
fcalc;
Fx=[Fx;fx];
Fy=[Fy;fy];
Flux=[Flux;flux];

