def volumes1d(self, sec):
        if not isinstance(sec, nrn.Section):
            sec = sec._sec
        arc3d = [h.arc3d(i, sec=sec)
                 for i in range(int(h.n3d(sec=sec)))]
        diam3d = [h.diam3d(i, sec=sec)
                  for i in range(int(h.n3d(sec=sec)))]
        vols = numpy.zeros(sec.nseg)
        dx = sec.L / sec.nseg
        for iseg in range(sec.nseg):
            # get a list of all pts in the segment, including end points
            lo = iseg * dx
            hi = (iseg + 1) * dx
            pts = [lo] + [x for x in arc3d if lo < x < hi] + [hi]
            
            diams = numpy.interp(pts, arc3d, diam3d)
            
            # sum the volume of the constituent frusta, hollowing out by the inside
            volume = 0
            for i in range(len(pts) - 1):
                diam0h, diam1h = self._hi * diams[i : i + 2]
                diam0l, diam1l = self._lo * diams[i : i + 2]
                pt0, pt1 = pts[i : i + 2]
                volume += numpy.pi * (pt1 - pt0) / 12. * ((diam0h ** 2 + diam0h * diam1h + diam1h ** 2) - (diam0l ** 2 + diam0l * diam1l + diam1l ** 2))
            vols[iseg] = volume

        return vols

def _volumes1d(sec):
    if not isinstance(sec, nrn.Section):
        sec = sec._sec
    arc3d = [h.arc3d(i, sec=sec)
             for i in range(int(h.n3d(sec=sec)))]
    diam3d = [h.diam3d(i, sec=sec)
              for i in range(int(h.n3d(sec=sec)))]
    vols = numpy.zeros(sec.nseg)
    dx = sec.L / sec.nseg
    for iseg in range(sec.nseg):
        # get a list of all pts in the segment, including end points
        lo = iseg * dx
        hi = (iseg + 1) * dx
        pts = [lo] + [x for x in arc3d if lo < x < hi] + [hi]
        
        diams = numpy.interp(pts, arc3d, diam3d)
        
        # sum the volume of the constituent frusta
        volume = 0
        for i in range(len(pts) - 1):
            diam0, diam1 = diams[i : i + 2]
            pt0, pt1 = pts[i : i + 2]
            volume += numpy.pi * (pt1 - pt0) / 12. * (diam0 ** 2 + diam0 * diam1 + diam1 ** 2)
        vols[iseg] = volume

    return vols

def get_pos_data_short():
    """
    Get positions of all segments currently loaded in Neuron in a simple matrix.
    Section position information is not available.

    :returns: 
        Matrix (3 x nSegments) With x,y,z positions. 
    :rtype: :class:`~numpy.ndarray`

    Example:
        .. code-block:: python

            data = get_pos_data_short()
    """
    n = 0
    for sec in h.allsec():
        n += int(h.n3d())
    data = np.zeros([4, n])
    cnt = 0
    for sec in h.allsec():
        for i in xrange(int(h.n3d())):
            data[0, cnt] = h.x3d(i)
            data[1, cnt] = h.y3d(i)
            data[2, cnt] = h.z3d(i)
            data[3, cnt] = h.diam3d(i)
            cnt += 1
    return data

 def result(sec):
     if not isinstance(sec, nrn.Section):
         sec = sec._sec
     arc3d = [h.arc3d(i, sec=sec)
             for i in range(int(h.n3d(sec=sec)))]
     diam3d = [h.diam3d(i, sec=sec)
             for i in range(int(h.n3d(sec=sec)))]
     sas = numpy.zeros(sec.nseg)
     dx = sec.L / sec.nseg
     for iseg in range(sec.nseg):
         # get a list of all pts in the segment, including end points
         lo = iseg * dx
         hi = (iseg + 1) * dx
         pts = [lo] + [x for x in arc3d if lo < x < hi] + [hi]
         
         diams = numpy.interp(pts, arc3d, diam3d)
         
         # sum the surface areas of the constituent frusta
         sa = 0
         for i in range(len(pts) - 1):
             diam0, diam1 = diams[i : i + 2]
             pt0, pt1 = pts[i : i + 2]
             sa += scale * 0.5 * (diam0 + diam1) * numpy.sqrt(0.25 * (diam0 - diam1) ** 2 + (pt1 - pt0) ** 2)
         sas[iseg] = sa
     return sas

