load_file("nrngui.hoc")


// --------------------------------------------------------------
// redefine some things in stdrun.hoc
// --------------------------------------------------------------

tstop = 5000
steps_per_ms = 50
dt = 1/steps_per_ms
secondorder = 0

// --------------------------------------------------------------
// cell geometry
// --------------------------------------------------------------

create soma 
soma { nseg = 1 L = 30 diam = 500/(PI*L) }

create soma1 
soma1 { nseg = 1 L = 30 diam = 500/(PI*L) }

//RNG
//objref rnd
//soma1 rnd = new RNG(0.5)
//rnd.seed = 365877 



// --------------------------------------------------------------
// passive & active membrane 
// --------------------------------------------------------------

ra        = 200
global_ra = ra
rm        = 400
c_m       = 0.75
v_init    = -70
celsius = 23

Ek = -90
Ena = 60
            

gna_dend = 300
gna_node = 30000
gk = 500
frac_gk_apical = 0
gk_node = 0


proc install_channels() {
  
  // passive
  forall {
    insert pas
    Ra = ra
    cm = c_m 
    g_pas = 1/rm
    e_pas = v_init
  }

  // Na+ channels

  soma insert na3
  soma.gbar_na3 = gna_dend

  forall if (ismembrane("na3")) ena = Ena     // per S&S ena = 90
/*
  soma1 insert sna
  soma1.gbar_sna = gna_dend
  soma1 { D_sna = int(gbar_sna/gama_sna) } 
*/
  forall if (ismembrane("sna")) ena = Ena         // per S&S ek = -100

}

install_channels()


// --------------------------------------------------------------
// stimulus
// --------------------------------------------------------------

objectvar st
soma  st = new IClamp(0.5)
st.amp = 0.003
st.del = 100
st.dur = 100
/*
objectvar st1
soma1  st1 = new IClamp(0.5)
st1.amp = 0.001
st1.del = 3
st1.dur = 400
*/
// --------------------------------------------------------------
// create useful graphs & panels
// --------------------------------------------------------------

access soma
nrnmainmenu()
nrncontrolmenu()
newPlotV()
//graphItem.addvar("soma1.v(.5)",7,1)


objref vv,vv1

vv = new Vector(tstop/dt+1)
vv.record(&soma.v(0.5))
tstop = 300
stdrun_quiet = 1
for (vi = -100;vi < -40;vi=vi+5) { 
  v_init = vi 
  soma.e_pas = vi
  init()
  run()
  vv1 = vv.at(((st.del+20)/dt)+1,(st.del+st.dur)/dt-2)
  print vv1.max() - soma.v(0.5)
}