% ======================================================
%  Simulazione di un processo di Poisson omogeneo
%   -- implementazione efficiente --
% ======================================================
%
%   [S,N]=proc_poisson(lambda,T)
%   lambda: parametro del processo (tasso di arrivo)
%   T: istante finale
%   S: vettore contenente i tempi di arrivo -- S(1) = 0,  S(end)=T
%   N: contatore degli eventi -- N(1) = 0, N(end): # eventi in [0,T]
% 

function [S,N]=proc_poisson(lambda,T)

% in media ci saranno lambda*T arrivi
S=0;
m=ceil(1.2*lambda*T);        % proviamo con 1.2*lambda*T eventi
test=0;

while ~test
%     X=exprnd(1/lambda,1,m);
    X=-log(rand(1,m))/lambda;   % campiona m esponenziali 
    S=[S, S(end)+cumsum(X)];    % calcola tempi di arrivo
    
    if S(end)<T          % dobbiamo campionare di piu'
        m=ceil(1.2*lambda*(T-S(end)));
    else
        test=1;
    end
end
% ora bisogna troncare la lista S includendo solo gli S_i<T

n=length(find(S<T));   % il numero di eventi verificatesi e' n-1
S=[S(1:n),T];
N=[0:n-1,n-1];

end
% subplot(2,1,1),plot(S,zeros(length(S)),'.')
% se volessimo conoscere i valori del processo ad es negli istanti
% t=0:0.1:ceil(S(end));
% for i=1:length(t);N(i)=sum(S<t(i));end
% subplot(2,1,2),plot(t,N)