Want more features on Pastebin? Sign Up, it's FREE!
Guest

Untitled

By: a guest on Mar 3rd, 2013  |  syntax: Python  |  size: 7.43 KB  |  views: 4  |  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. 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
clone this paste RAW Paste Data