flag=0
proc err () { xpanel("**** ERROR ****") xlabel($s1) xpanel() }
strdef tstr
if (unix_mac_pc() == 3) {
   chdir(neuronhome())
   if (chdir("fctb")==0) flag=1
   if (chdir("..\fctb")==0) flag=1
   if (flag==0) {
     sprint(tstr,"ERROR: FCTB directory not found; should be in %s or %s/..\n",neuronhome(),neuronhome()) err(tstr)
   } else if (default_dll_loaded_== 0) {
     flag = nrn_load_dll("nrnmech.dll")
     if (flag==1) { default_dll_loaded_ = 1 } else {
       sprint(tstr,"ERROR loading %s/nrnmech.dll; remake with mknrndll.\n",getcwd()) err(tstr) }
   }
}

objref helpbox, helpdeck
proc help_list() {
  helpbox=new HBox(2)
  helpbox.intercept(1)
  xpanel("")
  xbutton("General","helpdeck.flip_to(0)")
  xbutton("Button actions","helpdeck.flip_to(1)")
  xbutton("Exercises","helpdeck.flip_to(2)")
  xpanel()

 helpdeck=new Deck()
helpdeck.intercept(1)
  xpanel("")
    xlabel("Model for Fig. 11-10 shows how positive and negative stimuli add together.")
    xlabel("")
  xpanel()
  xpanel("")
    xlabel("Graph: remember that Right button->View->pull-right->View=plot can")
    xlabel("  be used to bring the full graph into view")
    xlabel("Run/Stop: run simulation")
    xlabel("Erase: erase graph")
    xlabel("tstop: duration for running simulation in milliseconds")
    xlabel("Change parameters: menu to choose set of parameters to alter")
    xlabel("  Membrane: allows change of membrane capacitance (cm) and conductance (g)")
    xlabel("  Stimulation #1 changes ")
    xlabel("    delay in millisec (del)")
    xlabel("    duration in millisec (dur)")
    xlabel("    maximum current amplitude (amp)")
    xlabel("    i -- not setable -- just gives current current value ")
    xlabel("  Stimulation #2 has parameters controlling the second current injection")
    xlabel("")
  xpanel()
  xpanel("")
  xpanel()
helpdeck.intercept(0)
  helpdeck.map()
  helpdeck.flip_to(0)
  helpbox.intercept(0)
  helpbox.map()
}
// Created 07/24/09 17:50:27 by "/usr/site/scripts/loadfiles -q -a BIO_neg.hoc"
//================================================================
// INSERTED BIO_neg.hoc
/* =Id=  BIO_neg.hoc,v 1.6 2002/10/06 01:42:03 billl Exp  */

// load_file("BIO_neg.hoc")
//================================================================
// INSERTED string2.hoc
// =Id=  string2.hoc,v 1.1 2002/09/07 12:55:06 billl Exp 

//*String2 template
begintemplate String2
public s,t
strdef s,t
proc init() {
    if (numarg() == 1) { s=$s1 }
    if (numarg() == 2) { s=$s1 t=$s2 }
  }
endtemplate String2
// END string2.hoc
//================================================================
//================================================================
// INSERTED bkutils.hoc
// =Id=  bkutils.hoc,v 1.73 2004/04/30 15:35:32 billl Exp 

// load_file("stdgui.hoc")
// load_file("setup.hoc")
load_file("stdrun.hoc")
if (name_declared("install_matrix")) execute("install_matrix()") 
if (name_declared("install_vecst")) execute("install_vecst()")

//* objects
objref  g[4], stim[2], nc, gp, sns, ind, vec0, vec[2], rdm[2], XO, YO, cvode
objref st,wt,scr,smat,box[4],insr,nil,sref,tmpfile,tmpobj,tmpvec,tmplist,sfunc,deck
strdef mesg,tstr,tstr2,temp_string_,filename,section
tmpfile = new File()
tmplist = new List()
fchooser_flag = 0
graph_flag = 1
double x[4],y[4]
rdm = new Random()
cvode = new CVode()
sfunc = new StringFunctions()
{ ind=new Vector() vec[0]=ind.c vec[1]=ind.c vec0=ind.c}

