clear all;close all;
%
% this matlab script reads ascii-format files (in column format, possibly
% some header lines with first character "#" or "%"), then plots them. 
%
% requires the function "read_skiplines.m" to be located on the current matlab-path.
%
% a.g. kidanemariam & m. uhlmann, oct. 2014
set(0,...
   'DefaultAxesFontSize',   12,...
   'DefaultAxesFontName',  'Helvetica',...
   'DefaultTextFontSize',   12,...
   'DefaultAxesVisible',   'on',...
   'DefaultLineLineWidth',  2,...
   'DefaultLineMarkerSize', 4,...
   'DefaultAxesBox',        'on',...
   'DefaultAxesTickDir',    'out'...
   );
%
basedir = '../tmp/'; %location of ascii data files
%
iprint  = 0;  % save plots to file
odir    = './'; %directory for output of graphs
%
domain  = 'pure_fluid'; % plot statistics of pure fluid only
%
ffluid = [basedir 'vertchan_fluid_statistics_' domain '.dat'];
pfluid = [basedir 'vertchan_particle_statistics.dat'];
%
f = read_skiplines(ffluid,6);  
p = read_skiplines(pfluid,8);
%
% column format in the *fluid* file is as follows:
% 1:y 2:<u_f> 3:<u'_f*u'_f> 4:<v'_f*v'_f> 5:<w'_f*w'_f> 6:<u'_f*v'_f>
%
% column format in the *particle* file is as follows:
%  1:y-cordi. 2:<u_p>  3:<u'_p*u'_p>  4:<v'_p*v'_p> 5:<w'_p*w'_p> 
%   6:<u'_p*v'_p>  7:<\omega_z_p>   8:\phi_s
%
% NOTE: all dimensional quantities are made dimensionless with the 
% following reference scales:
% length:   channel half-width "h"
% velocity: bulk velocity "ubulk"
% AND: the solid volume fraction is normalized with its global value,
% i.e. "<\phi_s>/\Phi_s"
%-------------------------------------------------
fig(1) = figure; ax(1) = axes; hold on
%
plot(f(:,1),f(:,2),'k')
plot(p(:,1),p(:,2),'ok')
%
set(ax(1),'xlim',[0 1],'ylim',[-0.3 1.2])
set(ax(1),'xtick',0:0.2:1,'ytick',-0.2:0.2:1.2)
%
xlabel('x/h');ylabel('mean velocity');
legend('fluid','particles');
%
if iprint==1
    nam = [odir 'mean_streamwise_velocity_' domain '.eps']
    print(nam,'-depsc'); 
end
%-------------------------------------------------
fig(2) = figure; ax(2) = axes; hold on
%
plot(f(:,1),sqrt(f(:,3)),'k')
plot(f(:,1),sqrt(f(:,4)),'r')
plot(f(:,1),sqrt(f(:,5)),'b')
%
plot(p(:,1),sqrt(p(:,3)),'ok')
plot(p(:,1),sqrt(p(:,4)),'or')
plot(p(:,1),sqrt(p(:,5)),'ob')
%
set(ax(2),'xlim',[0 1],'ylim',[0 0.3])
set(ax(2),'xtick',0:0.2:1,'ytick',0:0.1:0.3)
%
xlabel('x/h');ylabel('velocity variances');
legend('<uu>_f','<vv>_f','<ww>_f','<uu>_p','<vv>_p','<ww>_p');
%
if iprint==1
    nam = [odir 'rms_velocity_fluctuations_' domain '.eps']
    print(nam,'-depsc'); 
end
%-------------------------------------------------
fig(3) = figure; ax(3) = axes; hold on
%
plot(f(:,1),f(:,6),'k')
plot(p(:,1),p(:,6),'ok')
%
set(ax(3),'xlim',[0 1],'ylim',[-0.004 0.0005])
set(ax(3),'xtick',0:0.2:1,'ytick',-0.004:0.001:0)
%
xlabel('x/h');ylabel('velocity co-variance');
legend('<uv>_f','<uv>_p');
%
if iprint==1
    nam = [odir 'Reynolds_shear_stress_' domain '.eps']
    print(nam,'-depsc'); 
end
%-------------------------------------------------
fig(4) = figure; ax(4) = axes; hold on
%
plot(p(:,1),p(:,7),'ok')
%
set(ax(4),'xlim',[0 1],'ylim',[-0.3 0.05])
set(ax(4),'xtick',0:0.2:1,'ytick',-0.3:0.1:0)
%
xlabel('x/h');ylabel('mean angular velocity around z-axis');
legend('<\omega_z>_p');
%
if iprint==1
    nam = [odir 'mean_angular_velocity_of_particles.eps']
    print(nam,'-depsc'); 
end
%-------------------------------------------------
fig(5) = figure; ax(5) = axes; hold on
%
plot(p(:,1),p(:,8),'ok')
%
set(ax(5),'xlim',[0 1],'ylim',[0 1.5])
set(ax(5),'xtick',0:0.2:1,'ytick',0:0.5:1.5)
%
xlabel('x/h');ylabel('normalized mean solid volume fraction');
legend('<\phi_s>/\Phi_s');
%
if iprint==1
    nam = [odir 'mean_solid_volume_fraction.eps']
    print(nam,'-depsc'); 
end
%-------------------------------------------------
