neuron 7.5-1 / usr / share / nrn / lib / hoc / subiter.hoc

begintemplate SubsetDomainIterator
public sl, x, p, p0, p1, loop, update
objref sl, map, this, tobj

proc init() {
	nsc = 0
	sl = $o1
	x = 0
	p = 0
	p0 = p1 = 0
	map = new Vector()

	// distance_from_root, distance_from_origin
	// radial_from_origin, project_onto_line
	// position_from_line
	// only implemented for 0
	path_style_ = $2
	if (path_style_ == 1 || path_style_ > 3) {
		printf("path_style_ %d not implemented\n", path_style_)
		execerror("this path_style not implemented", "")
	}
	
	// no origin translation, 0 at most proximal end,
	// 0 at proximal disjoint ends
	// only implemented for 0 and 1
	proximal_style_ = $3
	if (proximal_style_ == 2) {
		printf("proximal_style_ %d not implemented\n", proximal_style_)
		execerror("this path_style not implemented", "")
	}

	// no distal end normalization, 1 at most distal end,
	// 1 at distal ends (longer path has uniform metric across branch point)
	// only implemented for 0 and 1
	distal_style_ = $4
	if (distal_style_ == 2) {
		printf("distal_style_ %d not implemented\n", distal_style_)
		execerror("this path_style not implemented", "")
	}
	if (path_style_ == 2) {
		xo_ = $5
		yo_ = $6
		zo_ = $7
	}
	if (path_style_ == 3) {
		origin_ = $5
		theta_ = $6
	}
}

iterator loop() {local i, usearg
	update()
	if (numarg() == 2) { usearg = 1 } else { usearg = 0 }
	i = 0	
	forsec sl for (x, 0) {
		p = map.x[i]  i += 1
		if (usearg) {
			$&1 = x
			$&2 = p
		}
		iterator_statement
	}
}

// everything following is map specific
proc update() {local i
	if (nsc == nrn_shape_changed_) {
		return
	}
//	printf("updating map for %s nsc=%d  nrn_shape_changed_=%d\n", this, nsc, nrn_shape_changed_)
	nsc = nrn_shape_changed_
	i = 0
	forsec sl for (x, 0) { i += 1 }
	map.resize(i)
	if (path_style_ == 0) { // path distance from root
		ps0()
	}else if (path_style_ == 2) {
		ps2()
	}else if (path_style_ == 3) {
		ps3()
	}
	if (proximal_style_ == 1) {
		map.sub(p0)
		p1 -= p0
		p0 = 0
	}
	if (distal_style_ == 1) {
		map.div(p1)
		p0 /= p1
		p1 = 1
	}
}

// path distance from root
// assume sectionlist is a subset of a single cell
proc ps0() {local i, x, d
	i = 0
	p0 = 1e9
	p1 = 1
	forsec sl for (x, 1) {
		if (i == 0) {
			tobj = new SectionRef()
			tobj.root { distance() }
			objref tobj
		}
		d = distance(x)
		if (x > 0 && x < 1) {
			map.x[i] = d
			i += 1
		}
		if (d < p0) { p0 = d }
		if (d > p1) { p1 = d }
	}
}

// radial distance from origin
proc ps2() {local i, a, x, y, z, d
	i = 0
	p0 = 1e9
	p1 = 1
	forsec sl {
		i3 = 1
		for (a, 1) {
			seg_from_3d(a)
			d = sqrt((x3 - xo_)^2 + (y3 - yo_)^2 + (z3 - zo_)^2)
			if (a > 0 && a < 1) {
				map.x[i] = d
				i += 1
			}
			if (d < p0) { p0 = d }
			if (d > p1) { p1 = d }
		}
	}
}

// projection onto line at angle theta with 0 defined at origin
proc ps3() {local i, a, c, s, d
	c = cos(theta_/DEG)  s = sin(theta_/DEG)
	i = 0
	p0 = 1e9
	p1 = 1
	forsec sl {
		i3 = 1
		for (a, 1) {
			seg_from_3d(a)
			d = x3*c + y3*s
			if (a > 0 && a < 1) {
				map.x[i] = d
				i += 1
			}
			if (d < p0) { p0 = d }
			if (d > p1) { p1 = d }
		}
	}
	map.sub(origin_)
	p0 -= origin_
	p1 -= origin_
}

proc seg_from_3d() {local th, a, i
    if ($1 == 0) {
	i = 0
	a = arc3d(i)
	x3 = x3d(i)
	y3 = y3d(i)
	z3 = z3d(i)
    }else if ($1 == 1) {
	i = n3d()-1
	a = arc3d(i)
	x3 = x3d(i)
	y3 = y3d(i)
	z3 = z3d(i)
    }else{
	$1 *= L
	while($1 > arc3d(i3)) {
		i3 += 1
	}
	i = i3 - 1
	a = arc3d(i)
	th = ($1 - a)/(arc3d(i3) - a)
	x3 = (1 - th)*x3d(i) + th*x3d(i3)
	y3 = (1 - th)*y3d(i) + th*y3d(i3)
	z3 = (1 - th)*z3d(i) + th*z3d(i3)
    }
}
endtemplate SubsetDomainIterator