//* iterators and templates
//** string iterator
iterator case() { local i
    i1 = 0
    for i = 2, numarg() {
        $&1 = $i
        iterator_statement
        i1+=1
      }
  }

iterator scase() { local i
    i1 = 0
    for i = 1, numarg() {
        temp_string_ = $si
        iterator_statement
        i1+=1
      }
  }

// eg for scase2("a","b","c","d","e","f") print temp_string_,temp_string2_
iterator scase2() { local i
    i1 = 0
    if (numarg()%2==1) {print "ERROR: scase2 needs even number of args" return }
    for i = 1, numarg() {
        tmpobj=new String2()
        tmpobj.s=$si i+=1  tmpobj.t=$si
        iterator_statement
        i1+=1
      }
  }

// like perl chop -- removes the last character
proc chop () { sfunc.left($s1,sfunc.len($s1)-1) }

//** list iterator ltr
// usage 'for ltr(XO, tmplist) { print XO }'
iterator ltr() { local i
    if (numarg()==3) {$&3=0} else {i1 = 0}
    for i = 0, $o2.count() - 1 {
        $o1 = $o2.object(i)
        iterator_statement
        if (numarg()==3) { $&3+=1 } else { i1+=1 }
      }
    $o1 = nil
  }

//** list pairwise iterator ltrp
// usage 'for ltrp(XO, YO, list) { print XO,YO }' takes them pairwise
iterator ltrp() { local i
    if (numarg()==4) {$&4=0} else {i1 = 0}
    for (i=0;i<$o3.count()-1;i+=2) {
        $o1 = $o3.object(i) $o2 = $o3.object(i+1)
        iterator_statement
        if (numarg()==4) { $&4+=1 } else { i1+=1 }
      }
    $o1=nil $o2=nil
  }

//** vector iterator vtr
// usage 'for vtr(&x, vec) { print x }'
iterator vtr() { local i
    if (numarg()==3) {$&3=0} else {i1 = 0}
    for i = 0, $o2.size() - 1 {
        $&1 = $o2.x[i]
        iterator_statement
        if (numarg()==3) { $&3+=1 } else { i1+=1 }
      }
  }

//* vlk(vec) -- display vector slice -- more flexible than vec.printf
// vlk(vec,max)
// vlk(vec,min,max)
// prints out a segment of a vector
vlk_width=20
proc vlk () { local i,j,min,max,dual,wdh,nonl,nl
    j=dual=0 nl=1 wdh=vlk_width
    if (numarg()==1) { min=0 max=$o1.size-1 }
    if (numarg()==2) if ($2==0) { nl=min=0 max=$o1.size-1 // vlk(vec,0) flag to suppress new lines
          } else if ($2>0) { min=0 max=$2-1 } else { min=$o1.size+$2 max=$o1.size-1 }
    if (numarg()==3) if ($3>-1) { min=$2 max=$3 } else { min=0 max=$o1.size-1 dual=1 }
    if (numarg()==4) { min=$3 max=$4 dual=1 }
    if (min<0) min=0
    if (max>$o1.size-1) max=$o1.size-1
    if (dual) if (max>$o2.size-1) max=$o2.size-1
    for i=min,max { 
        if (dual) printf("%g:%g ",$o1.x[i],$o2.x[i]) else printf("%g ",$o1.x[i]) 
        if ((j=j+1)%vlk_width==0 && nl) { print "" }
      }
    if (nl) print ""
  }

//** savevec([list,]vec1[,vec2,...]) add vector onto veclist or other list if given as 1st arg
objref veclist
veclist = new List()
proc savevec () { local i,flag,beg
    if (isobj($o1,"List")) beg=2 else beg=1
    for i=beg, numarg() { 
        tmpvec = new Vector($oi.size)
        tmpvec.copy($oi)
        if (beg==2) $o1.append(tmpvec) else veclist.append(tmpvec)
        tmpvec = nil
      }
  }

