if (name_declared("pkgversions") != 4 ) { execute("strdef pkgversions") }
sprint(pkgversions,"%sReferenceAxis = $Revision: 1.5 $, ",pkgversions)
load_file("Segment.hoc")
/*
ReferenceAxis
ra = new ReferenceAxis()
ra.setup_with_secrefs(SecRef base, SecRef apex)
Sets up a new ReferenceAxis with base at the location of the 0th line
segment of the section referred to by the SecRef base and the and apex
located at the 0th line segment of the section referred to by the
SecRef apex.
ra.apply_expr_to_segref(SegmentRef seg, String expr)
Returns the value of expr computed for the segment referred to by seg
with respect to the axis defined by ra. This expression is in fact
computed for each geometric line segment and then these values are
summed to give a total for the electrical segment. Quantities that
may appear in expr are:
H - the projection of the centre of the line segment onto the axis
contained in ra (its "height").
l - the length of the line segment
d - the mean diameter of the line segment
If a geometrical line segment overlaps the boundary between electrical
segments, the quantities are appropriately interpolated.
Example:
ra.apply_expr_to_segref(seg, "h*l")
returns the sum of the products of the height and the length for each
line segment.
*/
begintemplate ReferenceAxis
// Public variables
public base, apex // SectionRefs of base and apex
objref base, apex
public axis, normaxis // Vector joining base and apex
objref axis
public verb // Verbose output
public adjustment
// Public functions
public setup_with_segrefs
public apply_expr_to_segref
// Private variables
objref pt3dbase, pt3dapex // Coordinates of base and apex
objref normaxis // Normalised axis
objref r, r0 // Final and initial points of line segments
objref mean_r // Mean position of line segment
objref Delta_r // Direction of line segment
objref this
strdef tmpstr
// Function definitions
proc init() {
rp = 0 // Relative probablity of connection in segment
verb = 0
adjustment = 0
}
// proc setup_with_secrefs () {
// $o1.secref.sec base = new SegmentRef(0)
// $o2.secref.sec apex = new SegmentRef(0)
// }
proc setup_with_segrefs () {
base = $o1
apex = $o2
if (numarg() == 3) {
adjustment = $3
}
pt3dbase = new Vector()
pt3dapex = new Vector()
axis = new Vector()
// Get the location of the base section. We will assume that the
// 0th element of the line segment list is a good enough indicator
// of the location
base.secref.sec {
pt3dbase.append(x3d(base.x))
pt3dbase.append(y3d(base.x))
pt3dbase.append(z3d(base.x))
}
// Get the location of the apex segment
apex.secref.sec {
pt3dapex.append(x3d(apex.x))
pt3dapex.append(y3d(apex.x))
pt3dapex.append(z3d(apex.x))
}
// Find unormalised axis
// axis = apex - base
axis = pt3dapex.c // Copying syntax
axis.sub(pt3dbase)
normaxis = axis.c
normaxis.div(sqrt(axis.sumsq()))
r = new Vector(3)
r0 = new Vector(3)
mean_r = new Vector(3)
Delta_r = new Vector(3)
}
func apply_expr_to_segref() { local s_beg s_end l s0 p i d d0 s d1seg d0seg
// args: SegRef, probexpr
rp = 0 // rp can't be local since it must be accesed by execute
p = 0 // Probability of synapse in segment
// Find out the start and end points along the segment
if ( ( $o1.x != 0 ) && ($o1.x != 1) ) {
$o1.secref.sec s_beg = ($o1.x - 0.5/nseg) * L
$o1.secref.sec s_end = s_beg + L/nseg
if (verb==1) print "s_beg, s_end: ", s_beg, s_end
// Now look through each geometric segment in turn
s = 0
s0 = 0 // Old value of s
$o1.secref.sec for i=1,(n3d()-1) {
if (verb==1) print "coords: ", x3d(i), y3d(i), z3d(i)
r.x(0) = x3d(i)
r.x(1) = y3d(i)
r.x(2) = z3d(i)
r0.x(0) = x3d(i-1)
r0.x(1) = y3d(i-1)
r0.x(2) = z3d(i-1)
d1 = diam3d(i)
d0 = diam3d(i-1)
Delta_r = r.c
Delta_r.sub(r0)
s0 = s
ds = sqrt(Delta_r.sumsq())
s = s + ds
if (verb==1) print "s, area: ", s, s*PI*diam3d(i)
if (s <= s_beg) {
continue
}
if (s0 >= s_end) {
break
}
l = ds // Length of line segment
d0seg = d0
d1seg = d1
// Begining of segment
if ((s>s_beg) && (s0<s_beg)) {
r0.add(Delta_r.mul((s_beg-s0)/ds))
l = l - (s_beg - s0)
d0seg = d0 + (d1-d0)/ds*(s_beg - s0)
}
// End of segment
if ((s>s_end) && (s0<s_end)) {
r.sub(Delta_r.mul((s-s_end)/ds))
l = l - (s - s_end)
d1seg = d1 - (d1-d0)/ds*(s - s_end)
}
mean_r = r0.c
mean_r.add(r)
mean_r.div(2)
mean_r.sub(pt3dbase)
H = mean_r.dot(normaxis) + adjustment
d = (d0seg + d1seg)/2
if (verb==1) print "h, d, l, area: ", H, d, l, PI*d*l
sprint(tmpstr,"rp += %s",$s2)
execute(tmpstr,this)
if (verb==1) print rp
}
} else { // x == 0 or x == 1
l = 0
if ( $o1.x == 0 ) { i = 0 }
if ( $o1.x == 1 ) { i = n3d()-1 }
d = diam3d(i)
r.x(0) = x3d(i)
r.x(1) = y3d(i)
r.x(2) = z3d(i)
r.sub(pt3dbase)
H = mean_r.dot(normaxis) + adjustment
if (verb==1) print "h, d, l, area: ", H, d, l, PI*d*l
sprint(tmpstr,"rp = %s",$s2)
execute(tmpstr,this)
}
return rp
}
endtemplate ReferenceAxis
