• API
• FAQ
• Tools
• Trends
• Archive
daily pastebin goal
51%
SHARE
TWEET

# Untitled

a guest Mar 3rd, 2013 6 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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
RAW Paste Data
Top