clear all
clc

%to do
% - check Settling velocity ... DONE
% - try to have an adaptative t ...
% - eventually do a flag & while loop ... DONE
% - verify by decreasing dfrac (particles should get carried out by the
% fluid) ... DONE
% - Model entrance profile by linearizing the creation of the velocity
% profile ... DONE
% - Verify the values and calculations
%%%%


%% Code for assessing failure in plate/tube settlers %
%+ using a Verlet-Velocity Algorithm                +%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Some useful constants
in_to_m=0.0254; %conversion factor inches to meters
m_to_in=1/in_to_m;
g=9.81;
nuw = 9e-7; %nuwater in m^2/s
%establish conditions D= Dpipe, vb = linear velocity, alpha = angle = pi/3 for plate settlers and tubes experiments

%% geometric parameters
Din=1; %diameter of the tube in inches
Pi_V='sucess';

%Length of the tube in m
L=0.463;

D=Din*in_to_m;
R = D/2;
Vup = 1e-4; % upflow velocity
alpha = pi/3; %60deg
Valpha=Vup/sin(alpha);

%vratio = 1.5 for plates and 2 for tubes.
vratio=2;
Va=vratio*Valpha;

L_e=Re(D,Va/2,nuw)*0.06*D; %entrance region Length - Only for tubes !

%% particles = #part, Vset establishes array of velocities for each particle, Q = vol flow rate, 
particles = 10;

d0=1e-6;
dflocsp = logspace(log10(d0),-3,particles); %particles diameters
dfrac = 2.3;

rhof= 2.624*1e3; %kg/m^3
rhow= 1.0e3;

rhos= (rhof-rhow);

phi = 1;% 45/24; %shape factor of particles

mfloc =1/(18*phi*nuw.*rhow).*rhos.*d0.^(3-dfrac)./dflocsp.^(1-dfrac);
Vset = mfloc .* g./rhow;

%build matrix of settling velocities for inclined tube
VsetR = Vset .* sin(alpha);
VsetZ = Vset .* cos(alpha);
% Vup=vratio*2*Q/(pi*R^2.*sin(alpha));

% t=time step length, iterations = total length
% t = 1e-2;
iterations = 4e5;

iter=zeros(particles);
temp=min(VsetR./dflocsp,VsetZ./dflocsp).*10;

%establish storage parameters (time,particle,[r,z])
%initial conditions for all particles
position = zeros(iterations,particles,2);

% %% Determined position +R
startp=R-dflocsp/2;
position(1,1:particles,1) = startp;

%% Determined position 0
% startp=0;
% position(1,1:particles,1) = startp;

%% Random position for all particles
% position(1,1:particles,1)=  D*rand(particles,1)-R;

position(1,1:particles,2)=0;


CapturedTF=zeros(particles,1); %check if particles have been captured or not
%keep track is zero unless the particle has been captured
DoneTF=CapturedTF;
i_stop=ones(particles,1); %record when the DoneTF condition has been changed

%% floc rollup velocity
%Primitive test: 
Vrol = Vfluid (R-dflocsp, R, L/2, Va, L_e);
                                             
RollUpTF= VsetR < Vrol;
disp(['    VsetR    ' 'Vrol    ' '   RollupTF ' ' dflocsp (µm)'])
disp([VsetR' Vrol' RollUpTF' dflocsp'.*1e6]);

%% step through all times
for p = 1:particles
    
    if dflocsp(p) >= D %particles are not going to enter this tube
       DoneTF(p)=1;
       CapturedTF(p)=1;
       i_stop(p)=1;
    end
    
    % for i = 2:iterations or stop condition
    i=1;
    %for each time step, step through all particles
    
     %establish r, use it to calculate local velocity
     while ~DoneTF(p) && i<iterations
        i=i+1;
        r = position(i-1,p,1);
        z = position(i-1,p,2);
        t=temp(p);
        
        if r > -R+dflocsp(p)/2
            if  z < L
                  %velocity profile in the tube
                  Vf = Vfluid (r, R, z, Va, L_e);
                  
                  %velocity-verlet algorithm to calculate 
                  %new positions linearly
                  position(i,p,1) = r - VsetR(p) * t;
                  position(i,p,2) = z + (Vf-VsetZ(p)) * t;
                  
            else %  (position(i-1,p,2) >= L)
                  position(i,p,1) = r;
                  position(i,p,2) = L+1;
                  CapturedTF(p)=0;
                  DoneTF(p)=1;
            end
        else

            %force particles to stay at -R;
            position(i,p,1)=-R+dflocsp(p)/2;
            Vfedge = 0;% Vfluid (R-dflocsp(p)/2, R, z, Va, L_e);
            position(i,p,2) = z + (Vfedge-VsetZ(p)) * t ;
            
            if position(i,p,2) < z
                if position(i,p,2) < 0
                i=i+1;
                position(i,p,1)=r;
                position(i,p,2)=0;
                DoneTF(p)=1; CapturedTF(p)=1;
                else
                    CapturedTF(p)=0;
                end
            elseif position(i,p,2) >= L
                DoneTF(p)=1; CapturedTF(p)=0;
            elseif position(i,p,2) == z
                DoneTF(p)=1; CapturedTF(p)=1;
            end
        end
        i_stop(p)=i;
     end
end
%% plot positions along pipe
elems=iterations/10000;
range=1:floor(elems):iterations;

subplot(1,3,[1 3])
hold on
axis([0 L -R.*1.01 R.*1.01])

for np=1:particles
    range=1:i_stop(np);
    plot(position(range,np,2),position(range,np,1))
    disp([num2str(np,'%2.0d') ' - ' num2str(dflocsp(np).*1e6,'%1.3f')...
        ' - Captured : ' num2str(CapturedTF(np),'%2.0f')])
end

xlabel('z (m)', 'FontSize', 14)
ylabel('r (m)', 'Rotation', 0,'FontSize', 14)


x0=zeros(3,1);
xL=repmat(L,size(x0));
xTB=linspace(0,L,3);
y0=linspace(-R,R,3);
yR=repmat(R,size(x0));
graycolor=[0.502 0.502 0.502];
plot(x0,y0,'LineStyle','--','Color',graycolor)
hold on
plot(xL,y0,'LineStyle','--','Color',graycolor)
plot(xTB,yR,'LineStyle','--','Color',graycolor)
plot(xTB,-yR,'LineStyle','--','Color',graycolor)


title(['D=' num2str(Din) ' in - L=' num2str(L) ' m - \Pi_V ='...
num2str(Pi_V)],'FontSize',16);

% axis([0 max(max(position(range,1:particles,2))) -R R])
axis([0 L -R.*1.01 R.*1.01])

%particles not captured
partpassed=ones(particles,1);
partpassed=partpassed-CapturedTF;


cmon=cumsum(partpassed);
disp(['# of particles not captured: ' num2str(cmon(length(cmon)))...
    ' /' num2str(particles) ' (Pi_V=' Pi_V ')'])