function [clock,EL,NA]=arrival(clock,EL,NA,lambda0)
%    the function lambda(t) is arrival rate of non-homogeneous Poisson process.
%    lambda(t) <= lambda0=constant
   % advance clock to the time of the next arrival, 
   % update number of arrivals and number in system
   clock=EL(1);
   NA=NA+1;
      %schedule a new arrival
   TA=clock+exprnd(1/lambda0,1,1);
   while (lambda0*unifrnd(0,1,1,1)>lambda(TA))
      TA=TA+exprnd(1/lambda0,1,1);
   end
   EL(1)=TA;