// Init measurements
objref vm_s, vm_d0p, vm_d0m, vm_d0d, vm_d1p, vm_d1m, vm_d1d, vm_ap, vm_ad
vm_s = new Vector(tstop/dt+2,0)
vm_d0p = new Vector(tstop/dt+2,0)
vm_d0m = new Vector(tstop/dt+2,0)
vm_d0d = new Vector(tstop/dt+2,0)
vm_d1p = new Vector(tstop/dt+2,0)
vm_d1m = new Vector(tstop/dt+2,0)
vm_d1d = new Vector(tstop/dt+2,0)
vm_ap = new Vector(tstop/dt+2,0)
vm_ad = new Vector(tstop/dt+2,0)
vm_s.record(&maincell.soma.v(0.5), dt)
vm_d0p.record(&maincell.dend[0].v(0.1), dt)
vm_d0m.record(&maincell.dend[0].v(0.5), dt)
vm_d0d.record(&maincell.dend[0].v(0.9), dt)
vm_d1p.record(&maincell.dend[1].v(0.1), dt)
vm_d1m.record(&maincell.dend[1].v(0.5), dt)
vm_d1d.record(&maincell.dend[1].v(0.9), dt)
vm_ap.record(&maincell.axon.v(0.1), dt)
vm_ad.record(&maincell.axon.v(0.5), dt)


// Variables left over from obsolete functions
// TODO find obsolete variables
bin = 0.1 	 // [ms]  binwidth of APs
cvode.active(0) // constant integration time step

access maincell.axon

objref apc, v_ISI, v_voltage
objref pre_E, pre_I_contra, pre_I_ipsi

v_ISI = new Vector()
v_voltage = new Vector(tstop/dt+2,0)

pre_E = new Vector()
pre_I_contra = new Vector()
pre_I_ipsi = new Vector()

//==== counting APs
AP_threshold = -10 	   //[mV] AP detection threshold
apc = new APCount(0.5)   // measure action potenials at the midpoint of the axon
apc.thresh = AP_threshold

//==== saving event times and membrane potentials
apc.record(v_ISI) //record event time
v_voltage.record(&maincell.axon.v(0.5), dt)

//==== monitoring input events
netcon0[0].record(pre_E)		//record input event time
netcon_inhi0[0].record(pre_I_contra)  //record input event time
netcon_inhi1[0].record(pre_I_ipsi)    //record input event time

objref v_save,v_all,v_psth
objref g_psth,x_psth
objref d, xtime

bin_psth=bin

objref v_pth, g_pth, x_pth

objref hist_ISI,xtime,g_ISI,hist
objref mean_T, sq_T, num_T
objref temp,t1, v1,v2   // operating vectors

mean_ISI=0
OVF=0    // overflow
EI=0

mean_ISI2=0
std_ISI=0
sq_sum=0
CV_ISI=0

objref v_AP, mean_T, sq_T, num_T     		// rates and spike intervals at one ITD over N_run
objref r_mean,std_mean,vs_ITD, phase_ITD    	// rates at all ITDs
objref T_ITD,std_ITD,CV_ITD   	     	// spike interval statistics at all ITDs
objref ITD,g_Rate,g_VS,data

N_run = n_runs_total
ITD_max=1 //[ms] maximum ITD tested


// Functions
proc reset_membrane() {
	finitialize(v_init)
	fcurrent()
	
	maincell.soma {
		gl_na_VCN2003
		el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
		g_rest1=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
		cm/(g_rest1*1000)
	}
	maincell.axon {
		el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
		g_rest3=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
		cm/(g_rest3*1000)
	}
	maincell.dend[0] {
		e_pas=v_init
		g_rest2=g_pas
		cm/(g_rest2*1000)
	}
	maincell.dend[1] {
		e_pas=v_init
	}

	fcurrent()
}
// proc reset_membrane() {
// 	finitialize(v_init)
// 	fcurrent()  // make sure that leak current balance the non-leak current
// 	maincell.soma {
// 		print "gl_na_VCN2003=[S/cm^2]=", gl_na_VCN2003
// 		print "soma el=", el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
// 		print "g_rest=S/cm2  ", g_rest1=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
// 		print "tau_rest=", cm/(g_rest1*1000)
// 	}
// 	maincell.axon {
// 		print "axon el=", el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003
// 		print "g_rest=S/cm2  ", g_rest3=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003
// 		print "tau_rest=", cm/(g_rest3*1000)
// 	}
// 	maincell.dend[0] {
// 		print "dendrite el=", e_pas=v_init
// 		print " dend g_rest=S/cm2  ", g_rest2=g_pas
// 		print "tau_rest=", cm/(g_rest2*1000)
// 	}
// 	maincell.dend[1] {
// 		e_pas=v_init
// 	}