//* graphics; symbols and colors
{symnum = 7 colnum = 9}
objectvar cl[colnum], sym[symnum]
for ii=0,colnum-1 cl[ii]=new String()
{ cl[0].s ="white" cl[1].s ="black"  cl[2].s ="red"  cl[3].s ="blue"  cl[4].s ="green"  
  cl[5].s ="orange" cl[6].s ="brown"  cl[7].s ="violet"  cl[8].s ="yellow"  }
for ii=0,symnum-1 sym[ii]=new String()
{ sym[0].s = "o" sym[1].s = "t" sym[2].s = "s"  
    sym[3].s = "O" sym[4].s = "T" sym[5].s = "S"  sym[6].s = "+"}

gcl=0  // graph color
func color () {
    gcl = (gcl+1)%colnum
    if (gcl==0) gcl=1 // throw out white
    if (numarg()==1) $o1.color(gcl)
    return gcl
  }

//** qhp() basic buttons
proc qhp () {
    if (numarg()==1) xpanel("CONTROLS",1)
    xbutton("Quit","quit()")
    xbutton("Help","help_list()")
    xbutton("Print panel","pwman_place(50,50)")  
    xbutton("Reset","reset()")
    if (numarg()==1) xvarlabel(mesg)
    if (numarg()==1) xpanel(0,0)
    mesg="" 
  }

//** allgraphs()
proc remgrs () { local ii
    for ii=0,3 graphList[ii].remove_all()
    objref graphItem
  }
proc allgraphs () { for ii=0,3 for ltr(XO,graphList[ii]) XO.exec_menu($s1) }
proc vp() { allgraphs("View = plot") }
proc allgrop () { 
    for ii=0,3 for ltr(XO,graphList[ii]) {
        sprint(tstr,"%s.%s",XO,$s1)
        execute(tstr)
      }
  }

//** grrtsize() use view=plot and then pad a little
proc grrtsize () { local h,w,frac
    if (numarg()>=1) tmpobj=$o1 else tmpobj=graphItem
    if (numarg()>=2) frac = $2 else frac=.05
    tmpobj.exec_menu("View = plot")
    tmpobj.size(&x)
    w=frac*(x[1]-x[0])  h=frac*(x[3]-x[2])
    x[0]-=2*w x[1]+=w x[2]-=4*h x[3]+=h // need extra padding on bottom
    tmpobj.size(x[0],x[1],x[2],x[3])
  }

//** isobj(o1,s2) checks whether object $o1 is of type $s2
func isobj () {
    sprint(temp_string_,"%s",$o1)
    if (sfunc.substr(temp_string_,$s2)==0) {
        return 1
      } else {
        return 0
      }
  }

// cbin(min,max): binary colormap
proc cbin () {
    $o1.colormap(2)
    $o1.colormap(0, 255, 0, 0)
    $o1.colormap(1, 255, 255, 0)
    $o1.scale($2, $3)
  }

// my version of newPlotV() leaves off the label
proc newplotv () {
    newplot()
    graphItem.addvar("","v(.5)",1,3)
  }

proc redraw0 () { plrdrf=0 } // stub with flag
proc plcback () { } // callback
objref oxv,oyv // coordinates: 0,1 line start; 2,3 last loc
oxv=new Vector(2)
oyv=new Vector(2)
plrdrf=1 // flag will be unset if using redraw stub
plhorz=0 // 1 for horizontal, 2 for verticle lines

proc pl () {
    if ($1==2) { // mouse down -- new line
        redraw0() // unset plrdrf if redraw0() not overwritten
        oxv.x[0]=oxv.x[1]=$2 oyv.x[0]=oyv.x[1]=$3
      } else {
        if (plrdrf) graphItem.erase_all
        oyv.line(graphItem,oxv,0,2)
        oxv.x[1]=$2 oyv.x[1]=$3
        if (plhorz==1) oyv.x[1]=oyv.x[0]
        if (plhorz==2) oxv.x[1]=oxv.x[0]
        oyv.line(graphItem,oxv,2,2)
        redraw0() // will redraw if defined
        if ($1==3) plcback(oxv,oyv)
      }
  }

