begintemplate CrossCorrelation
public reference, test
public ydat, tdat, grid, cor, y, tcor
external hoc_obj_
objref ydat[2], tdat[2] // original full data
objref y[2] // subset of data interpolated onto power of 2 grid
objref grid
objref cor, tcor
objref this, b, gdat, gcor
strdef tstr


proc init() {
	if (numarg() < 2) {
		ydat[0] = new Vector(2)
		tdat[0] = new Vector(2)
		ydat[0].x[0] = 2
		tdat[0].indgen(100)
		reference(ydat[0], tdat[0])
		test(ydat[0], tdat[0])
	}else if (numarg() == 2) {
		reference($o1, $o2)
		test($o1, $o2)
	}else if (numarg() > 2) {
		reference($o1, $o2)
		test($o3, $o4)
	}
	nominal_dt = .025 //ms	
	usemean_ = 1
	cor = new Vector()
	tcor = new Vector()
	grid = new Vector()
	y[0] = new Vector()
	y[1] = new Vector()
	build()
	map()
	bcon(0)
	gdat.exec_menu("View = plot")
	gcor.exec_menu("View = plot")
	gdat.exec_menu("Select Boundaries")
}

proc map() {
	sprint(tstr, "%s", this)
	b.map(tstr)
}

proc reference() {
	ydat[0] = $o1.cl
	tdat[0] = $o2.c
	tbegin = 150
	tend = 300
	maxdat = ydat[0].max
	mindat = ydat[0].min
	maxshift = 10 //ms
}

proc test() {
	ydat[1] = $o1.cl
	tdat[1] = $o2.c
}

proc correl() {local ibegin, iend, i, x
	x = (tend  - tbegin)/nominal_dt
//print "tbegin=",tbegin, " tend=", tend, " x=", x
	x = int(log(x)/log(2))
	grid.resize(2^x)
	ddt = (tend  - tbegin)/grid.size
	grid.indgen(tbegin, ddt)
//print "grid size =", grid.size, " ddt=", ddt
	for i=0, 1 {
		ydat[i].line(gdat, tdat[i], 1, 1)
		ibegin=tdat[i].indwhere(">", tbegin)
		if (ibegin == -1) { ibegin = tdat[i].size-2 }
		iend = tdat[i].indwhere(">", tend)
		if (iend == -1) { iend = tdat[i].size-1 }
//print "ibegin=", ibegin, " iend=", iend
		y[i].copy(ydat[i], ibegin, iend)
		y[i].interpolate(grid, tdat[i].c(ibegin, iend))
		y[i].line(gdat, grid, 2, 1)
		if (usemean_) {
			y[i].sub(y[i].mean)
		}
		if (zcon_ == 1) {
			y[i].resize(y[i].size*2)
		}
	}
	cor.correl(y[0], y[1]).mul(ddt)
	actual_shift = maxshift
	if (maxshift*2 >= (tend - tbegin)) {
		actual_shift = (tend - tbegin)/2
	}
	cor.rotate(actual_shift/ddt)
	cor.resize(2*actual_shift/ddt+1)
	tcor.resize(2*actual_shift/ddt+1).indgen(-actual_shift, ddt)
}

proc build() {
	b = new VBox()
	b.ref(this)
	b.save("")
	b.intercept(1)
	gdat = new Graph()
	xpanel("")
	gdat.menu_tool("Select Boundaries", "bound")
	xcheckbox("Periodic Boundary Condition", &pcon_, "bcon(0)")
	xcheckbox("Zero  Boundary Condition", &zcon_, "bcon(1)")
	xcheckbox("Use mean = 0", &usemean_, "change()")
	xbutton("Reference from Clipboard", "bref()")
	xbutton("Test from Clipboard", "btes()")
	xpvalue("maxshift", &maxshift, 1, "change_maxshift()")
	xpvalue("nominal dt", &nominal_dt, 1, "change()")
	xpvalue("actual dt", &ddt)
	xpvalue("tbegin", &tbegin, 1, "region_change()")
	xpvalue("tend", &tend, 1, "region_change()")
	xpanel()
	gcor = new Graph()
	b.intercept(0)
}

proc bref() {
	reference(hoc_obj_, hoc_obj_[1])
	pl()
	gdat.exec_menu("View = plot")
	gcor.exec_menu("View = plot")
}

proc btes() {
	test(hoc_obj_, hoc_obj_[1])
	pl()
	gdat.exec_menu("View = plot")
	gcor.exec_menu("View = plot")
}

proc bcon() {
	pcon_ = ($1 == 0)
	zcon_ = ($1 == 1)
	pl()
}

proc change() { pl() }

proc change_maxshift() {
	pl()
	gcor.size(-maxshift, maxshift, gcor.size(3), gcor.size(4))
}

proc region_change() {
	if (tend > tdat[0].x[tdat[0].size-1]) {
		tend = tdat[0].x[tdat[0].size-1]
		if (tbegin + 4*nominal_dt > tend) {
			tbegin = tend - 4*nominal_dt
		}
	}
	if (tbegin < tdat[0].x[0]) {
		tbegin = tdat[0].x[0]
	}
	if (tbegin + 4*nominal_dt > tend) {
		tend = tbegin+4*nominal_dt
	}
	pl()
}

proc pl() {
	gdat.erase()
	gdat.vfixed(1)
	gdat.label(.5, .9)
	correl()
	bndry(tbegin)
	bndry(tend)
	gcor.erase()
	cor.line(gcor, tcor, 1,1)
}

proc bndry() {
	gdat.beginline(3,1)
	gdat.line($1, maxdat)
	gdat.line($1, mindat)
	gdat.flush()
}

proc bound() {
//	print $1, $2, $3
	if ($1 == 2) {
		if (abs($2 - tbegin) < abs($2 - tend)) {
			beg = 1
		}else{
			beg = 0
		}
	}
	if (beg) {
		tbegin = $2
	}else{
		tend = $2
	}
	region_change()
}

endtemplate CrossCorrelation

objref correl, ydat, tdat

/*
clipboard_retrieve("tmp1.dat")
ydat = hoc_obj_.cl
tdat = hoc_obj_[1].c
clipboard_retrieve("tmp2.dat")
correl = new CrossCorrelation(ydat, tdat, hoc_obj_, hoc_obj_[1])

ydat.resize(1000)
ydat.fill(0)
for i=0, 99 ydat.x[i] = 2
tdat.resize(1000)
tdat.indgen(.01)
correl = new CrossCorrelation(ydat, tdat, ydat, tdat)
*/