def _neighbor_areas1d(sec):
    arc3d = [h.arc3d(i, sec=sec._sec)
             for i in xrange(int(h.n3d(sec=sec._sec)))]
    diam3d = [h.diam3d(i, sec=sec._sec)
              for i in xrange(int(h.n3d(sec=sec._sec)))]
    area_pos = numpy.linspace(0, sec.L, sec.nseg + 1)
    diams = numpy.interp(area_pos, arc3d, diam3d)
    return numpy.pi * 0.25 * diams ** 2

 def result(sec):
     if not isinstance(sec, nrn.Section):
         sec = sec._sec
         arc3d = [h.arc3d(i, sec=sec) for i in xrange(int(h.n3d(sec=sec)))]
         diam3d = [h.diam3d(i, sec=sec) for i in xrange(int(h.n3d(sec=sec)))]
         area_pos = numpy.linspace(0, sec.L, sec.nseg + 1)
         diams = numpy.interp(area_pos, arc3d, diam3d)
         return scale * diams

 def neighbor_areas1d(self, sec):
     arc3d = [h.arc3d(i, sec=sec._sec)
              for i in xrange(int(h.n3d(sec=sec._sec)))]
     diam3d = [h.diam3d(i, sec=sec._sec)
               for i in xrange(int(h.n3d(sec=sec._sec)))]
     area_pos = numpy.linspace(0, sec.L, sec.nseg + 1)
     diams = numpy.interp(area_pos, arc3d, diam3d)
     if self._type == _lo_hi_shell:
         return numpy.pi * .25 * ((diams * self._hi) ** 2 - (diams * self._lo) ** 2)

def _neighbor_areas1d(sec):
    if not isinstance(sec, nrn.Section):
        sec = sec._sec
    arc3d = [h.arc3d(i, sec=sec)
             for i in range(int(h.n3d(sec=sec)))]
    diam3d = [h.diam3d(i, sec=sec)
              for i in range(int(h.n3d(sec=sec)))]
    area_pos = numpy.linspace(0, sec.L, sec.nseg + 1)
    diams = numpy.interp(area_pos, arc3d, diam3d)
    return numpy.pi * 0.25 * diams ** 2

 def rotateZ(self, theta):
     """Rotate the cell about the Z axis."""
     rot_m = numpy.array([[sin(theta), cos(theta)], [cos(theta), -sin(theta)]])
     for sec in self.all:
         for i in range(int(h.n3d())):
             xy = numpy.dot([h.x3d(i), h.y3d(i)], rot_m)
             h.pt3dchange(i, xy[0], xy[1], h.z3d(i), h.diam3d(i))

def lambda_f(section, freq):
    if h.n3d() < 2:
        return 1e5*math.sqrt(section.diam/(math.pi*4*freq*section.Ra*section.cm))
    else:
        x1 = h.arc3d(0)
        d1 = h.diam3d(0)
        lam = 0
        for i in range(int(h.n3d())):
            x2 = h.arc3d(i)
            d2 = h.diam3d(i)
            lam += (x2 - x1)/math.sqrt(d1 + d2)
            x1 = x2
            d1 = d2
        lam *=  math.sqrt(2) * 1e-5*math.sqrt(4*math.pi*freq*section.Ra*section.cm)

    return section.L / lam