// 	fcurrent()
// }

reset_membrane()

proc reset() {
	reset_membrane()

	//==== Reset all vec and var
	N_run = n_runs_total

	v_voltage=new Vector(tstop/dt+2,0)
	v_voltage.record(&maincell.axon.v(0.5),dt)   // measure the midpoint of the axon
	
	vm_s = new Vector(tstop/dt+2,0)
	vm_d0p = new Vector(tstop/dt+2,0)
	vm_d0m = new Vector(tstop/dt+2,0)
	vm_d0d = new Vector(tstop/dt+2,0)
	vm_d1p = new Vector(tstop/dt+2,0)
	vm_d1m = new Vector(tstop/dt+2,0)
	vm_d1d = new Vector(tstop/dt+2,0)
	vm_ap = new Vector(tstop/dt+2,0)
	vm_ad = new Vector(tstop/dt+2,0)
	vm_s.record(&maincell.soma.v(0.5), dt)
	vm_d0p.record(&maincell.dend[0].v(0.1), dt)
	vm_d0m.record(&maincell.dend[0].v(0.5), dt)
	vm_d0d.record(&maincell.dend[0].v(0.9), dt)
	vm_d1p.record(&maincell.dend[1].v(0.1), dt)
	vm_d1m.record(&maincell.dend[1].v(0.5), dt)
	vm_d1d.record(&maincell.dend[1].v(0.9), dt)
	vm_ap.record(&maincell.axon.v(0.1), dt)
	vm_ad.record(&maincell.axon.v(0.5), dt)

  	v_AP=new Vector(N_run,0)
  	mean_T=new Vector(N_run,0)
  	sq_T=new Vector(N_run,0)
  	num_T=new Vector(N_run,0)
}

objref fl
objref flPG
objref spkt_contra
objref spkt_ipsi
strdef flname
strdef cmd
fl = new File()
flPG = new File()

ispkt_contra = 0
ispkt_ipsi   = 0
spkt_contra  = new Vector()
spkt_ipsi    = new Vector()

proc rerun() {
	reset()

	//*********************************
	while (N_run>0) {
		v_ISI=new Vector()
		v_voltage=new Vector(tstop/dt+2,0)
		apc.record(v_ISI) 				//record event time
		v_voltage.record(&maincell.axon.v(0.5),dt)    // measure APs at the midpoint of the axon

		vm_s = new Vector(tstop/dt+2,0)
		vm_d0p = new Vector(tstop/dt+2,0)
		vm_d0m = new Vector(tstop/dt+2,0)
		vm_d0d = new Vector(tstop/dt+2,0)
		vm_d1p = new Vector(tstop/dt+2,0)
		vm_d1m = new Vector(tstop/dt+2,0)
		vm_d1d = new Vector(tstop/dt+2,0)
		vm_ap = new Vector(tstop/dt+2,0)
		vm_ad = new Vector(tstop/dt+2,0)
		vm_s.record(&maincell.soma.v(0.5), dt)
		vm_d0p.record(&maincell.dend[0].v(0.1), dt)
		vm_d0m.record(&maincell.dend[0].v(0.5), dt)
		vm_d0d.record(&maincell.dend[0].v(0.9), dt)
		vm_d1p.record(&maincell.dend[1].v(0.1), dt)
		vm_d1m.record(&maincell.dend[1].v(0.5), dt)
		vm_d1d.record(&maincell.dend[1].v(0.9), dt)
		vm_ap.record(&maincell.axon.v(0.1), dt)
		vm_ad.record(&maincell.axon.v(0.5), dt)
		
		flPG.ropen("eventtimes_contra.txt")
		spkt_contra.scanf(flPG)
		flPG.close()
		flPG.ropen("eventtimes_ipsi.txt")
		spkt_ipsi.scanf(flPG)
		flPG.close()

		for i = 0,9 {
			gen0[i].timeToNextEvt = spkt_contra.x[0] - dt
			gen1[i].timeToNextEvt = spkt_ipsi.x[0] - dt
		}

		finitialize(v_init)
		fcurrent()

		integrate()
		// get_pth()
		// get_ISI()
		v_AP.x[N_run-1] = apc.n
		N_run = N_run-1
		
		// Export to temporary files
		sprint(flname, "../data/_direct/run_%05i__vm_c.txt.tmp", iter)
		fl.wopen(flname)
		vm_s.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_d0p.txt.tmp", iter)
		fl.wopen(flname)
		vm_d0p.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_d0m.txt.tmp", iter)
		fl.wopen(flname)
		vm_d0m.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_d0d.txt.tmp", iter)
		fl.wopen(flname)
		vm_d0d.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_d1p.txt.tmp", iter)
		fl.wopen(flname)
		vm_d1p.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_d1m.txt.tmp", iter)
		fl.wopen(flname)
		vm_d1m.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_d1d.txt.tmp", iter)
		fl.wopen(flname)
		vm_d1d.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_ap.txt.tmp", iter)
		fl.wopen(flname)
		vm_ap.printf(fl)
		fl.close()
		
		sprint(flname, "../data/_direct/run_%05i__vm_ad.txt.tmp", iter)
		fl.wopen(flname)
		vm_ad.printf(fl)
		fl.close()

		// Rename files from temporary name to destined name
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_c.txt.tmp run_%05i__vm_c.txt", iter, iter)
		print "Cmd: ",cmd
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0p.txt.tmp run_%05i__vm_d0p.txt", iter, iter)
		print "Cmd: ",cmd
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0m.txt.tmp run_%05i__vm_d0m.txt", iter, iter)
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0d.txt.tmp run_%05i__vm_d0d.txt", iter, iter)
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1p.txt.tmp run_%05i__vm_d1p.txt", iter, iter)
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1m.txt.tmp run_%05i__vm_d1m.txt", iter, iter)
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1d.txt.tmp run_%05i__vm_d1d.txt", iter, iter)
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_ap.txt.tmp run_%05i__vm_ap.txt", iter, iter)
		system(cmd)
		sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_ad.txt.tmp run_%05i__vm_ad.txt", iter, iter)
		system(cmd)
	}

	N_run = n_runs_total

}

