/* Execute main.hoc to run the program. The program implements the model described in Steephen & Manchanda, 2009.
 * 
 * Steephen, J. E., & Manchanda, R. (2009). Differences in biophysical properties of nucleus accumbens medium spiny neurons emerging from inactivation of inward rectifying potassium currents. J Comput Neurosci, 
 * doi:10.1007/s10827-009-0161-7
 */
 
proc loadqueue(){
	if(!bSynput) return
	UniformRand = new Random(seed)
	NormRand = new Random(seed)
	if (eAct==6) { loadQ()
	} else {
		for i = 0, 250{
			Events(Dnf)
			Events(Upf)
			if(i>83) i+=1	
		}
	}
}

proc Events() {local start, Event
	start = UniformRand.uniform(0, 1000/$1)
	Event = start + NormRand.normal(0, 1e6/(16*$1*$1))
	for j=1, 3+int(StateNo*SI*$1/1000) { //3 is for generating some extra spikes as a safety margin 
		if(Event>=0) {
			if((abs((int(int(Event/SI)/2)==(int(Event/SI)/2))-($1==Upf))&&!eAct)||(abs((Event<300)-($1==Upf))&&eAct)){
				ns.o[i].event(Event)
				if(i>83){ns.o[i+1].event(Event)}
			}
		}
		Event = start + 1000 * j/$1 + NormRand.repick()
	}
}
objref mStmType
mStmType = new Matrix(4,3)
proc loadQ(){local p localobj /*mStmType,*/ mStmPat, vEvents
	
	mStmType.x[0][0] = 300
	mStmType.x[0][1] = Dnf
				
	mStmType.x[1][0] = stim.dur
	mStmType.x[1][1] = Upf
			
	mStmType.x[2][0] = 300
	mStmType.x[2][1] = Dnf
						
	mStmType.x[3][0] = stim.dur
	mStmType.x[3][1] = Upf
		
	mStmPat = GenerateStimPattern()
				
	for p=0, 250 {
		for q=1, mStmPat.nrow()-1 {
			vEvents = events(mStmPat.x[q][1])
			FilterAdd(mStmPat, q, vEvents, ns.o[p])
			if(p>83)FilterAdd(mStmPat, q, vEvents, ns.o[p+1])				
		}
		if(p>83) p+=1		
	}	
}

proc FilterAdd() {local i,begin,end,q,dutycycle localobj mStmPat, vEvents, netcon
	mStmPat = $o1
	q =$2
	vEvents = $o3
	netcon = $o4
		
	begin = mStmPat.x[q-1][0]
	end = mStmPat.x[q][0]
	dutycycle=mStmPat.x[mStmPat.nrow()-1][0]
	if (vEvents.x[0]<0) vEvents.remove(0)
	
	while (begin<=tstop) {
		for i=0, vEvents.size()-1 {
			if((vEvents.x[i]>=begin)&&(vEvents.x[i]<end)) {
				netcon.event(vEvents.x[i])
			}
		}
		begin += dutycycle 
		end += dutycycle
	}	
}

obfunc GenerateStimPattern() {local p, q, steps, dF, nrowid localobj mStmPat
	mStmPat = new Matrix(1,2)
	nrowid=0
	for p=0, mStmType.nrow()-1 {
		steps = abs(mStmType.x[p][2]-1) // No. of actual steps is one less than that specified except for 0 for which steps is 1.
		dF = (mStmType.x[p][1]-mStmPat.x[nrowid][1])/mStmType.x[p][2]^!!mStmType.x[p][2] //avoid divide by 0
		dT =mStmType.x[p][0]/steps
		
		for q=1, steps {
			nrowid += 1
			mStmPat.resize(nrowid+1, 2)
			mStmPat.x[nrowid][0] = mStmPat.x[nrowid-1][0] + dT
			mStmPat.x[nrowid][1] = mStmPat.x[nrowid-1][1] + dF
		}
	}
	
	return mStmPat
}

obfunc events() {local start, ISI localobj vEvents
	vEvents = new Vector()
	ISI = 1000/$1
	start = UniformRand.uniform(0, ISI)
	vEvents.append(start + NormRand.normal(0, (ISI/4)^2))
	while(vEvents.max()<tstop) {vEvents.append(start + ISI * vEvents.size() + NormRand.repick())}
	return vEvents
}