// newplot([width,height,pointer]) OR newplot([pointer])
proc newplot () { local wd,ht
    graphItem = new Graph(0)
    if (numarg()==1) $o1=graphItem // make a pointer
    if (numarg()==3) $o3=graphItem
    if (numarg()>1) { wd=$1 ht=$2 } else { wd=500 ht=300 }
    graphItem.save_name("graphList[0].")
    graphList[0].append(graphItem)
    graphItem.view(0,-90,tstop,150,600,200,wd,ht)
  }

//* cbw(min,max): b/w binary colormap
proc cbw () {
    $o1.colormap(2)
    $o1.colormap(0, 0, 0, 0)
    $o1.colormap(1, 255, 255, 255)
    $o1.scale($2, $3)
  }

//* I/O
//** svrd(0/1,name,file_ext[,fflag]) save-0 and read-1
// eg svrd(1,"Image","img") will look for file.img files to read
// fflag -- return filename only -- don't call parser
proc svrd () { local ii,num,cnt,fflag
    rdflag=$1
    if (numarg()==4) fflag=$4 else fflag=0
    sprint(tstr2,"*.%s",$s3)
    if (rdflag==0) { // write
        sprint(tstr,"Write %s",$s2)
        tmpfile.chooser("w",tstr,tstr2,"WRITE","","")
        if (tmpfile.chooser()==1) {
            tmpfile.getname(filename)
            if (!fflag) parsr(rdflag)
          }
      } else if (rdflag==1) { // read
        sprint(tstr,"Read %s",$s2)
        tmpfile.chooser("r",tstr,tstr2,"READ","","")
        if (tmpfile.chooser()==1) {
            tmpfile.getname(filename)
            if (!fflag) parsr(rdflag)
          }
      }
    mesg=filename
    tmpfile.close
  }

//* xgetargs, d2b
//** range(val,min,max) -- return val only if in proper range
// eg stim.noise=range(stim.noise,0,1)
func range () {
    if ($1<$2) return $2 
    if ($1>$3) return $3
    return $1
  }

//** err()
proc err () {
    if (numarg()==1) mesg=$s1
    print mesg
  }

//** xvarstr(): display a list of strings in an xpanel using xvarlabel
proc xvarstr () { local ii,cnt
    // ivoc_style("*font", "fixed")
    if (unix_mac_pc() == 1) {
        ivoc_style("*font", "*helvetica-bold-r-normal*--14*")
      }
    xpanel($s1)
    for ltr(XO,$o2) xvarlabel(XO.s)
    xpanel()
  }

//** xgetargs(panel_name,command,arg1[,arg2,...],defaults)
// xgetargs("Random session","newrand","# of patts","patt size  ","overlap   ","5,33,7")
objref argv
argv = new Vector()
proc xgetargs () { local i,args
    args=numarg()-3  i=numarg()
    argv.resize(0)
    sprint(temp_string_,"argv.append(%s)",$si)
    execute(temp_string_)
    if (argv.size!=args) argv.resize(args)
    xpanel($s1)
    mesg=$s1
    xvarlabel(mesg)
    for i=3,numarg()-1 {
        sprint(temp_string_,"argv.x[%d]",i-3)
        xvalue($si,temp_string_)
      }
    sprint(temp_string_,"xgetexec(\"%s\",%d)",$s2,args)
    xbutton("Execute",temp_string_)
    xpanel()
  }

proc xgetexec () { local i,args
    args = $2
    if (argv.size!=args) { mesg="Error-close & relaunch panel" return }
    sprint(temp_string_,"%s(",$s1)
    for i=0,args-2 sprint(temp_string_,"%s%g,",temp_string_,argv.x[i])
    sprint(temp_string_,"%s%g)",temp_string_,argv.x[i])
    print temp_string_
    execute(temp_string_)
  }
  