def coarse(hoc_filename,cube_length,save_filename):
	# INPUT: NEURON .hoc filename to import (str), voxel cube side length (fl/str for gcd), name of file to create for mesh output (str)
	# 		>> cube_length: 'gcd' (gcd of box dims) OR floating-point (must be common factor of all box dims)
	# This function reads in NEURON data and passes their info to
	# coarse_gen(), then associates the tets of the STEPS Tetmesh object returned
	# to the NEURON sections which exist inside of them. 
	# Returns a tet_hoc dictionary -- tet_hoc[tet_index] = [encapsulated hoc section references] -- as well as the Tetmesh object

	## GET HOC SECTION INFO ##
	h.load_file(hoc_filename)

	allp = [[],[],[]]
	for s in h.allsec():
		for j in range(int(h.n3d())):
			allp[0].append(h.x3d(j))
			allp[1].append(h.y3d(j))
			allp[2].append(h.z3d(j))

	maxl = [max(allp[0]),max(allp[1]),max(allp[2])]
	minl = [min(allp[0]),min(allp[1]),min(allp[2])]
	bdim = [ maxl[0] - minl[0], maxl[1] - minl[1], maxl[2] - minl[2] ]
	print "dims: ", bdim
	print "mins: ", minl

	## CREATE COARSE MESH ##
	if (cube_length == 'gcd'):
		gcd = fractions.gcd(fractions.gcd(bdim[0],bdim[1]),fractions.gcd(bdim[2],bdim[1]))
		print "GCD: ", gcd
		cube_length = gcd

	sm = coarse_gen(cube_length,bdim,minl,save_filename)

	## ASSOCIATE HOC SECTIONS WITH THEIR TETS ##
	tet_hoc = tet_associate(sm[0])
	return tet_hoc, sm[0]

def morphology_to_dict(sections, outfile=None):
    section_map = {sec: i for i, sec in enumerate(sections)}
    result = []
    h.define_shape()

    for sec in sections:
        my_parent = parent(sec)
        my_parent_loc = -1 if my_parent is None else parent_loc(sec, my_parent)
        my_parent = -1 if my_parent is None else section_map[my_parent]
        n3d = int(h.n3d(sec=sec))
        result.append({
            'section_orientation': h.section_orientation(sec=sec),
            'parent': my_parent,
            'parent_loc': my_parent_loc,
            'x': [h.x3d(i, sec=sec) for i in xrange(n3d)],
            'y': [h.y3d(i, sec=sec) for i in xrange(n3d)],
            'z': [h.z3d(i, sec=sec) for i in xrange(n3d)],
            'diam': [h.diam3d(i, sec=sec) for i in xrange(n3d)],
            'name': sec.hname()           
        })

    if outfile is not None:
        with open(outfile, 'w') as f:
            json.dump(result, f)

    return result

    def retrieve_coordinate(self, sec):
        """Retrieve the coordinates of the section avoiding duplicates"""
        
        sec.push()
        x, y, z, d = [],[],[],[]

        tot_points = 0
        connect_next = False
        for i in range(int(h.n3d())):
            present = False
            x_i = h.x3d(i)
            y_i = h.y3d(i)
            z_i = h.z3d(i)
            d_i = h.diam3d(i)
            # Avoiding duplicates in the sec
            if x_i in x:
                ind = len(x) - 1 - x[::-1].index(x_i) # Getting the index of last value
                if y_i == y[ind]:
                    if z_i == z[ind]:
                        present = True
                    
            if not present:
                k =(x_i, y_i, z_i)
                x.append(x_i)
                y.append(y_i)
                z.append(z_i)
                d.append(d_i)                
        h.pop_section()
        #adding num 3d points per section
        self.n3dpoints_per_sec[sec.name()] = len(d)
        return (np.array(x),np.array(y),np.array(z),np.array(d))

def get_pos_data():
    """
    Get positions x, y, z for all segments and their diameter. 

    :returns: 
        4 lists: x,y,z,d. One element per section where each element is
        a :class:`~numpy.ndarray`.

    Example:
        .. code-block:: python

            x,y,z,d = get_pos_data()
            for sec in xrange(len(x)):
                for seg in xrange(len(x[sec]):
                    print x[sec][seg], y[sec][seg], z[sec][seg]
    """
    x = []
    y = []
    z = []
    d = []
    for sec in h.allsec():
        n3d = int(h.n3d())
        x_i, y_i, z_i = np.zeros(n3d), np.zeros(n3d), np.zeros(n3d),
        d_i = np.zeros(n3d)
        for i in xrange(n3d):
            x_i[i] = h.x3d(i)
            y_i[i] = h.y3d(i)
            z_i[i] = h.z3d(i)
            d_i[i] = h.diam3d(i)
        x.append(x_i)
        y.append(y_i)
        z.append(z_i)
        d.append(d_i)
    return x, y, z, d