proc integrate() {
	// while the simulation has not yet ended
	while (t < tstop) {
		// CONTRA
		// if t is earlier than the last spiketime
		if (t < spkt_contra.x[spkt_contra.size()-1]) {
			// if t is right after a spiketime has been issued, issue the next one
			if (t >= spkt_contra.x[ispkt_contra-1]+dt && t < spkt_contra.x[ispkt_contra-1]+dt*2) {
				if (ispkt_contra == spkt_contra.size()-1) {
					for i = 0,9 {
						gen0[i].timeToNextEvt = 0
					}
				} else {
					for i = 0,9 {
						gen0[i].timeToNextEvt = spkt_contra.x[ispkt_contra+1] - spkt_contra.x[ispkt_contra]
					}
				}

				ispkt_contra = ispkt_contra + 1
			} else if (t >= dt && t < dt*2) { // if t is right the beginning of the simulation
				for i = 0,9 {
					gen0[i].timeToNextEvt = spkt_contra.x[ispkt_contra+1] - spkt_contra.x[ispkt_contra] - dt
				}
				ispkt_contra = ispkt_contra + 1
			}
		}

		// IPSI
		// if t is earlier than the last spiketime
		if (t < spkt_ipsi.x[spkt_ipsi.size()-1]) {
			// if t is right after a spiketime has been issued, issue the next one
			if (t >= spkt_ipsi.x[ispkt_ipsi-1]+dt && t < spkt_ipsi.x[ispkt_ipsi-1]+dt*2) {
				if (ispkt_ipsi == spkt_ipsi.size()-1) {
					for i = 0,9 {
						gen1[i].timeToNextEvt = 0
					}
				} else {
					for i = 0,9 {
						gen1[i].timeToNextEvt = spkt_ipsi.x[ispkt_ipsi+1] - spkt_ipsi.x[ispkt_ipsi]
					}
				}

				ispkt_ipsi = ispkt_ipsi + 1
			} else if (t >= dt && t < dt*2) { // if t is right the beginning of the simulation
				for i = 0,9 {
					gen1[i].timeToNextEvt = spkt_ipsi.x[ispkt_ipsi+1] - spkt_ipsi.x[ispkt_ipsi] - dt
				}
				ispkt_ipsi = ispkt_ipsi + 1
			}
		}

		fadvance() // advance solution
	}
	print "N_run=",N_run   // show run time
	print "spike number=",apc.n
	// print "-------------------------"
	//print "Number_input_E=",pre_E.size
	//print "Number_input_I_contra=",pre_I_contra.size
	//print "Number_input_I_ipsi=",pre_I_ipsi.size
	//print "*************"

	if(v_ISI.size<=2) { v_ISI=new Vector(3,0) }
	L_ISI=v_ISI.size
	//print "L_ISI=",L_ISI
	L_PTH=v_voltage.size
	//print "L_PTH=",L_PTH
}

rerun()