1. # disclaimer and manual
  2. print
  3. print ("Application created by Piotr Stuglik.")
  4. print
  5. 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.")
  6. print
  7. print ("This application allows one to plot the concentration of reagents as a function of time for consecutive reactions A->B->C->...")
  8. print
  9. print ("n - number of reagents, k_n - kinetic constants, c_0_n - initial concentrations, c_n(t) - concentration, t - time.")
  10. print
  11. print ("NOTE: kinetic constant k_n values cannot be equal to one another!")
  12. print
  13. print ("Input 'run' to start.")
  14. print
  15.  
  16. # define initial values
  17. n = 0
  18. k_n = 0
  19. c_0_n = 0
  20.  
  21. # import necessary modules
  22. import matplotlib.pyplot as plt
  23.  
  24. import numpy as np
  25.  
  26. import decimal as deci
  27. deci.getcontext().prec = 28
  28.  
  29. # define manual command
  30. def manu():
  31.     print
  32.     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.")
  33.     print
  34.     print ("This application allows one to plot the concentration of reagents as a function of time for consecutive reactions A->B->C->...")
  35.     print
  36.     print ("n - number of reagents, k_n - kinetic constants, c_0_n - initial concentrations, c_n(t) - concentration, t - time.")
  37.     print
  38.     print ("NOTE: kinetic constant k_n values cannot be equal to one another!")
  39.     print
  40.     print ("Input 'run' to start.")
  41.     print
  42.     return
  43.  
  44. # define run command
  45. def run():
  46.     scope()
  47.     conc()
  48.     kin()
  49.     summary()
  50.     graphit()
  51.     return
  52.  
  53. # define default testing parameters
  54. def test():
  55.     global n, scope_n, k_n, c_0_n
  56.     n = 5
  57.     scope_n = range(1, n + 1)
  58.     k_n = [1,1.5,2,2.5,0]
  59.     c_0_n = [10,0,0,0,0]
  60.     return
  61.  
  62. # displays the list of commands
  63. def commands():
  64.     print
  65.     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.")
  66.     print
  67.     return
  68.  
  69. # number of reagents query
  70. def scope():
  71.     global n, scope_n, c_0_n, k_n
  72.     while True:
  73.         try:
  74.             n = raw_input("Define the number of n reagents: ")
  75.             if (int(n) < 1):
  76.                 print
  77.                 print "INVALID INPUT! Provide integer greater than 0."
  78.                 print
  79.                 continue
  80.             elif (int(n) > 0):
  81.                 n = int(n)
  82.                 scope_n = range(1, n + 1)
  83.                 c_0_n = [int(0)] * n
  84.                 k_n = [int(0)] * n
  85.                 break
  86.         except (ValueError):
  87.             print
  88.             print "INVALID INPUT! Provide integer greater than 0."
  89.             print
  90.     return
  91.  
  92. # initial concentrations query
  93. def conc():
  94.     if (n == 0):
  95.         scope()
  96.     while True:
  97.         try:
  98.             y = raw_input("Define the value of c_0_n for n equal to (press 0 to break): ")
  99.             if int(y) == 0:
  100.                 break
  101.             elif (int(y) > n or int(y) < 1):
  102.                 print
  103.                 print "INVALID INPUT! Argument out of range."
  104.                 print
  105.                 continue
  106.             y = int(y)
  107.             x = raw_input("Define the value of c_0_" + str(y) + ": ")
  108.             if "." in x:
  109.                 c_0_n[y - 1] = deci.Decimal(x)
  110.             else:
  111.                 c_0_n[y - 1] = int(x)
  112.         except (ValueError):
  113.             print
  114.             print "INVALID INPUT! Provide integer for n and integer/float for c_0_n."
  115.             print
  116.     return
  117.  
  118. # kinetic constants query
  119. def kin():
  120.     if (n == 0):
  121.         scope()
  122.     while True:
  123.         try:
  124.             q = raw_input("Define the value of k_n for n equal to (press 0 to break): ")
  125.             if int(q) == 0:
  126.                 break
  127.             elif (int(q) > n or int(q) < 1):
  128.                 print
  129.                 print "INVALID INPUT! Argument out of range."
  130.                 print
  131.                 continue
  132.             q = int(q)
  133.             p = raw_input("Define the value of k_" + str(q) + ": ")
  134.             if "." in p:
  135.                 k_n[q - 1] = deci.Decimal(p)
  136.             else:
  137.                 k_n[q - 1] = int(p)
  138.         except (ValueError):
  139.             print
  140.             print "INVALID INPUT! Provide integer for n and integer/float for k_n."
  141.             print
  142.     return
  143.  
  144. # display the table with the initial data
  145. def summary():
  146.     if n == 0:
  147.         scope()
  148.         print
  149.         print "n:     ", scope_n
  150.         print "c_0_n: ", c_0_n
  151.         print "k_n:   ", k_n
  152.         print
  153.     else:
  154.         print
  155.         print "n:     ", scope_n
  156.         print "c_0_n: ", c_0_n
  157.         print "k_n:   ", k_n
  158.         print
  159.     return
  160.  
  161. # define a secondary function
  162. def m(x,y,w,z):
  163.     if(z > x):
  164.         return "deci.Decimal(-k_n[%i - 1] * t).exp()" % w
  165.     else:
  166.         return ("deci.Decimal(k_n[%i - 1])/(deci.Decimal(k_n[%i - 1]) - deci.Decimal(k_n[%i - 1])) * " % (y,z,w)) + "(" + m(x,y+1,w,z+1) + " - " + m(x,y+1,z,z+1) + ")"
  167.  
  168. # x is the concentration index
  169. # y is the yth term of c_x(t)
  170. def q(x,y):
  171.     return "c_0_n[%i - 1] * " % y + m(x,y,y,y+1)
  172.  
  173. # define the final function
  174. def c_(u):
  175.     global c_n
  176.     c_n = ""
  177.     for el in range(1,u+1):
  178.         c_n = c_n + q(u,el) + " + "
  179.     if c_n.endswith(" + "):
  180.         c_n = c_n[:-3]
  181.     return c_n
  182.  
  183. # graphit subdefinition 1
  184. def reag():
  185.     global reag_var
  186.     while True:
  187.         try:
  188.             reag_var = map(int, raw_input("Provide the reagents to plot (separate with spacebar): ").split(" "))
  189.             if all(i in scope_n for i in reag_var):
  190.                 return
  191.             else:
  192.                 print
  193.                 print "INVALID INPUT! Provide integers between 1 and n."
  194.                 print
  195.                 continue
  196.         except (ValueError):
  197.             print
  198.             print "INVALID INPUT! Provide integers between 1 and n."
  199.             print
  200.     return
  201.  
  202. # graphit subdefinition 2
  203. def t_k():
  204.     global t_k_var
  205.     while True:
  206.         try:
  207.             t_k_var = deci.Decimal(raw_input("Define the time range from 0 to: "))
  208.             if (deci.Decimal(t_k_var) <= 0):
  209.                 print
  210.                 print "INVALID INPUT! Provide integer/float greater than 0."
  211.                 print
  212.                 continue
  213.             elif (deci.Decimal(t_k_var) > 0):
  214.                 return
  215.         except (ValueError):
  216.             print
  217.             print "INVALID INPUT! Provide integer/float greater than 0."
  218.             print
  219.     return
  220.  
  221. # graphit subdefinition 3
  222. def t_d():
  223.     global t_d_var
  224.     while True:
  225.         try:
  226.             t_d_var = deci.Decimal(raw_input("Define the precision of the time axis (e.g. 0.01): "))
  227.             break
  228.         except (ValueError):
  229.             print
  230.             print "INVALID INPUT! Provide integer/float."
  231.             print
  232.     return
  233.  
  234. # plot the requested functions
  235. def graphit():
  236.     if (n == 0):
  237.         scope()
  238.         conc()
  239.         kin()
  240.         summary()
  241.     elif (all(ele1 == 0 for ele1 in k_n) or all(ele2 == 0 for ele2 in c_0_n)):
  242.         conc()
  243.         kin()
  244.         summary()
  245.     try:
  246.         global t
  247.         reag()
  248.         t_k()
  249.         t_d()
  250.         t = np.arange(0,t_k_var,t_d_var)
  251.         p = []
  252.  
  253.         plt.xlabel("time t")
  254.         plt.ylabel("concentration c_n(t)")
  255.         for i in reag_var: # the actual plot command
  256.             p += plt.plot(t,eval(c_(i)),label= "c_" + str(i) + "(t)")
  257.             max_c = max(eval(c_(i)))
  258.             max_t = t[eval(c_(i)).argmax()]
  259.             print "t_%i_max =" % (i), max_t
  260.             print "c_%i_max =" % (i), max_c
  261.         plt.legend(loc="center right")
  262.         plt.show()
  263.     except (ZeroDivisionError):
  264.         print
  265.         print "INVALID INPUT! Kinetic constant k_n values cannot be equal to one another."
  266.         print
  267.     return
  268.  
  269. # interface below; stand-by mode
  270. while True:
  271.     command = raw_input("Provide command: ")
  272.     allowed = ["scope", "commands", "graphit", "kin", "conc", "summary", "test", "run", "manu"]
  273.     if command != "q":
  274.         comm = "%s()" % (command)
  275.         if command not in allowed:
  276.             print
  277.             print "INVALID COMMAND! Provide command from the list."
  278.             print
  279.         else:
  280.             exec(comm)
  281.     elif command == "q":
  282.         break