# disclaimer and manual
print
print ("Application created by Piotr Stuglik.")
print
print ("Commands: q - quit, manu - manual, run - will run commands essential to plot the functions, commands - list of commands, test - sets initial parameters, scope - asks for the number of reagents, summary - displays the data table, conc - asks for the initial concentrations, kin - asks for the kinetic constants, graphit - plots.")
print
print ("This application allows one to plot the concentration of reagents as a function of time for consecutive reactions A->B->C->...")
print
print ("n - number of reagents, k_n - kinetic constants, c_0_n - initial concentrations, c_n(t) - concentration, t - time.")
print
print ("NOTE: kinetic constant k_n values cannot be equal to one another!")
print
print ("Input 'run' to start.")
print

# define initial values
n = 0
k_n = 0
c_0_n = 0

# import necessary modules
import matplotlib.pyplot as plt

import numpy as np

# define manual command
def manu():
	print
	print ("Commands: q - quit, run - will run commands essential to plot the functions, commands - list of commands, test - sets initial parameters, scope - asks for the number of reagents, summary - displays the data table, conc - asks for the initial concentrations, kin - asks for the kinetic constants, graphit - plots.")
	print
	print ("This application allows one to plot the concentration of reagents as a function of time for consecutive reactions A->B->C->...")
	print
	print ("n - number of reagents, k_n - kinetic constants, c_0_n - initial concentrations, c_n(t) - concentration, t - time.")
	print
	print ("NOTE: kinetic constant k_n values cannot be equal to one another!")
	print
	print ("Input 'run' to start.")
	print
	return

# define run command
def run():
	scope()
	conc()
	kin()
	summary()
	graphit()
	return

# define default testing parameters
def test():
	global n, scope_n, k_n, c_0_n
	n = 5
	scope_n = range(1, n + 1)
	k_n = [1,1.5,2,2.5,0]
	c_0_n = [10,0,0,0,0]
	return

# displays the list of commands
def commands():
	print
	print ("Commands: q - quit, manu - manual, run - will run commands essential to plot the functions, commands - list of commands, test - sets initial parameters, scope - asks for the number of reagents, summary - displays the data table, conc - asks for the initial concentrations, kin - asks for the kinetic constants, graphit - plots.")
	print
	return

# number of reagents query
def scope():
	global n, scope_n, c_0_n, k_n
	while True:
		try:
			n = raw_input("Define the number of n reagents: ")
			if (int(n) < 1):
				print
				print "INVALID INPUT! Provide integer greater than 0."
				print
				continue
			elif (int(n) > 0):
				n = int(n)
				scope_n = range(1, n + 1)
				c_0_n = [int(0)] * n
				k_n = [int(0)] * n
				break
		except (ValueError):
			print
			print "INVALID INPUT! Provide integer greater than 0."
			print
	return

# initial concentrations query
def conc():
	if (n == 0):
		scope()
	while True:
		try:
			y = raw_input("Define the value of c_0_n for n equal to (press 0 to break): ")
			if int(y) == 0:
				break
			elif (int(y) > n or int(y) < 1):
				print
				print "INVALID INPUT! Argument out of range."
				print
				continue
			y = int(y)
			x = raw_input("Define the value of c_0_" + str(y) + ": ")
			if "." in x:
				c_0_n[y - 1] = float(x)
			else:
				c_0_n[y - 1] = int(x)
		except (ValueError):
			print
			print "INVALID INPUT! Provide integer for n and integer/float for c_0_n."
			print
	return

# kinetic constants query
def kin():
	if (n == 0):
		scope()
	while True:
		try:
			q = raw_input("Define the value of k_n for n equal to (press 0 to break): ")
			if int(q) == 0:
				break
			elif (int(q) > n or int(q) < 1):
				print
				print "INVALID INPUT! Argument out of range."
				print
				continue
			q = int(q)
			p = raw_input("Define the value of k_" + str(q) + ": ")
			if "." in p:
				k_n[q - 1] = float(p)
			else:
				k_n[q - 1] = int(p)
		except (ValueError):
			print
			print "INVALID INPUT! Provide integer for n and integer/float for k_n."
			print
	return

