clear all
close all
clc

%-------------------------------------------------------------------------
% Physical parameters
%-------------------------------------------------------------------------
rho_c = 1000;  % kg/m3
mu_c = 0.001;  % Pa.s
rho_d = 680;   % kg/m3
mu_d = 5e-4;   % Pa.s
sigma = 0.03;  % N/m

% Range of diameters (mm)
dmin = 0;
dmax = 3;

% number of classes
Nc = 15;
delta = (dmax-dmin)/Nc;

% vector of lower diameter of the classes (di in mm)
di = dmin:delta:dmax-delta;
di_full = dmin:delta:dmax;    % with 1 element more (the last one)

% Initial distribution of diameters (before passage through the orifice)
Ni = [0  0  0  0  0  0  25  50  100  200  250  200  100  50  25];

%-------------------------------------------------------------------------
% Operating parameter of the process
%-------------------------------------------------------------------------
U = 
CD0 = 0.85;
beta = 0.5;
DeltaP = 

%-------------------------------------------------------------------------
% Breakup modelling
%-------------------------------------------------------------------------
% Weber number We_di for each class (vector)
We = 
% Breakup probability for each class (vector)
P = 
%-------------------------------------------------------------------------
% Number of fragments for each class (vector)
Nf = 
% Size of fragments appearing from each class (size of daughter droplets) (vector)
df = 
%-------------------------------------------------------------------------
% Search of the class number of the fragments 
% cl_f is a vector that contains the class number of fragments of size df(i)
for i=1:Nc
    look_for = df(i);
    for j=1:i
        if look_for >= di_full(j) && look_for < di_full(j+1)
        cl_f(i) = j;
        end
    end
end
%-------------------------------------------------------------------------
% Number of unbroken droplets in each class (vector)
not_bkup = 
% Number of broken droplets in each class (vector)
bkup = 
%-------------------------------------------------------------------------
% Total number of daughter droplets originating from mother droplets of higher classes 
for i=1:Nc
   
    for k=1:Nc-i
        
    end
    
end
%-------------------------------------------------------------------------
% Final state distribution
Nf_final = 
%-------------------------------------------------------------------------

% Figure
figure(1)
%bar(di, Ni, delta, 'r')