strdef cmd, isomer
objref volt_cl			  				// variable for Point-Process
objref vbox1, vbox2, vbox3, vbox4		// graphs variables
objref g_v, g_c, g_r1, g_r2		 		// graphs variables
objref i1, t1, t2					    // vectors for calculated graphs

i1 = new Vector()
t1 = new Vector()
t2 = new Vector()

create soma
soma {
	diam = 50        // micron
	L    = 63.66198  // micron, so that area = 10000 micron2
	nseg = 1 	     // dimensionless
	cm   = 1         // uF/cm2
	Ra   = 70        // ohm-cm
	
	volt_cl = new VClamp_plus(.5)
}

displ = 0 // current and voltage curves not showed (0), or 
		  // showed (1) during the simulation. Showing them will
		  // considerably slow down the simulation
incr = 0.1
st_gr = -1

proc Clamp() {

     finitialize(hold_pot)

     g_v.erase()
     g_v.beginline()
     g_c.erase()
     g_c.beginline(2,1)

	 pre_i1 = 0
	 pre_i2 = 0
	 dens   = 0    // variables for peak calculation    
     
     peak_curr1 = 0     // initializing the peak current value

     while (t<tstop) {

             dens=volt_cl.i/area(.5)*100-soma.i_cap(0.5) // clamping current in mA/cm2,
             											 // subtracted the capacitive current

             if (t>volt_cl.dur[0] && t<(volt_cl.dur[0]+10)) {
                   if (pre_i1<abs(dens)) {
                   		peak_curr1 = abs(dens)
                   	}
            	    pre_i1=abs(dens)
             }

             if (t>(volt_cl.dur[0]+volt_cl.dur[1]+$1) && t<(volt_cl.dur[0]+volt_cl.dur[1]+$1+10)) {
                   if (pre_i2<abs(dens)) {
                   		peak_curr2 = abs(dens)
                   	}
            	    pre_i2=abs(dens)
             }
             
             if (displ==1) {
				g_v.size(0, tstop, -120, 0)
				g_v.line(t, soma.v(.5))
				g_v.flush()
				g_c.size(volt_cl.dur[0]+volt_cl.dur[1]+$1, volt_cl.dur[0]+volt_cl.dur[1]+$1+20, y_min, 0.1)
				g_c.line(t, dens)
				g_c.flush()
			 }

        fadvance()
     }

     i1.append(peak_curr2/peak_curr1)
     t1.append($1)
     t2.append(log10($1))
     print peak_curr2/peak_curr1, $1

     doEvents()
}


// procedure for graph erasing

proc erase() {

     g_v.erase(0)
     g_c.erase(0)
     g_v.size(0, tstop, hold_pot-10, amp_st_test+10)
     g_c.size(volt_cl.dur[0]+volt_cl.dur[1]+$1, volt_cl.dur[0]+volt_cl.dur[1]+$1+20, y_min, 0.1)
     g_v.beginline()
     g_c.beginline(2,1)

}

proc avvio() {

    i1.resize(0)
    t1.resize(0)
    t2.resize(0)
    
    i = min_inter
    while (i <= max_inter) {
    
	    Tstop(i)
		Clamp(i)
		
		if (i>=10000 && i<=100000) {
			i = i+10000
		}
	    
		if (i>=1000 && i<10000) {
			i = i+1000
		}
		
		if (i>=100 && i<1000) {
			i = i+100
		}
		
		if (i>=10 && i<100) {
			i = i+10
		}
		
		if (i>=1.0 && i<10) {
			i = i+1
		}
		
		if (i>=0 && i<1.0) {
			i = i+0.1
		}
		
    }
    
    g_r1.erase()
    g_r1.size(0, max_inter, 0, 1)
    g_r1.begin()
    i1.line(g_r1,t1,1,1)

    for i=0, t1.size()-1 {
		g_r1.flush()
		doNotify()
    }
    
    mil=t2.max()

    g_r2.erase()
    g_r2.size(st_gr, mil, 0, 1)
    g_r2.begin()
    i1.line(g_r2,t2,1,1)

    for i=0, t2.size()-1 {
		g_r2.flush()
		doNotify()
    }
}

proc change_isomer() {
	num_to_strg(num_iso)
	sprint (cmd, "%s %s", "uninsert", isomer)
	execute(cmd)
	
	num_to_strg(nw_num_iso)
	sprint (cmd, "%s %s", "insert", isomer)
	execute(cmd)
	
	num_iso = nw_num_iso
}	
	
