% ** coded by Christian Weiland **
% ** christian.weiland@tuhh.de **
%
% This script "metric_angle_sim.m" creates data for figure 14.
% It calculates the metric d(u,v) and angle phi(u,v) for each particle
% according to equation 30 and equation 29, respectively. These values
% are averaged in selected bins over time.
%
% ADDITIONAL INFORMATION
% 1. The input files must conform to the standard output structure of the
% particle data by the option "Moving Probes" of M-Star CFD, manipulated by
% the .bash script "skriptSortieren.bash".

clear variables
close all
tic

%Radius of the reactor
R = 0.065;
%Height of the reactor
H = 0.25;

disp('Loading data...')

dataidGS    = importdata('probeGS.txt');
dataidKL    = importdata('probeKL.txt');

disp('... ready with loading data!')

GS          = dataidGS.data;
KL          = dataidKL.data;

%Particle-ID starts at 0, should be 1
GS(:,1)     = GS(:,1) + 1;
KL(:,1)     = KL(:,1) + 1;

%Sort
GS          = sortrows(GS);
KL          = sortrows(KL);

%Program puts (0,0,0) at outer point, should be at the stirrer shaft
% -> linear translation in x- und z direction
GS(:,3)     = GS(:,3) - R;
GS(:,5)     = GS(:,5) - R;
KL(:,3)     = KL(:,3) - R;
KL(:,5)     = KL(:,5) - R;

%Calculation of radial position of particles
rGS         = sqrt(GS(:,3).^2 + GS(:,5).^2);
rKL         = sqrt(KL(:,3).^2 + KL(:,5).^2);

%Radius is 14th entry
GS          = [GS rGS];
KL          = [KL rKL];

clear data*

%How many particles? Which are the time steps?
schritte     = find(GS(:,1)==1, 1, 'last');
numPart      = max(GS(:,1));
zeit = GS(1:schritte,2)-GS(1,2);

%Velocity particle
uGS               = [GS(:,7) GS(:,8) GS(:,9)];    
uKL               = [KL(:,7) KL(:,8) KL(:,9)];
    
%Velocity fluid
vGS               = [GS(:,11) GS(:,12) GS(:,13)];    
vKL               = [KL(:,11) KL(:,12) KL(:,13)];
    
disp('Calculate the angles...')

%Calculation of angle between particle and fluid velocity (not used for paper)
winkelGS     = acos( sum(uGS.*vGS,2) ./ ( sqrt(sum(uGS.*uGS,2)) .* sqrt(sum(vGS.*vGS,2)) ));
winkelKL     = acos( sum(uKL.*vKL,2) ./ ( sqrt(sum(uKL.*uKL,2)) .* sqrt(sum(vKL.*vKL,2)) ));
 
disp('Calculate the velocity metric')

% Modified sign function
vorzeichenGS         = sign(sum(uGS.*vGS,2));
vorzeichenGS(vorzeichenGS==0) = 1;
vorzeichenKL         = sign(sum(uKL.*vKL,2));
vorzeichenKL(vorzeichenKL==0) = 1;

% Calculation of metric (eq 30) for each particle at each recorded time step
differenzGS        = vorzeichenGS.*sqrt(sum((uGS - vGS).*(uGS - vGS),2));
differenzKL        = vorzeichenKL.*sqrt(sum((uKL - vKL).*(uKL - vKL),2));
   
radius = 0:0.001:R;
hoehe  = 0:0.001:H;

phiGS  = zeros(length(hoehe), length(radius));
phiKL  = zeros(length(hoehe), length(radius));
metrikGS  = zeros(length(hoehe), length(radius));
metrikKL  = zeros(length(hoehe), length(radius));
    
disp('... now for loop...')
    for r = 1:length(radius)-1
        for h = 1:length(hoehe)-1
            % Assign positions to metric and angles. Average over time.
            reihe           = find( GS(:,14) >= radius(r) & GS(:,14) < radius(r+1) & GS(:,4) >= hoehe(h) & GS(:,4) < hoehe(h+1) );
            phiGS(h, r)     = mean(winkelGS(reihe), 'omitnan');
            metrikGS(h,r)   = mean(differenzGS(reihe), 'omitnan');
            
            reihe           = find( KL(:,14) >= radius(r) & KL(:,14) < radius(r+1) & KL(:,4) >= hoehe(h) & KL(:,4) < hoehe(h+1) );
            phiKL(h, r)     = mean(winkelKL(reihe), 'omitnan');
            metrikKL(h,r)   = mean(differenzKL(reihe), 'omitnan');
        end
        disp(['Ready with: radius = ', num2str(radius(r)), ' m'])
    end
    
dauer = toc;
disp(['... ready! The script took ', num2str(dauer), ' s.'])




%CW%%#%%%%%%%%%%%%%%%%#%%%%%%%####
%%%%%%%%%%%%#%v~~~~~~\%%%#%%%%%%%%
%%%%%%%#%%%%v'        ~~~~\%%%%%%%
%%%%#%%%%%%v'dHHb      a%%%#%%%%%%
%%%%%%%#%%v'dHHHA     :%%%%%%#%%%%
%%%%%#%%%v' VHHHHaadHHb:%#%%%%%%%%
%%%%%%%#v'   `VHHHHHHHHb:%%%%%#%%%
%%%#%%%v'      `VHHHHHHH:%%%#%%%%%
%%%%%%%'        dHHHHHHH:%%#%%%%%%
%%%%#%%        dHHHHHHHH:%%%%%%%%%
%%%%%%%       dHHHHHHHHH:%%#%%%%%%
%%#%%%%       VHHHHHHHHH:%%%%%#%%%
%%%%%%#   b    HHHHHHHHV:%%%#%%#%%
%%%%%%%   Hb   HHHHHHHV'%%%%%%%%%%
%%%%#%%   HH  dHHHHHHV'%%%#%%%%%%%
%%%#%%%   VHbdHHHHHHV'#%%%%%%%%%%%
%%%%%#%    VHHHHHHHV'%%%%%#%%#%%%%
%%%#%%%%    VHHHHHHH:%%#%%#%%#%%#%
%%#%%#%%#%%#%%#%%#%%#%%#%%#%%#%%#%
%%#%%#%%#%%#%%#%%#%%#%%#%%#%%#%IMS