close all;
clear all;


%% Poisson distribution of droplets for different labda's

LabdaArray = 0:0.001:5;                                  %Mean number of cells

kansZero = poisspdf(0,LabdaArray);                       %Poisson distribution 0 cells
kansOne = poisspdf(1,LabdaArray);                        %Poisson distribution 1 cell
kansOverig = size(length(kansZero));                     
verhoudingEenOfMeer = size(length(kansZero));       
verhoudingEenOfNul = size(length(kansZero));       

for i = 1:length(kansZero)
    kansOverig(i) = 1 - kansZero(i) - kansOne(i);       %Poisson distribution 2 or more cells
    verhoudingEenOfMeer(i) = kansOne(i)/kansOverig(i);  %Ratio droplets with 1 cell compared to droplets with more cells
    verhoudingEenOfNul(i) = kansOne(i)/kansZero(i);     %Ratio droplets with 1 cell compared to droplets with 0 cells
end

fig1=figure;
plot(LabdaArray, kansZero, LabdaArray, kansOne, LabdaArray, kansOverig)
ylabel('Poisson distribution')
xlabel('Labda')
title('Poisson distribution for multiple labdas', 'FontSize', 16)

legend('Poisson distribution 0 cells', 'Poisson distribution 1 cell', 'Poisson distribution 2 or more cells')

fig2=figure;
plot(LabdaArray, verhoudingEenOfMeer, LabdaArray, verhoudingEenOfNul)
ylabel('Ratio')
xlabel('Labda')
title('Ratio between different number of cells','FontSize', 16)

axis([min(LabdaArray) max(LabdaArray) 0 40])
legend('Ratio 1 cell / more cells', 'Ratio 1 cell / 0 cells')

%% Law of large numbers

Labda=0.3;                                              %Mean number of cells
POne=kansOne(find(LabdaArray==Labda));                  %Probability of success in a single trial
QOne=1-POne;                                            %Probability of failure in a single trial
Epsilon=0.05;                                           %Tolerance
ConfLvl=0.95;                                           %Confidence level
alfa=1-ConfLvl;

numberTrials=(POne*QOne)/((Epsilon.^2)*alfa);           %Number of trials needed

%% Poisson distribution for specific MEAN number of cells (labda).
meanCells=Labda;                                        %Mean number of cells (labda)

number=0:5;                                             %Calculate poisson distribution for different number of cells
arrayPoisson=poisspdf(number,meanCells);
figure;                                                 %Graph of Poisson distribution for differen number of cells
bar(number,arrayPoisson,0.5);
ylabel('Poisson distribution')
xlabel('Number of cell(s)')
title(sprintf('Poisson distribution cells with labda = %d', meanCells),'FontSize', 16);
for i1=1:numel(arrayPoisson)
    text(number(i1),arrayPoisson(i1),num2str(arrayPoisson(i1),'%0.2f'),...
               'HorizontalAlignment','center',...
               'VerticalAlignment','bottom')
end

kansOne = poisspdf(1,meanCells);                        %Poisson distribution 1 cell

%% Basic values for input based on experiments FACS device
usedCells=(10^9)/5;                                      %cells used for reaction DBCO PEG
volDBCO= 1.21;                                           %amount DBCO PEG (uL)
concDBCO=5;                                              %concentration DBCO PEG (milliMolair)
numberAvo=6.02214129*10^23;                              %Number of avogadro (mol^-1)

numberDroplets=10^9;

numberDBCO=((volDBCO*10^-6)*concDBCO*10^-3*numberAvo)/usedCells; %Number of DBCO molecules used per cell.

%% Poisson distribution for specific MEAN number of molecules PEG (labda).
meanMol=numberDBCO;                                     %Mean number of molecules (labda)

arrayMol=[poisscdf(meanMol,meanMol)];                   %Calculate cumulative poisson distribution for DBCO PEG
index=meanMol+100;
index1=meanMol-100;
while arrayMol(1)~=0 || arrayMol(size(arrayMol,2))~=1;
    if arrayMol(1)~=0;
        arrayMol=[poisscdf(index1,meanMol),arrayMol];
        index1=index1-100;
    end
    if arrayMol(size(arrayMol,2))~=1;
        arrayMol=[arrayMol,poisscdf(index,meanMol)];
        index=index+100;
    end;
