COMMENT
This mod file is based on netstim.mod in NEURON
modified by Yi Zhou to generate Gaussian interval distributions
ENDCOMMENT

NEURON {
	POINT_PROCESS GaussStim
	RANGE x
	RANGE interval, number, start, factor, refrac
	RANGE N_backward, N_forward, N_normal, N_total
	RANGE rnd
	RANGE rand
}

PARAMETER {
	interval	= 10 (ms) <1e-9,1e9> 	: ideal time between spikes (ms)
	number	 	= 10000 <0,1e9>			: maximum number of spikes
	start		= 10 (ms)				: delay before first spiketime
	factor  	= 4  					: fraction of mean to compute std (= mean / factor)
	refrac		= 0.5 <0,3>    			: absolute refractory period (upper limit of firing = 1 kHz)
										: refrac >= dt, otherwise x can't be reset to 0 
}

ASSIGNED {
	x			: if x > 0, cell is in refractory period
	event (ms)	: spiketime of next event in absolute time
	on
	end (ms)
	m (ms)		: mean of Gaussian
	diff (ms)
	N_forward	: swap spike whose value exceed forwardly one interval 
	N_backward	: swap spike whose value exceed backwardly one interval 
	N_normal
	N_total
	rnd
	rand
}

PROCEDURE seed(x) {
	set_seed(x)
}

INITIAL {
	on			= 0
	x			= 0
	diff		= 0
	m 			= interval / 2	: each T has normal distribution N(interval/2, interval/2/factor)
	N_forward	= 0
	N_backward	= 0
	N_normal	= 0
	N_total		= 0

	if (start >= 0 && number > 0) {
		:: generate the first spiketime
		event = start + invl(m) 
		
		:: make sure first spiketime is not before 0
		if (event < 0) {
			event = 0
		}

		:: skip directly to next event
		net_send(event, 3)
	}
}

PROCEDURE init_sequence(t(ms)) {
:: Turn on spike generation if n_spikes_max > 0
::	- "on" is set to 1
::	- "event" is set to "t" (input of this call)
::	- "end" is set to "interval * n_spikes_max"

	if (number > 0) {
		on = 1
		event = t
		end = t + 1e-6 + interval * (number-1)
	}
}

PROCEDURE event_time() {
:: Generate the next spiketime
::	- call invl which calls normrand()

	:: define local variables
	LOCAL diff2, T

	:: BAD diff gets overwritten in invl :facepalm:
	diff2 = diff

	:: if n_spikes_max is larger than 0
	if (number > 0) {
		:: generate a new interval until next spiketime
		T = invl(m)

		:: BAD useless
		rnd = T

		:: if new spiketime T is equal to an existing spiketime
		if (T == 0 && diff2 == 0) {
			T = T + dt
		}

		:: compute absolute time of next event (relative to 0 ms)
		event = event + diff2 + T

		:: augment total counter
		N_total = N_total + 1
 	}

	:: if next event is after the end, turn off spike generation
	if (event > end) {
		on = 0 :: BAD not the role of this function
	}
}

FUNCTION invl(mean (ms)) (ms) {
:: Returns a new interval invl
::	Warning: also sets the global variable diff

	:: define local variable
	LOCAL std

	:: compute standard deviation
	if (mean <= 0.) {
		mean = .01 (ms) 
	}
	std = mean / factor  

	:: generate new interval invl
	invl = normrand(mean, std)

 	:: make sure invl is within 0 and max interval
	if (invl >= interval) {
		:: decrease invl by n*interval
		invl = fmod(invl, interval)
		:: augment forward counter
		N_forward = N_forward + 1
	} else if (invl < 0) {
		:: increase invl by n*interval
		invl = fmod(invl, interval) + interval
		:: augment backward counter
		N_backward = N_backward + 1
	} else {
		:: augment normal counter
		N_normal = N_normal + 1
	}

	:: overwrite global variable diff
	diff = interval - invl
}

NET_RECEIVE(w) {
	:: EXTERNAL CALL
	if (flag == 0) {
		if (w > 0 && on == 0) {
			:: turn on spike generation
			init_sequence(t)
			net_send(0, 1)
		} else if (w < 0 && on == 1) {
			:: turn off spike generation
			on = 0
		}
	}

	:: MAIN CALL
	if (flag == 1 && on == 1) {
		if (x == 0) { :: time to spike (refractory period has passed)
			:: set x to random uniform number (useless, just needs to > 0)
			rand = scop_random() : [0, 1] uniform
			x = rand

			:: output event to synapse
			net_event(t)

			:: generate next spiketime
			event_time()

			:: if next spiketime is between refrac and refrac+dt
			if (event-t >= refrac && event-t <= refrac+dt) {
				:: move next spiketime up by one dt (just after refrac)
				event = event + dt
			}

			:: if spike generation is turned on
			if (on == 1) {
				:: skip directly to next event
				net_send(event - t, 1)
			}

			:: refractory period call
			net_send(refrac, 2)
		} else if (x != 0) { :: time to spike (within refractory period)
			:: output event to synapse
			net_event(t) :: BAD even still emitted even though refractory period, should be deleted

			:: generate next spiketime
			event_time() : although this spike will be ignored, still call next spike (independent cycles)

			:: if spike generation is turned on
			if (on == 1) {
				:: skip directly to next event
				net_send(event - t, 1)
			}
		}
	}

	:: REFRACTORY PERIOD CALL (END OF REFRACTORY PERIOD)
	if (flag == 2) {
		x = 0
	}

	:: INITIAL CALL
	if (flag == 3) {
		:: if spike generation is not yet turned on
		if (on == 0) {
			:: turn on spike generation
			init_sequence(t)

			:: call net_receive directly
			net_send(0, 1)
		}
	}
}