/* btree6.hoc
   creates several binary dendritic trees
     and attaches them to a core of branches
   the trees are stylized, the core is pt3d
     (a trick to "zigzag" shape plot so trees don't overlap)

   revised 2/5/2010
*/

load_file("nrngui.hoc")

NUMTREES = 4	// how many to build
CORE_L = 300	// dimensions of the core sections
CORE_D = 16

objref norders, mxL, mnL, mxdiam, kexp
norders = new Vector(NUMTREES)	// how many orders of branching in a subtree
mxL = new Vector(NUMTREES)	// max length of a section in a subtree
mnL = new Vector(NUMTREES)	// min length of a section
mxdiam = new Vector(NUMTREES)	// max diam of a section
kexp = new Vector(NUMTREES)	// dparent^k = 2*ddaughter^k

ii = 0
norders.x[ii] = 3
mxL.x[ii] = 250
mnL.x[ii] = 50
mxdiam.x[ii] = 10
kexp.x[ii] = 1.5	// like a cylinder with const diam

ii = 1
norders.x[ii] = 2
mxL.x[ii] = 250
mnL.x[ii] = 50
mxdiam.x[ii] = 10
kexp.x[ii] = 2.0	// like a "flaring cylinder"

ii = 2
norders.x[ii] = 2
mxL.x[ii] = 250
mnL.x[ii] = 50
mxdiam.x[ii] = 10
kexp.x[ii] = 1.0	// like a "tapering cylinder"

ii = 3
norders.x[ii] = 3
mxL.x[ii] = 250
mnL.x[ii] = 50
mxdiam.x[ii] = 10
kexp.x[ii] = 1.0	// a bigger "tapering cyl" tree

/* arguments:
   $1	how many orders of branching
   $2	length of 0th order branch
   $3	length of highest order branches
   $4	diam of 0th order branch
   $5	exponent for branching rule

   all branches are public
*/

begintemplate Subtree
	public dend

	create dend[1]

	/* functions that specify dimensions */

	// arg == order of branch
	func Lfunc() { local normorder, Lspan
		normorder = $1/(numorders - 1)
		Lspan = maxL - minL

		// linear
//		return (maxL - normorder * Lspan)

		// parabolic
		return (minL + Lspan * (normorder - 1)^2)
	}

	// arg == order of branch
	func dfunc() {
		return (maxdiam / (2^($1/k)))
	}


	proc init() { local ii
		numorders = $1
		numbranches = -1 + 2^numorders
		maxL = $2
		minL = $3
		maxdiam = $4
		k = $5

		create dend[numbranches]

		// index of the first of the highest-order branches
		temp = -1 + 2^(numorders-1)
		ii = 0
		jj = 0
		while (ii<temp) {
			dend[ii] {
				jj+=1
				connect dend[jj](0), 1
				jj+=1
				connect dend[jj](0), 1
			}
			ii+=1
		}

		// set dimensions
		for (order=0; order<numorders; order+=1) {
			temp = 2^order
			istart = temp - 1
			ifin = 2*istart
			for (ii=istart; ii<=ifin; ii+=1) {
				dend[ii] {
					L = Lfunc(order, numorders)
					diam = dfunc(order)
				}
			}
		}
	}
endtemplate Subtree



create core[3]
access core[1]

pt3dclear()

xtemp=0
ytemp=0
sx = 1

for ii = 0,2 {
	core[ii] {
		pt3dadd(xtemp,ytemp,0,CORE_D)
		ytemp += CORE_L / sqrt(2)
		xtemp += sx*CORE_L / sqrt(2)
		pt3dadd(xtemp,ytemp,0,CORE_D)
		sx = -sx
	}
}


objref stree[NUMTREES]

for ii = 0,NUMTREES-1 {
	stree[ii] = \
	new Subtree(norders.x[ii], mxL.x[ii], mnL.x[ii], \
			 mxdiam.x[ii], kexp.x[ii])
}

connect core[1](0), core[0](1)
connect core[2](0), core[1](1)

connect stree[0].dend[0](0), core[0](0)
connect stree[1].dend[0](0), core[0](1)
connect stree[2].dend[0](0), core[1](1)
connect stree[3].dend[0](0), core[2](1)


topology()  // reveal branched architecture in TTY format

objref sh
sh = new PlotShape(0)
sh.size(-481.198,692.596,-269.066,904.728)
sh.variable("v")
sh.view(-481.198, -269.066, 1173.79, 1173.79, 100, 105, 200.64, 200.32)


// just for the heck of it

forall {
	insert pas
	e_pas = -65
}