end; 

figure;                                                 %Graph of cumulative poisson distribution.
plot(index1+100:100:index-100,arrayMol);
xlabel('Number of Molecules');
ylabel('Cumulative Poisson distribution')
title('Cumulative Poisson distribution for DBCO PEG','FontSize',16)
grid on;

index2=ceil(meanMol)-10000:500:ceil(meanMol)+10000; 
poisMol=poisspdf(index2,ceil(meanMol));                 %Array for Poisson distribution different number of cells

figure;                                                 %Graph for Poisson distribution different number of cells
bar(index2,poisMol,0.4);
axis([min(index2)-1000 max(index2)+1000 min(poisMol)-(max(poisMol)-min(poisMol)) max(poisMol)+(max(poisMol)-min(poisMol))])
grid on;
ylabel('Poisson distribution')
xlabel('Number of molecules')
title(sprintf('Poisson distribution DBCO PEG with labda = %d', meanMol),'FontSize', 16);
for i2=1:find(poisMol==max(poisMol));
    text(index2(i2),poisMol(i2),num2str(poisMol(i2)*10^5,'%0.2f'),...
               'HorizontalAlignment','center',...
               'VerticalAlignment','bottom')
end

kansMol=1-poisscdf(meanMol-1,meanMol);                  %Calculate poisson distribution desired number of molecules (or more)

%% Calculate the ratio of how many droplets have the correct form

successratio=kansOne*kansMol;

%% Calculate the volume of the droplet.

dropletDM=15;                                          %Diameter of droplet (um)

dropletVol=((4*pi*(dropletDM/2)^3)/3)*10^(-9);         %Volume of the droplet (uL)

%% Calculate the concentration needed for microfluidics if speed of solutions is the same.

% numberCell=meanCells*numberDroplets;                   %Number of cells needed
% concCell=(10^3/(dropletVol/2))*meanCells;              %Concentration of cells (cells/mL)
% 
% numberMol=meanMol*numberDroplets;                      %Number of molecules needed
% concMol=((meanMol/numberAvo)*10^3)/((dropletVol/2)*10^-6);         %Concentration molecules (milliMolair)

%% Calculate the concentration needed for microfluidics if speed of solutions is different.

vCell= 1;                                            %Speed of solution with cells (ul/min)
vDBCO= 1;                                            %Speed of solution with DBCO  (ul/min)

time=dropletVol/(vCell+vDBCO);

volCell=vCell*time;                                    %Volume of cell solution/droplet
volDBCO=vDBCO*time;                                    %Volume of DBCO solution/droplet

concCell=(10^3/volCell)*meanCells;                     %Concentration of cells (cells/mL)

numberMol=meanMol*numberDroplets;
concMol=((meanMol/numberAvo)*10^3)/((volDBCO)*10^-6);         %Concentration molecules (milliMolair)


%% Export important data to table (and .txt file but this is not

Data=[meanCells;numberDBCO;meanMol;kansOne;kansMol;concCell;concMol;successratio];
Units={'(-)';'Number of DBCO/cell';'(-)';'(-)';'(-)';'Cells/mL';'milliMolair';'(-)'};
name={'Mean number of cells (labda cells)';'DBCO molecules needed per cell';'Mean number of molecules (labda molecules)';...
    'Poisson distribution 1 cell (rate of success)';'Poisson distribution DBCO PEG (rate of success)';...
    'Concentration of cell solution';'Concentration of DBCO PEG solution';'Rate of success 1 droplet'};

tabel=table(Data,Units,'Rownames',name)
%writetable(tabel,'MicroData.txt','Delimiter','\t','WriteRowNames',true)

tabel1=table(name,Data,Units);
dataTotaal=table2cell(tabel1);
xlswrite('Data.xls',dataTotaal);