# display the table with the initial data
def summary():
	if n == 0:
		scope()
		print
		print "n:     ", scope_n
		print "c_0_n: ", c_0_n
		print "k_n:   ", k_n
		print
	else:
		print
		print "n:     ", scope_n
		print "c_0_n: ", c_0_n
		print "k_n:   ", k_n
		print
	return

# define a secondary function
def m(x,y,w,z):
    if(z > x):
        return "np.exp(-k_n[%i - 1] * t)" % w
    else:
        return ("float(k_n[%i - 1])/(float(k_n[%i - 1]) - float(k_n[%i - 1])) * " % (y,z,w)) + "(" + m(x,y+1,w,z+1) + " - " + m(x,y+1,z,z+1) + ")"

# x is the concentration index
# y is the yth term of c_x(t)
def q(x,y):
	return "c_0_n[%i - 1] * " % y + m(x,y,y,y+1)

# define the final function
def c_(u):
	global c_n
	c_n = ""
	for el in range(1,u+1):
		c_n = c_n + q(u,el) + " + "
	if c_n.endswith(" + "):
		c_n = c_n[:-3]
	return c_n

# graphit subdefinition 1
def reag():
	global reag_var
	while True:
		try:
			reag_var = map(int, raw_input("Provide the reagents to plot (separate with spacebar): ").split(" "))
			if all(i in scope_n for i in reag_var):
				return
			else:
				print
				print "INVALID INPUT! Provide integers between 1 and n."
				print
				continue
		except (ValueError):
			print
			print "INVALID INPUT! Provide integers between 1 and n."
			print
	return

# graphit subdefinition 2
def t_k():
	global t_k_var
	while True:
		try:
			t_k_var = float(raw_input("Define the time range from 0 to: "))
			if (float(t_k_var) <= 0):
				print
				print "INVALID INPUT! Provide integer/float greater than 0."
				print
				continue
			elif (float(t_k_var) > 0):
				return
		except (ValueError):
			print
			print "INVALID INPUT! Provide integer/float greater than 0."
			print
	return

# graphit subdefinition 3
def t_d():
	global t_d_var
	while True:
		try:
			t_d_var = float(raw_input("Define the precision of the time axis (e.g. 0.01): "))
			break
		except (ValueError):
			print
			print "INVALID INPUT! Provide integer/float."
			print
	return

# plot the requested functions
def graphit():
	if (n == 0):
		scope()
		conc()
		kin()
		summary()
	elif (all(ele1 == 0 for ele1 in k_n) or all(ele2 == 0 for ele2 in c_0_n)):
		conc()
		kin()
		summary()
	try:
		global t
		reag()
		t_k()
		t_d()
		t = np.arange(0,t_k_var,t_d_var)
		p = []

		plt.xlabel("time t")
		plt.ylabel("concentration c_n(t)")
		for i in reag_var: # the actual plot command
			p += plt.plot(t,eval(c_(i)),label= "c_" + str(i) + "(t)")
			max_c = max(eval(c_(i)))
			max_t = t[eval(c_(i)).argmax()]
			print "t_%i_max =" % (i), max_t
			print "c_%i_max =" % (i), max_c
		plt.legend(loc="center right")
		plt.show()
	except (ZeroDivisionError):
		print
		print "INVALID INPUT! Kinetic constant k_n values cannot be equal to one another."
		print
	return

# interface below; stand-by mode
while True:
	command = raw_input("Provide command: ")
	allowed = ["scope", "commands", "graphit", "kin", "conc", "summary", "test", "run", "manu"]
	if command != "q":
		comm = "%s()" % (command)
		if command not in allowed:
			print
			print "INVALID COMMAND! Provide command from the list."
			print
		else:
			exec(comm)
	elif command == "q":
		break