load_file("nrngui.hoc")


///////////////////////// /* model specification */ /////////////////////////
/////// topology ////////
create soma, apical, basilar, axon
connect apical(0), soma(1)
connect basilar(0), soma(0) 
connect axon(0), soma(0)

access soma

/////// geometry ////////
soma { L = 30
diam = 30 
nseg = 1
}
apical { L = 600
diam = 1
nseg = 23
}

basilar { L = 200
diam = 2 
nseg = 5
}

axon { L = 1000
diam = 1 
nseg = 37
}

/////// biophysics //////
forall { 
Ra = 100
cm = 1
}

soma {
 insert hh
}
apical { 
insert pas 
g_pas = 0.0002 
e_pas = -65
} 

basilar {
insert pas 
g_pas = 0.0002 e_pas = -65
}

axon { 
     insert hh
}

///////////////////////// /* instrumentation */ /////////////////////////
///// synaptic input ////

objref syn 
soma syn = new AlphaSynapse(0.5) 
syn.onset = 0.5 
syn.tau = 0.1 
syn.gmax = 0.05 
syn.e = 0


/// graphical display ///
objref g 
g = new Graph()
//addplot(g, 0)
g.size(0,5,-80,40) 
g.addvar("soma.v(0.5)", 1, 1, 0.6, 0.9, 2)






///////////////////////// /* simulation control */ /////////////////////////
dt = 0.025
tstop = 5
v_init = -65

proc initialize() { 
t=0
finitialize(v_init) 
fcurrent()
}

proc integrate() { 
g.begin()

while (t<tstop) {
 fadvance()
 g.plot(t)
}
 g.flush()

}

objref recv, rect, savdata, tempmatrix
recv = new Vector()
rect = new Vector()
recv.record(&soma.v(0.5))
rect.record(&t)

proc savedatafile() {
savdata = new File()
savdata.wopen("voltage.dat")

//recv.printf(savdata)

/*
for i=0,rect.size()-1 {
   savdata.printf("t=%g   v= %g\n", rect.x(i), recv.x(i))
}

*/

tempmatrix = new Matrix()
tempmatrix.resize(recv.size(),2)
tempmatrix.setcol(0, rect)
tempmatrix.setcol(1, recv)
tempmatrix.fprint(savdata, " %g")

savdata.close()
}


proc go() {
 t1 = startsw()
 initialize() 
 integrate()
 savedatafile()
 t2 = startsw()
 print "run time= ", t2-t1
}