//** d2b(): decimal to binary, gives a binary appearing number for showing
// eg for ii=0,100 printf("%010d\n",d2b(ii))
func d2b () { local num,ii,rem,res
    num=$1 ii=0 res=0
    while (num>0) {
        rem = int(2*(num/2-int(num/2))) // right digit
        num=int(num/2) // remove right digit
        res+=rem*10^ii
        ii+=1
      }
    return res
  }

//* matrix stuff

//** cvec(),rvec() put out col bzw row vecs in latex format
proc cvec () { local i
    printf("\\bpm ")
    for vtr(&x,$o1) if (i1<$o1.size-1) printf("%g \\\\ ",x)
    printf("%g \\epm\n",x)
  }

proc rvec () { local i
    printf("\\bpm ")
    for vtr(&x,$o1) if (i1<$o1.size-1) printf("%g & ",x)
    printf("%g \\epm\n",x)
  }

//** matp(mat,rows,cols) prints out matrix in latex format
proc matp () { local i,j,rows,cols
    rows=$2 cols=$3
    if ($o1.size!=rows*cols) {print "ERR: matp - mat.size != rows*cols" return}
    printf("\\bpm\n")
    for (i=0;i<rows;i+=1) {
        for (j=0;j<cols-1;j+=1) {
            printf("%4g & ",$o1.x[i*cols+j])
          }
        printf("%4g \\\\ \n",$o1.x[i*cols+j])
      }
    printf("\\epm\n")
  }

//** same as matp but not latex format
proc mpf () { local i,j,rows,cols
    rows=$2 cols=$3
    if ($o1.size!=rows*cols) {print "ERR: - mat.size != rows*cols" return}
    for (i=0;i<rows;i+=1) {
        for (j=0;j<cols-1;j+=1) {
            printf("%12.3g  ",$o1.x[i*cols+j])
          }
        printf("%12.04g \n",$o1.x[i*cols+j])
      }
  }

//** print out a matrix*vec multiplication
proc mpfv () { local i,j,rows,cols
    cols=$o2.size rows=$o1.size/cols
    for (i=0;i<rows;i+=1) {
        printf("|  ")
        for (j=0;j<cols;j+=1) {
            printf("%-12.03g  ",$o1.x[i*cols+j])
          }
        if (i<$o2.size) printf("  |   %12.04g \n",$o2.x[i]) else printf("  |\n")
      }
    if (cols>rows) for i=rows,cols-1 printf("%70.04g \n",$o2.x[i])
  }

//** mmlt(mdest,m1,m2,r2) matrix mult between 2 mats
// r2==c1 gives shared dimension
proc mmlt () { local i,j,r1,c1,r2,c2
    c1 = r2 = $4 // shared dimension (cols of 1st/rows of second)
    r1=$o2.size/c1 c2=$o3.size/r2
    printf("%dx%d * %dx%d = %dx%d\n",r1,c1,r2,c2,r1,c2)
    $o1.resize(r1*c2)
    vec[0].resize(c1) // length of a row of mat A
    vec[1].resize(r2) // length of a col of mat B
    for i=0,r1-1 for j=0,c2-1 {
        vec[0].mrow($o2,i,c1) vec[1].mcol($o3,j,c2)
        $o1.x[i*c2+j]=vec[0].dot(vec[1])
      }
    $o1.mprintf(r1,c2)
    dealloc(a)
  }

//* membrane and pp's
//** runbutton
running_ = 1
proc runbutton () {
    if (running_ == 0) { stoprun = 1  return }
    if (t>0 && t<tstop-0.1) continuerun(tstop) else run()
    stoprun=1 running_=0
  }

//** proc print_pp_location(PP), from doc/html/refman/nocmodl.html
proc print_pp_location() { local x //arg1 must be a point process
     x = $o1.get_loc()
     sectionname(section)
     printf("%s located at %s(%g)\n", $o1, section, x)
     pop_section()
  }

func get_pp () { local x
     x = $o1.get_loc()
     sectionname(section)
     section sref = new SectionRef()
     pop_section()
     return x
  }
// END bkutils.hoc
//================================================================
//================================================================
// INSERTED rc.hoc
// =Id=  rc.hoc,v 1.32 2002/10/11 19:49:17 billl Exp 

