function [t_vect, syn_vect] = syn_act(t_max,t_step,rho,nsyn)
%
% [t_vect,syn_vect] = syn_act(t_max,t_step,rho,nsyn)
%
%syn_act generates a train of synaptic activations with rate rho
%(activation per second per synapse) for nsyn independent synapses 
%(total rate: nsyn*rho in units of activation/sec)

%time vector
t_vect = 0:t_step:t_max; %in msec

%prepare space for synaptic events
syn_vect = zeros(size(t_vect));

rho_tot = nsyn*rho; %total rate in activations/sec
rho_tot_ms = rho_tot/1000; %total rate in activations/msec

%number of expected events in the interval 0;t_max
n_act = ceil(t_max*rho_tot_ms); 

mu = 1/rho_tot_ms; %interevent interval in msec

%generate twice as many events as 
%expected to be on the safe side
r = exprnd(mu,1,2*n_act); 
                          

%compute spike times from spike intervals
r2 = cumsum(r); 

r2 = r2(find(r2<=t_max)); %keep only the spikes times <= t_max

%convert to indices in the syn_vect vector since
%a spike time of 0 gives an index of 0 we add
%1 in the end
r2_ind = round(r2/t_step)+1; 
                             
%sets the appropriate indices to 1     
syn_vect(r2_ind) = 1;