// Modified from interpxyz.hoc and setpointers.hoc of Ted Carnevale
// $Id: interpxyz.hoc,v 1.2 2005/09/10 23:02:15 ted Exp $
// 
// Computes xyz coords of nodes in a model cell whose topology & geometry are 
// defined by pt3d data. Expects sections to already exist, and that the xtra 
// mechanism has been inserted


objref xx, yy, zz, length       // original, irregularly spaced
objref xint, yint, zint, range  // interpolated, spaced at regular intervals


proc grindaway() { local ii, nn, kk, xr
    forall {
        if (ismembrane("xtra")) {
            // get the data for the section
            nn = n3d()
            xx = new Vector(nn)
            yy = new Vector(nn)
            zz = new Vector(nn)
            length = new Vector(nn)

            for ii = 0,nn-1 {
                xx.x[ii] = x3d(ii)
                yy.x[ii] = y3d(ii)
                zz.x[ii] = z3d(ii)
                length.x[ii] = arc3d(ii)
            }

            // to use Vector class's .interpolate() 
            // must first scale the independent variable
            // i.e. normalize length along centroid
            length.div(length.x[nn-1])

            // initialize the destination "independent" vector
            range = new Vector(nseg+2)
            range.indgen(1/nseg)
            range.sub(1/(2*nseg))
            range.x[0]=0
            range.x[nseg+1]=1

            // length contains the normalized distances of the pt3d points 
            // along the centroid of the section.  These are spaced at 
            // irregular intervals.
            // range contains the normalized distances of the nodes along the 
            // centroid of the section.  These are spaced at regular intervals.
            // Ready to interpolate.

            xint = new Vector(nseg+2)
            yint = new Vector(nseg+2)
            zint = new Vector(nseg+2)
            xint.interpolate(range, length, xx)
            yint.interpolate(range, length, yy)
            zint.interpolate(range, length, zz)

            // for each node, assign the xyz values to x_xtra, y_xtra, z_xtra
            // for ii = 0, nseg+1 {
            // don't bother computing coords of the 0 and 1 end also avoid writing 
            // coords of the 1 end into the last internal node's coords
            for ii = 1, nseg {
                xr = range.x[ii]
                x_xtra(xr) = xint.x[ii]
                y_xtra(xr) = yint.x[ii]
                z_xtra(xr) = zint.x[ii]
            }
        }
    }
}

proc setpointers() { local done
    // determines interpolated locations of nodes
    grindaway()
    forall {
        if (ismembrane("xtra")) {
	        for (x, 0) {
		        setpointer im_xtra(x), i_membrane(x)
		        setpointer ex_xtra(x), e_extracellular(x)
	        }
        }
    }
}


/////////////////////////////////////////////////////////////////////

print "INFO: Changes to geometry or nseg should be followed by setpointers()"

setpointers()