// load_file("rc.hoc")

//* create a cylindrical compartment 25x25 microns
create soma
access soma
L=20
diam=30
tstop=5

//** conformed units
cm = 1    // 1 uF/cm2 capacitance 
gm  = 1    // 1 mS/cm2 conductance
r =  1/gm  // 1 kOhm-cm2 resistance
Iin= 1    // 1 uA/cm2 injEcted current
dt=0.001  // 1 us, time step (ms)
v_init=0  // no battery

//*** derived measures
Vmax=Iin*r   // mV, maximal voltage (from Ohm's law)
tau=r*cm    //  ms, time constant (the time it will take to reach 63% of final value)

//*** absolute values of units
C = cm*1e-8*area(0.5) // total capacitance in uF ls

G = gm*1e-8*area(0.5)   // total conductance in mS
R = 1/G                // total resistance in kOhms
I = Iin*1e-8*area(0.5) // injected current in uA
//*** derived measures
Vmax=I*R   // maximal voltage (from Ohm's law)
tau=R*C    // time constant (the time it will take to reach 63% of final value)

//** put a passive resistor and a stimulating electrode in
insert pas
e_pas= v_init  // no battery in series with resistor
g_pas= gm*1e-3 // S/cm2
objref stim, stim1
stim = new IClamp(0.5)
stim.amp = I*1e3   // uA->nA
stim.dur = 5 // s, forever

stim1= new IClamp(0.5)
stim1.amp= I*1e3
stim1.dur= 5

//* 2 simulation methods
//** basic run
objref box[4],insr,nil
strdef name
cvode.active()

//* showmech
objectvar deck, v1, mt,mt1, deck, mlist, mm, sl
objectvar showmechseclist
strdef s
//** shminit()
proc shminit () {
    sl = $o1
    item=0
    build()
    if (numarg(2)) {
        if ($2 == 1) {
            v1.map("ShowMechanism")
          }
      }
  }

proc build () {
    mt = new MechanismType(0)
    mt.action("change_mech()")
    mt1= new MechanismType(1)
    mt1.action("change_mech()")
    deck = new Deck()
    mlist = new List()
    build_deck()
    build_box()
    mt.select(item)
    hoc_ac_=item
    change_mech()
  }

proc build_deck () { local i
    deck.intercept(1)
    xpanel("")
    xlabel("Membrane parameters:")
    xvalue("Capacitance","cm",1)
    xvalue("Conductance","g_pas",1)
    xlabel("(Conductance=1/Resistance; g=1/R)")
    xpanel()
    nrnpointmenu(stim)
    nrnpointmenu(stim1)
    deck.intercept(0)
  }

proc build_box () {
    v1 = new VBox()
    v1.save("")
    v1.intercept(1)
    // list the option available for mechtype
    xpanel("ShowMechanism", 1)
    xmenu("Change parameters")
    xbutton("Membrane","deck.flip_to(0)")
    xbutton("Stimulation #1","deck.flip_to(1)")
    xbutton("Stimulation #2","deck.flip_to(2)")
    xmenu()
    xpanel()
    deck.map("")
    v1.intercept(0)
  }
proc change_mech () {
    item = hoc_ac_
    mt.select(item)
    item = mt.selected(s)
    deck.flip_to(1)
  }
proc change_prop () {
    forsec sl if (ismembrane(s)) {
        mlist.object(item).out()
      }
  }
proc makeshow () {
    showmechseclist = new SectionList()
    forall showmechseclist.append()
    shminit(showmechseclist, 1)
  }

next=0
proc init () {
    next=0.1
    finitialize(v_init)
    fcurrent()
  }

proc resgr () {
    gcl=0
    gcl=color(graphItem)
    graphItem.size(0,tstop,0+v_init,1+v_init)
    graphItem.family(1)
    graphItem.brush(3)
    run()
  }

proc advance () {
    fadvance()
    if (t>next) {
        graphItem.flush
        doEvents()
        next+=1
      }
  }

