################ FN model ################

# declare ode system for FN model
# dv/dt = ...
# dr/dt = ...
# input: 
#	t: current time point
#	var: c(v,r) current state
#	p: params c(a,b,c,I)
# output: list(c(dv/dt,dr/dt))
FN.ode.system <- function(t, var, p) {
	v <- var[1]; r <- var[2]
	a <- p[1]; b <- p[2]; c <- p[3]; I <- p[4]

	return(list(c(
		# write the function dv/dt after this line
		c * (v - 1/3*v^3 + r + I),	
		# write the function dr/dt after this line
		-1/c * (v - a + b*r)
	)))
}

# given v, find r such that 
# dv/dt = r - f(v) = 0
# dr/dt = r - g(v) = 0
FN.nullcline <- function(x, p) {
	v <- x
	a <- p[1]; b <- p[2]; c <- p[3]; I <- p[4]

	return(c(
		# write the function f(v) after this line
		-v + 1/3*v^3 - I,
		# write the function g(v) after this line
		(-v+a)/b
	))
}

########### run experiment ##########

	test.ode.system	= FN.ode.system 
	test.nullcline	= FN.nullcline
	test.params		= c(0.7,0.8,3,0)
	test.xrange		= c(-3,3)
	test.yrange		= c(-3,3) 
	test.xstep		= 0.25
	test.ystep		= 0.25
	test.xlabel		= "v"
	test.ylabel		= "r"
	test.x.ini		= c(3)
	test.y.ini		= c(-3)
	test.times		= seq(1,100,by=0.01)
	test.ode.root	= NULL
	test.ode.event	= NULL

draw.phase.plane(
			ode.sys = test.ode.system, ode.event = test.ode.event, ode.root = test.ode.root,
			nullcline = test.nullcline, 
			params = test.params,
			xrange = test.xrange, yrange = test.yrange, 
			xstep = test.xstep, ystep = test.ystep,
			xlabel = test.xlabel, ylabel = test.ylabel,
			x.ini = test.x.ini, y.ini = test.y.ini, times = test.times)

#draw.trajectory.over.time(
#						ode.sys = test.ode.system, ode.event = test.ode.event, ode.root = test.ode.root, 
# 						params = test.params,
#						x.ini = test.x.ini, y.ini = test.y.ini, 
#						times = test.times,
#						xlabel = test.xlabel, ylabel = test.ylabel)