function [ e1, e2 ] = decoyprotocol(lambda,  eta, y0, n)
% Decoy state protocol simulation
% Author: Michel Barbeau
% Version: March 7, 20176
% lambda = parameter of Poisson distribution
% eta = theoretical quantum efficiency
% y0 = vacuum pulse detection (dark count rate)
% n = number of generated k-photon random pulses
% e1 = quantum efficiency under normal operation
% e2 = quantum efficiency under attack
%
% generate n random k-photon pulses
k=poissrnd(lambda,1,n);
% generate n random numbers in the interval (0,1)
r=rand(1,n);
%%% Normal operation
% detection probability
y=y0+1-(1-eta).^k;
% measure the gain
G=sum(r<=y)/n;
% using measured gain, determine quantum efficiency
e1=-log(1+y0-G)/lambda;
%%% Under attack
% find indices of single-photon pulses
i=find(k==1);
% switch single photon pulses to null
k(i)=0; 
% find indices of non null-photon pulses
i=find(k);
% split the pulse
k(i)=k(i)-1; 
% detection probability
y=y0+1-(1-eta).^k;
% measure the gain
G=sum(r<=y)/n;
% using measured gain, determine quantum efficiency
e2=-log(1+y0-G)/lambda;
end