proc runbutton () {
    if (runstate_ == 0) {
        stoprun=1
        return
      }
    gcl=color()
    graphItem.color(gcl)
    execute1("run()")
    stoprun=0
    runstate_ = 0
  }

proc boxit () { 
    runstate_=0
    mapped_nrnmainmenu_ = 0
    name=$s1
    for ii=0,3 box[ii]=nil
    box[0] = new VBox() 
    box[1] = new HBox() 
    box[2] = new HBox() 
    box[3] = new VBox()
    box[1].intercept(1) // open HBox
    box[0].intercept(1) // open VBox#1
    qhp(1)
    newplotv()
    graphItem.size(0,tstop,0+v_init,2+v_init)
    box[0].intercept(0)  // close VBox#1
    box[0].map("")
    box[3].intercept(1) // open VBox#2
    xpanel("")
    xstatebutton("Run/Stop",&runstate_,"runbutton()")
    xbutton("Erase","graphItem.erase")
    xvalue("tstop", "tstop", 1, "tstop_changed()", 0, 1)
    xpanel()
    makeshow()
    box[3].intercept(0) // close VBox#2
    box[3].map("")
    box[1].intercept(0) // close HBox
    box[1].map($s1)
  }

runstate_=1

//* integrations and explicit solution
proc printeqns  () {
    printf("Explicit integration: v = v*(1-dt*gm/cm) + Iin*(dt/cm) ->\n\tv = %g*v + %g\n",(1-dt*gm/cm),Iin*(dt/cm))
    printf("Implicit integration:  v = (v + Iin*dt/cm)/(1+dt*g/cm) ->\n\tv = %g*v + %g\n",1/(1+dt*gm/cm),Iin*dt/cm/(1+dt*gm/cm))
  }

proc integ  () {   // explicit
    if (t==0) v=0
    v = v*(1-dt*gm/cm) + Iin*(dt/cm)
    t += dt
  }

proc integ2 () {  // implicit
    if (t==0) v=0
    v = (v + Iin*dt/cm)/(1+dt*gm/cm)
    t += dt
  }

func soln () {
    if (numarg()==1) t=$1
    v = (v_init - I/G) * exp(-t*G/C) + I/G
    return v
  }

proc accuracy () { local imp,exp,n,sol
    tstop=$1 dt=$2 n=int(tstop/dt)
    t=0 v=0  
    for ii=0,n integ2()
    imp=v
    t=0 v=0 
    for ii=0,n integ()
    exp=v
    sol = soln(tstop) 
    print tstop,imp,exp,sol," Diffs: ",abs(sol-imp),abs(sol-exp)
  }

//* bring up graphers
dt=0.001 // reset time step (ms) in case cvode used

objref grpr

// grpr.info( indep, xexpr, generate, g.size($4, $5, $6, $7), steps, x1, x2 )
proc newgrapher () { local flag
    // grpr=new Grapher()  // need to rewrite
    dt=0.001 t=0 v=0
    flag=$1
    grpr.g.erase_all()
    grpr.vbox.unmap()
    doNotify()
    if (flag==0) {
        grpr.vbox.map("Explicit integration")
        grpr.info("t","t","integ()",0,tstop,0,Vmax,tstop/dt,0,tstop)
      } else if (flag==1) {
        grpr.vbox.map("Implicit integration")
        grpr.info("t","t","integ2()",0,tstop,0,Vmax,tstop/dt,0,tstop)
      } else if (flag==2) {
        grpr.vbox.map("Analytic solution")
        grpr.info("t","t","soln()",0,tstop,0,Vmax,tstop/dt,0,tstop)
      } else print "Unknown flag"
    grpr.g.addvar("soma.v(0.5)")
    print "Press 'plot' then 'run' and compare simulations"
  }
// END rc.hoc
//================================================================

proc reset () {
    tstop=20
    cm=1
    g_pas=0.001
    stim.del=0
    stim.amp=0.01885
    stim.dur=5
    stim1.dur=5
    stim1.amp=-0.01885
    stim1.del=6
    resgr()
  }

boxit("Negative current injection")
reset()
// END BIO_neg.hoc
//================================================================