Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Untitled

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