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