// Created by Brandon Thio 7 February 2018
// create c-fibers based upon literature defined models
// set these parameters in MATLAB:
//		tstop, len, D, segdensity, type, initialdel, del
load_file("cFiberBuilder.hoc")
load_file("balanceTigerholm.hoc")
objref fiber, f,g,c,p
objectvar stim,stim2,conduct,InFile
strdef name
proc stimulate(){
// Bring the voltages to steady state
	// dtsav = dt						// Store value of dt for later
	// dt = dt_initSS					// Set dt to large dt to speed simulation in initialization phase
	// t = -inittime
	// while (t <= -dt) {		// Advance until t = 0
		// fadvance()
	// }
	// dt = dtsav						// Reset dt to simulation dt
	// t = 0							// Set t to 0, but with initialized values
	// fcurrent()						// Sets all currents and other assigned values in accordance with current values of state variables.
// Run the simulation	
	first=0
	past=-100
	time=0
	while (t<tstop){		       	   		 // Simulation Loop starts here
		time=time+1
		sec=0
		forsec fiber.sl{
			if(sec==int(2)){
				
					//print int(fiber.nsegments/4*3)
				if(v>0 && first==1 && past<0){
					fire = 1// time of ap at 0.75 down fiber which is 3.75mm down
				}
				if(v>0 && past<0 && first==0){
					first = 1// time of ap at 0.25 down fiber which is 1.25mm down
				}
				past = v
			}
			sec=sec+1
		}
		fadvance()	      	 	// Advance simulation one time step	    
		
	}
	print time
}
proc run_nrn_script(){
//////////////////// VARIABLES /////////////////////
celsius = 37
dt_initSS = 1
inittime = 1000
rec=0
//////////////////// Stimulus //////////////////////
dura=0.1 //ms
insert97na=1
conductances97=1
///////////////// Create Fibers ///////////////////////
if(type==4){
	fiber = new cFiberBuilder(D,len,type,celsius,50/6,1,1)
	// parameters diameter, length, type, temp, segment density, insert97na flag, conductnaces97 flag
	// flags = 1 if you want them in your fiber
	// type 1:Sundt 2:Tigerholm 3:Rattay 4:Sundt
}else{
	if(type==5){
		fiber = new cFiberBuilder(D,len,4,celsius,50/6,0,0)
	}else{
		fiber = new cFiberBuilder(D,len,type,celsius)
	// parameters diameter, length, type, temp
	// type 1:Sundt 2:Tigerholm 3:Rattay 4:Sundt
	}
}
///////////////////////// STIMULATE AXON /////////////////////
v_init = fiber.v
Vrest = v_init
t=0	     			    // Set time variable (t) to zero
fiber.section[300].sec stim2 = new IClamp(0.5)		    // Define stim as a current clamp (IClamp) at position 0 on second section
fiber.section[300].sec{
	stim2.del = 0			    // Stimulus delay
	stim2.dur = dura		    	// Duration of the stimulus
	stim2.amp = start * 1.5			    // Amplitude of the stimulus starts as zero.
}

high=start*1.5
fire=0
ct=0
low=0
print "Really start"
fiber.section[300].sec stim = new IClamp(0.5)		    // Define stim as a current clamp (IClamp) at position 0 on second section
fiber.section[300].sec{
	stim.del = initialdel+delay			    // Stimulus delay
	stim.dur = dura		    	// Duration of the stimulus
	stim.amp = mid			    // Amplitude of the stimulus starts as zero.
}
finitialize(v_init)
fcurrent()
if(type==2){
 balance()
}
fire=0
while((high-low)>0.001){
t=0
ct=ct+1
mid=(high+low)/2
fiber.section[300].sec{
	stim.amp = mid			    // Amplitude of the stimulus starts as zero.
}
fire=0
finitialize(v_init)
fcurrent()
stimulate()
if(fire==1){
	high=mid
}else{
	low=mid
}
}
//////////////////// Print Threshold to File //////////////////
sprint(name,"PairPulse/c_fiber_Thresh_%d_D_%d_Delay_%d_Type_%d.dat",D*1000,delay,type,ParticleID)
f=new File()
f.wopen(name)
f.printf("%f\n",high)
f.close()
}