#ifndef _SPIKEGEN_H
#define _SPIKEGEN_H
#include <math.h>
#include <stdlib.h>
#include <time.h>

#define MAXSPIKES 100000

/* Pass the output of Synapse model through the Spike Generator */

int SpikeGenerator(double *synout, double binwidth, int totalstim, double nrep, double *sptime) 
{     
  
    time_t seed;
    double c0, s0, c1, s1, dead, rtime, rsptime, rint, prob;
    int nspikes, isp, i, j, out;
    
    seed = time(NULL);
    srand48(seed);
	
    c0      = 0.5;
    s0      = 0.001;
    c1      = 0.5;
    s1      = 0.0125;
    dead    = 0.00075;
    rtime   = 0.0; 
    rsptime = 0.0;
 
    nspikes = 0;   	/* synout[] contains prob of spike vs. time - run through it nrep times to look for spikes */
    isp = 0;
    sptime[0] = 0.;     /* start out with a spike at t=0 on 1st rep */
 
    for(i = 0; i < nrep; i++)
    {		
        if   (synout[0]>0) rsptime = rtime - drand48()*1/synout[0];
	    else { rsptime = rsptime - rtime; rtime = 0; };
	
	    for(j = 0; j < totalstim; j++)
		{
		    rint,prob;
            out = 0;
            rtime  += binwidth;            /* running time */
            rint    = rtime - rsptime;
            if(rint > dead )
			{
	            prob = synout[j]*(1.0-(c0*exp(-(rint-dead)/s0) + c1*exp(-(rint-dead)/s1)))*binwidth;
	    	    if(prob > drand48())
				{
	               rsptime = rtime;
	               out = 1;
				}
			};
      
		  
		    if(out==1)
			{
				sptime[isp] = rtime;      /* leave sptime in sec */
				isp++;
				if(isp>MAXSPIKES) mexErrMsgTxt("\nToo many spikes generated: increase MAXSPIKES\n");
			}
		}
	  nspikes = isp;
    };
	return (nspikes);
};

#endif
