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

#include "gsl_rng.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) 
{     
    double c0, s0, c1, s1, dead, rtime, rsptime;
    double rint, prob;
    int nspikes, isp, i, j, out;
    
    const gsl_rng_type * T;
    gsl_rng * r;
    time_t seed;

    T = gsl_rng_rand48;
    r = gsl_rng_alloc (T);
  
    seed = time(NULL);
    gsl_rng_set(r, (unsigned long int)time((time_t *)NULL));
	
    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 - gsl_rng_uniform(r)*1/synout[0];
	    else { rsptime = rsptime - rtime; rtime = 0; };
	
	    for(j = 0; j < totalstim; j++)
		{
            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 > gsl_rng_uniform(r))
				{
	               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;
    };
    gsl_rng_free(r);
	return (nspikes);
}

#endif