proc num_to_strg() {
	if ($1 == 1) {
		isomer = "na11a"
		y_min=-0.9
		hold_pot=-120
		dur_st_cond = 100
		amp_st_cond = -10
		min_inter = 1
		max_inter = 10000
		incr = 0.1
		dur_st_test = 20
		amp_st_test = -10
		st_gr = -1
	}
	if ($1 == 2) {
		isomer = "na12a"
		y_min=-0.9
		hold_pot=-120
		dur_st_cond = 100
		amp_st_cond = -10
		min_inter = 1
		max_inter = 5000
		incr = 0.1
		dur_st_test = 20
		amp_st_test = -10
		st_gr = 0
	}
	if ($1 == 3) {
		isomer = "na13a"
		y_min=-0.9
		hold_pot=-100
		dur_st_cond = 100
		amp_st_cond = -10
		min_inter = 1
		max_inter = 5000
		incr = 0.1
		dur_st_test = 20
		amp_st_test = -10
		st_gr = 0
	}
	if ($1 == 4) {
		isomer = "na14a"
		y_min=-0.8
		hold_pot=-120
		dur_st_cond = 100
		amp_st_cond = -10
		min_inter = 1
		max_inter = 1000
		incr = 0.1
		dur_st_test = 20
		amp_st_test = -10
		st_gr = 0
	}
	if ($1 == 5) {
		isomer = "na15a"
		y_min=-1.1
		hold_pot=-120
		dur_st_cond = 1000
		amp_st_cond = -20
		min_inter = 0.1
		max_inter = 5000
		incr = 0.1
		dur_st_test = 20
		amp_st_test = -20
		st_gr = -1
	}
	if ($1 == 6) {
		isomer = "na16a"
		y_min=-1.7
		hold_pot=-90
		dur_st_cond = 100
		amp_st_cond = 0
		min_inter = 0.1
		max_inter = 200
		incr = 0.1
		dur_st_test = 20
		amp_st_test = 0
		st_gr = -1
	}
	if ($1 == 7) {
		isomer = "na17a"
		y_min=-1.5
		hold_pot=-140
		dur_st_cond = 50
		amp_st_cond = -20
		min_inter = 0.1
		max_inter = 2000
		incr = 0.1
		dur_st_test = 20
		amp_st_test = -20
		st_gr = -1
	}
	if ($1 == 8) {
		isomer = "na18a"
		y_min=-1.1
		hold_pot=-70
		dur_st_cond = 100
		amp_st_cond = 0
		min_inter = 1
		max_inter = 1000
		incr = 0.1
		dur_st_test = 10
		amp_st_test = 0
		st_gr = 0
	}
	if ($1 == 9) {
		isomer = "na19a"
		y_min=-3.3
		hold_pot=-120
		dur_st_cond = 300
		amp_st_cond = -40
		min_inter = 1
		max_inter = 1000
		incr = 0.1
		dur_st_test = 50
		amp_st_test = -40
		st_gr = 0
	}
}

proc Tstop() {

     volt_cl.dur[0] = 10	  		  // ms
     volt_cl.amp[0] = hold_pot   	  // mV
     volt_cl.dur[1] = dur_st_cond     // ms
     volt_cl.amp[1] = amp_st_cond	  // mV
     volt_cl.dur[2] = $1     		  // ms
     volt_cl.amp[2] = hold_pot  	  // mV
     volt_cl.dur[3] = dur_st_test 	  // ms
     volt_cl.amp[3] = amp_st_test	  // mV
     volt_cl.dur[4] = 10     		  // ms
     volt_cl.amp[4] = hold_pot  	  // mV
     
	tstop = 0
	for ii=0, 4 {
		tstop = tstop + volt_cl.dur[ii]
	}
}

access soma
num_iso = 1
nw_num_iso = num_iso
num_to_strg(num_iso)
sprint (cmd, "%s %s", "insert", isomer)
execute(cmd)
variable_domain("nw_num_iso", 1, 9)
	
{
xpanel("RunControl", 0)
hold_pot = -120
v_init = hold_pot
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 131
xvalue("t","t", 2 )
tstop = 131
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.025
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 0.00999999
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(3,102)
}

{
xpanel("Avvio", 0)
xbutton("Avvio","avvio()")
celsius=22
xvalue("celsius")
ena=65
xvalue("ena")
xvalue("isomer (1-9)","nw_num_iso",1,"change_isomer()")
dur_st_cond = 100
xvalue("dur_st_cond", "dur_st_cond")
min_inter = 1
xvalue("min_inter", "min_inter")
xvalue("max_inter","max_inter")
xvalue("incr", "incr")
hold_pot=v_init
xvalue("holding_pot", "hold_pot")
xstatebutton("display", &displ)
xpanel(280,102)
}

// graph for soma voltage

vbox1=new VBox()
vbox1.intercept(1)
g_v=new Graph()
g_v.size(0, tstop, hold_pot-10, amp_st_test+10)
vbox1.intercept(0)
vbox1.map("Membrane voltage", 3, 530, -1, 0)

// graph for clamping current (clamp.i)

vbox2=new VBox()
vbox2.intercept(1)
g_c=new Graph()
g_c.size(111, 131, -0.9, 0.1)
vbox2.intercept(0)
vbox2.map("Soma clamp current", 530, 10, -1, 0)

// graph for normalized repriming

vbox3=new VBox()
vbox3.intercept(1)
g_r1=new Graph()
g_r1.size(0, 5000, 0, 1)
vbox3.intercept(0)
vbox3.map("Normalized repriming", 530, 320, 500, 350)

// graph for normalized log_repriming

vbox4=new VBox()
vbox4.intercept(1)
g_r2=new Graph()
g_r2.size(st_gr, 4, 0, 1)
vbox4.intercept(0)
vbox4.map("Normalized log_repriming", 1050, 320, 500, 350)