function output_trains=random_train(pulse_duration_presyn,pulse_duration_postsyn,samples_per_ms,halfsecond,trial,hy,cs,samples_per_trigger,...
    OVSF_presyn,OVSF_postsyn,output_channel_to_presynaptic,evoke_postsynpatic_AP_this_trial,end_postsyn_current,postsyn_AP_delay)

%this function gets information about how many IPIs, what range (if any) IPIs can vary over, and how many steps are permitted within the range (this info comes in the form of a matrix, hy)
%it then randomly generates a first train, in accord with the specifications, and creates a tally matrix
%recording which values were assigned to each interval each dimension of the tally matrix corresponds to one IPI, for example rows corresponds to IPIa, cols to IPIb
%using this tally matrix ensures that possible trains are played once before any is played twice, and similarly for subsequent repetitions
samples_per_pulse_presyn=pulse_duration_presyn*samples_per_ms
samples_per_pulse_postsyn=pulse_duration_postsyn*samples_per_ms;
postsyn_AP_delay=postsyn_AP_delay*samples_per_ms;
%-----------------
if hy
    global tally
    au=hy(1,:);                                 %the first row of hy is the number of possible values for each IPI
    %cannot_occur=[3 5 6 7 8 10 11 12 13 14 15 16 17 18]   %THIS LINE is part of the effort to reduce the possible values of IPIa
    if trial==1
        %na=cannot_occur(1)                     %THIS LINE is part of the effort to reduce the possible values of IPIa
        %while find (na(1)== cannot_occur)      %THIS LOOP is part of the effort to reduce the possible values of IPIa 
            na=rand(1,length(au));
            na=na.*au;
            na=ceil(na);                        %then na goes and provides the info needed to create the output train
            %end
        tally=zeros(au);
        if length(au)<2                         %the eventuality where au is only one element long, that is, there is only one IPI, that is, there are only two presynaptic APs
            tally=zeros(au,1);
        end
        dummest=length(au);
        nummest=prod(au);
        while dummest                           %this loop exists so that this emmfile can handle hy of any size (that is, any number of IPIs)
            nummest=nummest-(((au(dummest)-na(dummest))/au(dummest))*prod(au(1:dummest)));
            dummest=dummest-1;    
        end
        tally(nummest)=tally(nummest)+1;
        %tally(cannot_occur)=1                  %THIS LINE is part of the effort to reduce the possible values of IPIa
    end

    if trial>1
        %tally
        ca_before=min(tally);
        for iii=1:length(size(tally))-1         %this tomfoolery is necessary bc 'min' takes the min along each dimension, so need 'number of [non-singleton] dimensions of tally' times
            ca_before=min(ca_before);           %the first time is outside the for loop (if one or more of tally's dimensions are singleton, this won't change the loop's behavior)  
        end
        ca_after=ca_before-1;
        while ca_before~=ca_after
            na=rand(1,length(au));
            na=na.*au;
            na=ceil(na); 
            dummest=length(au);
            nummest=prod(au);
            while dummest                       %this loop exists so that this emmfile can handle hy of any size (that is, any number of IPIs)
                nummest=nummest-(((au(dummest)-na(dummest))/au(dummest))*prod(au(1:dummest)));
                dummest=dummest-1;    
            end
            ca_after=tally(nummest);            %escape the while loop if the randomly generated value's spot in the tally matrix has the same value as the min of the tally matrix
        end
        tally(nummest)=tally(nummest)+1;
        %if min(tally(:,1))==max(tally(:,1))    %THIS LOOP was added for when each possible IPI has happened exactly as many times as each other possible IPI -- part of the "THIS LINE" modification oct 29 2003
            %tally(cannot_occur)=tally(cannot_occur)+1   %THIS LINE is part of the effort to reduce the possible values of IPIa
        %end
    end
    %then, after the while loop, na goes and provides the info needed to create the output train
    %--------------------
    ba=na-1;                                    %correcting, for indexing (see next few lines)
    intervuls=ba.*hy(3,:);
    intervuls=intervuls+hy(2,:);
    intervuls=intervuls*samples_per_ms;         %translating the intervals from milliseconds into samples

    train_pre=zeros(samples_per_trigger,1);
    %ones=size(train_pre)
    train_pre(3001:3000+(100*samples_per_ms))=-1/(4*OVSF_presyn);   %assuming Multiclamp Commander is set at 20mv/V, after multiplication by OVSF_presyn below, %feb 17 2004 switched from 10ms to 50ms duration
    %twos=size(train_pre)                                          %this will give a hyperpolarizing pulse of -5mv for 10ms starting at ms300
    train_pre(halfsecond+1:halfsecond+samples_per_pulse_presyn)=1; %the first pulse, which occurs beginning at halfsecond and lasts (samples_per_ms*pulse_duration) samples
    for jjj=1:length(ba)                     %the length of ba = the number of IPIs -- with the initial pulse @ halfsecond, the # of IPIs = the # of remaining pulses
        fe=halfsecond+sum(intervuls(1:jjj)); %now interval is beginning of one pulse to beginning of the next (old way was end of one to beginning of the next -- 
                                             %halfsecond+samples_per_pulse+sum(intervuls(1:jjj))+((jjj-1)*samples_per_pulse)
        train_pre(fe+1:fe+samples_per_pulse_presyn)=1;
    end 
elseif isempty(hy)
    train_pre=zeros(samples_per_trigger,1);
    %ones=size(train_pre)
    train_pre(3001:3000+(100*samples_per_ms))=-1/(4*OVSF_presyn);   %assuming Multiclamp Commander is set at 20mv/V, after multiplication by OVSF_presyn below, %feb 17 2004 switched from 10ms to 50ms duration
    %twos=size(train_pre)                                          %this will give a hyperpolarizing pulse of -5mv for 10ms starting at ms300
    train_pre(halfsecond+1:halfsecond+samples_per_pulse_presyn)=1;
end                                                           
train_pre=OVSF_presyn*train_pre;
%champ=find(train_pre)

%making the pulse train for the postsynaptic neuron
%the pulse train will be in the output train from the assigned time + one sample to the assigned time + a pulse's worth of samples
%example: if the pulse is set for 600ms, sampling 10kHz, and pulse_duration 1ms, the pulse will last from samples 6001-6010, inclusive
%this keeps it consistent with the presynaptic pulses
train_post=zeros(samples_per_trigger,1);
train_post(3001:3000+(100*samples_per_ms))=-1/(4*OVSF_postsyn);
if evoke_postsynpatic_AP_this_trial
    postsynaptic_evoke_time=end_postsyn_current+postsyn_AP_delay;
    train_post(postsynaptic_evoke_time+1:postsynaptic_evoke_time+samples_per_pulse_postsyn)=1;
    %train_post(halfsecond+postsyn_AP_delay+1:halfsecond+postsyn_AP_delay+samples_per_pulse_postsyn)=1;
end

if  cs
    %this has not been created b/c a choice must be made as to whether presyn and postsyn stim times should vary independently or in concert
    %if both are present and both are random it would be more efficient to have them vary in concert
    %(i.e., s.t. no combinatorial possibility of presyn and postsyn times can occur more than one time more than any other possibility)
end
train_post=OVSF_postsyn*train_post;

%threes=size(train_pre)
%fours=size(train_post)
output_trains=[train_pre train_post];

if output_channel_to_presynaptic==2  %switching things around when the presynaptic neuron is patched through the second output channel
        output_trains=fliplr(output_trains);
end

%train_spritz is the train determining when the picospritzer is picospritzing
%train_cam is the train determining when the camera shutter is gated