    {load_file("nrngui.hoc")}

// toy model cell

create soma, dend, axon
access soma

connect axon(0), soma(0)
connect dend(0), soma(1)

soma {
  diam = 30
  L = 30
//      insert traub
  insert pas
  g_pas /= 20
  e_pas = -65  // same as HH resting potential
}

dend {
  diam(0:1) = 2:0.5
  L = 600
  insert pas
  g_pas /= 20
  e_pas = -65  // same as HH resting potential
}

axon {
  diam = 1
  L = 1000
  insert hh
}

// specify reasonable values for nseg--see
//   http://www.neuron.yale.edu/neuron/static/docs/d_lambda/d_lambda.html

freq = 100      // Hz, frequency at which AC length constant will be computed
d_lambda = 0.1

func lambda_f() { // currently accessed section, $1 == frequency
  return 1e5*sqrt(diam/(4*PI*$1*Ra*cm))
}

proc geom_nseg() {
  soma area(0.5) // make sure diam reflects 3d points
  forall { nseg = int((L/(d_lambda*lambda_f(freq))+0.9)/2)*2 + 1  }
}

geom_nseg()

// prepare to use Impedance class

// always a good idea to finitialize before computing impedance
v_init=-65
{finitialize(v_init)}

// demonstrate use of impedance class
objref zz
zz = new Impedance()

FREQ = 0 // Hz
WHERE = 0.5 // location in the soma that is the reference point

soma distance(0, WHERE)  // sets origin for distance calculations

proc calcZ() {
  soma zz.loc(WHERE)  // sets origin for impedance calculations
  zz.compute(FREQ, 1) // takes the impedance contributions of
                      // gating state differential equations into account
                      // but requires mechanisms to be compatible with CVODE

  print "x distance(x) input(x) input_phase(x) transfer(x) transfer_phase(x) ratio(x)"
  forall {
    print secname(), nseg
    for (x) print x, distance(x), zz.input(x), zz.input_phase(x), zz.transfer(x), zz.transfer_phase(x), zz.ratio(x)
  }
}

objref Ri_vec, x_vec

{Ri_vec = new Vector()}
{x_vec = new Vector()}
proc calcZ_g() {
  soma zz.loc(WHERE)  // sets origin for impedance calculations
  zz.compute(FREQ, 1) // takes the impedance contributions of
                      // gating state differential equations into account
                      // but requires mechanisms to be compatible with CVODE

  // print "x distance(x) input(x) input_phase(x) transfer(x) transfer_phase(x) ratio(x)"
  forall {
    print secname(),", ", nseg
    for (x) {
      // print x, distance(x), zz.input(x), zz.input_phase(x), zz.transfer(x), zz.transfer_phase(x), zz.ratio(x)
      x_vec.append(distance(x))
      Ri_vec.append(zz.input(x))
    }
  }
}

calcZ_g()

objref g
g = new Graph()

{Ri_vec.mark(g, x_vec, "+")}
{g.exec_menu("View = plot")}
{g.exec_menu("10% Zoom out")}