def get_nseg( mysec ):
    """
    return the number of segments necesary to sample the smallest distance 
    between 2 consequtive 3D coordinates in a segment according to the
    Nysquid criterum (in our case 2 times smaller than the smallest 
    distance).
    """
    ncoord = int( h.n3d( sec = mysec) )

    # calculate distance between 0 and 1 coordinate
    x = h.x3d(0, sec = mysec) - h.x3d(1, sec = mysec)
    y = h.y3d(0, sec = mysec) - h.y3d(1, sec = mysec)
    z = h.z3d(0, sec = mysec) - h.z3d(1, sec = mysec)
    dist =  sqrt(x*x + y*y + z*z) 

    for i in range(1, ncoord-1):
        x = h.x3d(i, sec = mysec) - h.x3d(i+1, sec = mysec)
        y = h.y3d(i, sec = mysec) - h.y3d(i+1, sec = mysec)
        z = h.z3d(i, sec = mysec) - h.z3d(i+1, sec = mysec)
        
        value = sqrt(x*x + y*y + z*z)
        dist = value if value < dist else dist 

    nseg =  int((mysec.L/dist)*2.0) + 1 # odd number
    return( nseg ) 

    def draw_mayavi(self, x, y, z, d, edges):
        "Draw the surface the first time"

        # rendering disabled
        self.mayavi.visualization.scene.disable_render = True

        points = mlab.pipeline.scalar_scatter(x, y, z, d / 2.0)
        dataset = points.mlab_source.dataset
        dataset.point_data.get_array(0).name = "diameter"
        dataset.lines = np.vstack(edges)
        dataset.point_data.update()
        self.dataset = dataset

        # The tube
        src = mlab.pipeline.set_active_attribute(points, point_scalars="diameter")
        stripper = mlab.pipeline.stripper(src)
        tube = mlab.pipeline.tube(stripper, tube_sides=6, tube_radius=1)
        tube.filter.capping = True
        #        tube.filter.use_default_normal = False
        tube.filter.vary_radius = "vary_radius_by_absolute_scalar"
        self.tube = tube

        # Setting the voltage
        # Making room for the voltage
        v = []
        for sec in h.allsec():
            sec.push()
            v.extend(np.repeat(0.0, h.n3d()))
            h.pop_section()

        v = np.array(v)
        self.draw_surface(v, "v")

        # ReEnable the rendering
        self.mayavi.visualization.scene.disable_render = False

def get_section_path(h,sec):
    n3d = int(h.n3d(sec=sec))
    xyz = []
    for i in range(0,n3d):
        xyz.append([h.x3d(i,sec=sec),h.y3d(i,sec=sec),h.z3d(i,sec=sec)])
    xyz = np.array(xyz)
    return xyz

def retrieve_coordinate(sec):
    sec.push()
    x, y, z, d = [],[],[],[]
    area = 0
    tot_points = 0
    connect_next = False
    for i in range(int(h.n3d())):
        present = False
	x_i = h.x3d(i)
	y_i = h.y3d(i)
	z_i = h.z3d(i)
	d_i = h.diam3d(i)
	a_i = h.area(0.5)
	if x_i in x:
            ind = len(x) - 1 - x[::-1].index(x_i) # Getting the index of last value
	    if y_i == y[ind]:
                if z_i == z[ind]:
                    present = True
                    
	if not present:
            k =(x_i, y_i, z_i)
	    x.append(x_i)
	    y.append(y_i)
	    z.append(z_i)
	    d.append(d_i)                
	    area += np.sum(a_i)
    h.pop_section()
        #adding num 3d points per section
    n3dpoints[sec.name()] = [np.array(x),np.array(y),np.array(z),np.array(d)]
    return (np.array(x),np.array(y),np.array(z),np.array(d),area)

def accumulate_density(sec, density, domain):
  sec.push()
  for i in range(int(h.n3d())):
    x,y = (h.x3d(i), h.y3d(i))
    r = (round(x, domain[0]),round(y, domain[1]))
    if not False in r:
      density[r] += 1
  h.pop_section()

        def append_v(sec, v):
            """ Append data to v """
            sec.push()
            for ii in xrange(1, int(nrn.n3d())):
                v.append(sec.v)
            nrn.pop_section()

            return v