% Drop Encapsulation Simulator
% Jon F Edd
% Todd P Lagus
% Vanderbilt University
% February 14, 2012
% This MATLAB script simulates encapsulation of particles in a set number 
% of drops (N_drops) by calculating downstream positions of a particle and 
% an associated drop.  The simulation models, effectively, a very long 
% channel with a volume equal to the total volume of aqueous fluid 
% dispensed (N_drops times the average drop volume). The channel length is 
% split into sections of an effective drop length (drop volume divided by 
% the aqueous flow channel cross section). If the particle downstream 
% position intersects with a drop downstream position limits, then that 
% particle is assigned to that particular drop. The downstream particle and 
% drop positions are assigned using average particle spacing and drop 
% volumes, each with offsets determined by the respective standard 
% deviation multiplied by a random number (randn) to simulate a 
% random experiment. Polydispersity in drop size is accounted for in the 
% standard deviation of the drop volume, altering the effective lengths of 
% the drop sections, with monodisperse drops exhibiting low standard 
% deviations. The effect of ordering is included with the standard 
% deviation of particle spacing. For high degrees of ordering, the standard 
% deviation will be small, and for no ordering, the standard deviation will 
% be at its largest. For monodisperse drops with high ordering, both the 
% particles and drops will be spaced at regular intervals, indicating more 
% controlled ordering. Note that very high standard deviations may cause 
% problems if regions overlap. The simulation counts the number of 
% particles in each drop and provides a histogram matrix for the number of 
% particles, as well as a plot of the percentage of drops with 1,2,3, … Nk 
% particles. 

% Inputs:
N_drops=1e3; % number of drops (drop sample size)
channel_height=52; %height in microns
channel_width=27; %width in microns
dx_part_avg=17.5; % average particle spacing (observed, in microns)
dx_part_std=5; % standard deviation of the particle spacing (degree of ordering)
drop_vol_avg=24.39*1000; % average drop volume, cubic micron (shown here - picoliters*1000)
drop_vol_std=5*1000; % standard deviation of drop volume, cubic micron (shown here - picoliters*1000) (degree of monodispersity)

% Calculate effective drop diameter
dx_drop_avg=drop_vol_avg/(channel_height*channel_width); % average effective droplet spacing
dx_drop_std=drop_vol_std/(channel_height*channel_width); % standard deviation of the effective droplet spacing
 
Nk=zeros(100,1); % histogram counter for number of particles per drop
x_part=0; x_drop=0; m=0; % initialize x_distance down channel
h=waitbar(0,'Encapsulating particles...');
for i=1:N_drops
    x_drop=x_drop+(dx_drop_avg+dx_drop_std*randn); % generate particle x_position
    k=0; % Start the count of particles for the upcoming drop
    while x_part<x_drop %generate drop x_position
        k=k+1;
        x_part=x_part+(dx_part_avg+dx_part_std*randn);
    end
    Nk(k+1)=Nk(k+1)+1;
    m=m+1;
    if m>(N_drops/100)
        waitbar(i/N_drops,h);
        m=0;
    end
end
close(h);
 
Pk=Nk/sum(Nk);
figure; bar(0:(length(Pk)-1),Pk);