function stored_jitters=montecarlo_jittered_events(measured_size,measured_riseplateau,successful_sims_required,maximum_possible_simulation_loops)
%descended from montecarlo_noisy_events, which runs either single or double event waveforms and checks the longest plateau
%montecarlo_jittered_events works like this:
%input the height and the longest rise-phase plateau of a measured event
%assume the event reflects the release of two vesicles
%we want to estimate how jittered those vesicles could have been
%so we simulate events due to two vesicles having a range of size and rise/decay time constants, added at a range of jitters [these ranges are obviously a matter of some importance]
%if the summed simulated event's size is within x percent of the size of the measured event
%we add noise [an additional matter of some delicacy] to the summed event and find the longest plateau
%if the plateau's length is within x percent of that of the real event's plateau, we record the jitter between the release of the two vesicles in the simulation
%the number of successful simulations (those that meet all criteria) is a variable
%a vector containing all the jitters from succesful simulations will be returned
%nb:  a plateau is defined one of two ways, see below

%% variance reduction, antithetic variates, control variates, importance sampling, stratified sampling

%close all, tic

tauriserange=[7 25];
taudecrange1=[10 79];
taudecrange2=[80 300];
%amprange=[floor(.45*measured_size) floor(.55*measured_size)];  %in pA
amprange=[1 ceil(-1*(measured_size))];  %the individual waveform templates are initially positive-going

jitterrange=[0 50];  %in 0.1 ms steps

noiserange=[3 6];    %in pA

taper_length=10;

time_vector=[1:200]; 
stored_jitters=[];
succesful_sims=0;
loopz=0;

while succesful_sims<successful_sims_required 
    %clf
    %making the template
	tau_rise=randint(1,1,tauriserange);
	tau_decay1=randint(1,1,taudecrange1);
    tau_decay2=randint(1,1,taudecrange2);
	template1=(1-exp(-time_vector/tau_rise)).*exp(-time_vector/tau_decay1).*exp(-time_vector/tau_decay2);
    %scaling the template
    template_amplitude=randint(1,1,amprange);
    scaler=template_amplitude/max(template1);
    template1=template1*scaler;
    %tapering the beginning
	a=repmat(template1(1)/template1(2),1,taper_length); b=[length(a):-1:1]; c=a.^b;d=c*template1(1); template1=[d template1];
    template1=-1*template1;
    %tapering the end 
	%ratte=(template1(end)/template1(end-1));
	%a=repmat(ratte,1,taperlength); b=[2:(length(a)+1)]; c=a.^b; d=c*template1(end); template1=[template1 d];
	%plot(template1,'r'),hold on
       
    vesicle_jitter=randint(1,1,jitterrange);
    
    %making the second template
    tau_rise=randint(1,1,tauriserange);
	tau_decay1=randint(1,1,taudecrange1);
    tau_decay2=randint(1,1,taudecrange2);
	template2=(1-exp(-time_vector/tau_rise)).*exp(-time_vector/tau_decay1).*exp(-time_vector/tau_decay2);
    %scaling the second template
    template_amplitude=randint(1,1,amprange);
    scaler=template_amplitude/max(template2);
    template2=template2*scaler;
    %tapering the beginning
	a=repmat(template2(1)/template2(2),1,taper_length); b=[length(a):-1:1]; c=a.^b;d=c*template2(1); template2=[d template2];
    template2=-1*template2;
    
    template2=[zeros(1,vesicle_jitter) template2];
    template2=template2(1:length(template1));
    %plot(template2,'b'),hold on
    
    template=template1+template2;
    %plot(template,'k'), hold on
		
	[min_val,min_ind]=min(template);    
    
    if min_val<.95*measured_size & min_val>1.05*measured_size
        %make some noise
		noise=rand(1,length(template));noise=noise*randint(1,1,noiserange);noise=noise-mean(noise); %correlated noise would be better
		%add the noise to the templates (note that this is *after* the templates are added together, maybe closer to how noise comes in to patch recordings)
		ntemplate=template+noise;
		%plot(ntemplate,'g')
        [min_val,min_ind]=min(template);
		
		[plateau_onset, longest_plateau]=find_riseplateau(1, min_ind, ntemplate);
        
		if longest_plateau>.95*(measured_riseplateau) & longest_plateau<1.05*(measured_riseplateau)
            stored_jitters=[stored_jitters vesicle_jitter];
            succesful_sims=succesful_sims+1;
            %loopz=0;
            
            %figure
            %plot(-1*ntemplate,'k', 'linewidth',3), hold on
            %abb=plateau_onset
            %bbbb=longest_plateau
            %plot([plateau_onset plateau_onset+longest_plateau], -1*[ntemplate(plateau_onset) ntemplate(plateau_onset+longest_plateau)],'r','linewidth',3)
        end
    end
    loopz=loopz+1;
    if loopz>maximum_possible_simulation_loops
        jitter_packer=repmat(inf,1,successful_sims_required-length(stored_jitters));
        stored_jitters=[stored_jitters jitter_packer]; 
        break
    end
end
% a=size(stored_jitters)
% b=size(loopz)
% c=loopz
stored_jitters=[stored_jitters loopz];

%toc