• API
• FAQ
• Tools
• Archive
SHARE
TWEET

Planks constant

a guest Feb 14th, 2020 77 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import numpy as np
2. import matplotlib.pyplot as plt
3. from scipy.optimize import least_squares
4.
5.
6.
7. def fitFunc(p, x):
8.     '''
9.     Fit function
10.     '''
11.     f = (p[0]+ p[1]*x) + p[2]*(np.exp(p[3]*(x-p[4])) - 1)
12.     return f
13.
14. def fitFuncDiff(p, x):
15.     '''
16.     Differential of fit function
17.     '''
18.     df = p[1] + p[3]*p[2]*np.exp(p[3]*(x-p[4]))
19.     return df
20.
21. def calcChiSq(p, x, y, xerr, yerr):
22.     '''
23.     Error function for fit
24.     '''
25.     e = (y - fitFunc(p, x))/(np.sqrt(yerr**2 + fitFuncDiff(p, x)**2*xerr**2))
26.     return e
27.
28. def fitStdError(jacMatrix):
29.
30.     # Compute covariance
31.     jMat2 = np.dot(jacMatrix.T, jacMatrix)
32.     detJmat2 = np.linalg.det(jMat2)
33.
34.     # Prepare output
35.     output = np.zeros(jMat2.shape[0])
36.     if detJmat2 < 1E-32:
37.         print("Value of determinat detJmat2",detJmat2)
38.         print("Matrix singular, error calculation failed.")
39.         return output
40.     else:
41.         covar = np.linalg.inv(jMat2)
42.         for i in range(len(output)):
43.             output[i] = np.sqrt(covar[i, i])
44.
45.         return output
46.
47. planck_data = np.genfromtxt("Planck.csv", dtype=None, delimiter=',')
48.
49. V1_405nm = planck_data[:, 0]
50. V1_405nm = V1_405nm[~np.isnan(V1_405nm)]
51.
52. i1_405nm = planck_data[:, 1]
53. i1_405nm = i1_405nm[~np.isnan(i1_405nm)]
54.
55. V1_436nm = planck_data[:, 2]
56. V1_436nm = V1_436nm[~np.isnan(V1_436nm)]
57.
58. i1_436nm = planck_data[:, 3]
59. i1_436nm = i1_436nm[~np.isnan(i1_436nm)]
60.
61. V1_546nm = planck_data[:, 4]
62. V1_546nm = V1_546nm[~np.isnan(V1_546nm)]
63.
64. i1_546nm = planck_data[:, 5]
65. i1_546nm = i1_546nm[~np.isnan(i1_546nm)]
66.
67. V1_578nm = planck_data[:, 6]
68. V1_578nm = V1_578nm[~np.isnan(V1_578nm)]
69.
70. i1_578nm = planck_data[:, 7]
71. i1_578nm = i1_578nm[~np.isnan(i1_578nm)]
72.
73. V1_365nm = planck_data[:, 8]
74. V1_365nm = V1_365nm[~np.isnan(V1_365nm)]
75.
76. i1_365nm = planck_data[:, 9]
77. i1_365nm = i1_365nm[~np.isnan(i1_365nm)]
78.
79.
80. yerror = i1_405nm*0.01+0.005*6
81. xerror = V1_405nm*0.01+0.005*6
82.
83.
84.
85.
86.
87. # Set initial values of fit parameters
88. #          0         1        2
89. #          A         b         c      d         v0
90. pInit =   [0,        0,       10,     0.1,      -1.25]
91. lBounds = [-10,      -10,     0,     0,      -1.5]
92. uBounds = [10,       10,       100,      10,      -1]
93. nPoints = len(V1_405nm)
94. nPars = 5
95.
96. # Run fit
97. output = least_squares(calcChiSq, pInit, args = (V1_405nm, i1_405nm, xerror, yerror),
98.                 bounds = (lBounds, uBounds))
99.
100.
101.
102. # Get least_squares output, stored in array output.x[]
103. A = output.x[0]
104. b = output.x[1]
105. c = output.x[2]
106. d = output.x[3]
107. V0 = output.x[4]
108.
109. # Get errors from our fits using fitStdError(), defined above
110. pErrors = fitStdError(output.jac)
111. d_A = pErrors[0]
112. d_b = pErrors[1]
113. d_c = pErrors[2]
114. d_d = pErrors[3]
115. d_V0 = pErrors[4]
116. # Output fit parameters
117. print("Fitted parameters a: {0:.2f}, b: {1:.2f}, c: {2:.2f}, d: {3:.2f}, V0: {4:.2f}".format(A, b, c, d ,V0))
118. print("Parameter errors: a: {0:.2f}, b: {1:.2f}, c: {2:.2f},d: {3:.2f}, V0: {4:.2f}" .format(d_A, d_b, d_c, d_d, d-V0))
119.
120.
121.
122. # Calculate chis**2 per point, summed chi**2 and chi**2/NDF
123. chiarr = calcChiSq(output.x, V1_405nm, i1_405nm, xerror, yerror)**2
124. chisq = np.sum(chiarr)
125. NDF = nPoints - nPars
126. chisqndf = chisq/NDF
127. print("ChiSq = {:5.2e}, ChiSq/NDF = {:5.2f}.".format(chisq, chisqndf))
128.
129. # Calculate fitted y-values using our fit parameters and the original fit function
130. xPlot = np.linspace(0.9*np.min(V1_405nm), 1.1*np.max(V1_405nm), 100)
131. fitData = fitFunc(output.x, xPlot)
132.
133.
134.
135.
136. fig = plt.figure(figsize = (8, 6))
137. plt.title('405nM')
138. plt.xlabel('V')
139. plt.ylabel('I')
140. plt.grid(color = 'g')
141. plt.errorbar(V1_405nm, i1_405nm, yerror, xerror, marker='.', linestyle='', label='data')
142. plt.plot(xPlot, fitData, color = 'r')
143. plt.show()
144.
145.
146. ###################################################################################
147. del yerror
148. del xerror
149.
150. yerror = i1_436nm*0.01+0.005*6
151. xerror = V1_436nm*0.01+0.005*6
152.
153. pInit =   [0,        0,       10,     0.1,      -1.25]
154. lBounds = [-10,      -10,     0,     0,      -1.5]
155. uBounds = [10,       10,       100,      10,      -1]
156. nPoints = len(V1_436nm)
157. nPars = 5
158.
159.
160.
161. # Run fit
162. output = least_squares(calcChiSq, pInit, args = (V1_436nm, i1_436nm, xerror, yerror),
163.                     bounds = (lBounds, uBounds))
164.
165.
166.
167. # Get least_squares output, stored in array output.x[]
168. A = output.x[0]
169. b = output.x[1]
170. c = output.x[2]
171. d = output.x[3]
172. V0 = output.x[4]
173.
174. # Get errors from our fits using fitStdError(), defined above
175. pErrors = fitStdError(output.jac)
176. d_A = pErrors[0]
177. d_b = pErrors[1]
178. d_c = pErrors[2]
179. d_d = pErrors[3]
180. d_V0 = pErrors[4]
181. # Output fit parameters
182. print("Fitted parameters a: {0:.2f}, b: {1:.2f}, c: {2:.2f}, d: {3:.2f}, V0: {4:.2f}".format(A, b, c, d ,V0))
183. print("Parameter errors: a: {0:.2f}, b: {1:.2f}, c: {2:.2f},d: {3:.2f}, V0: {4:.2f}" .format(d_A, d_b, d_c, d_d, d-V0))
184.
185.
186.
187. # Calculate chis**2 per point, summed chi**2 and chi**2/NDF
188. chiarr = calcChiSq(output.x, V1_436nm, i1_436nm, xerror, yerror)**2
189. chisq = np.sum(chiarr)
190. NDF = nPoints - nPars
191. chisqndf = chisq/NDF
192. print("ChiSq = {:5.2e}, ChiSq/NDF = {:5.2f}.".format(chisq, chisqndf))
193.
194. # Calculate fitted y-values using our fit parameters and the original fit function
195. xPlot = np.linspace(0.9*np.min(V1_436nm), 1.1*np.max(V1_436nm), 100)
196. fitData = fitFunc(output.x, xPlot)
197.
198.
199.
200.
201. fig = plt.figure(figsize = (8, 6))
202. plt.title('436 nM')
203. plt.xlabel('x')
204. plt.ylabel('y')
205. plt.grid(color = 'g')
206. plt.errorbar(V1_436nm, i1_436nm, yerror, xerror, marker='.', linestyle='', label='data')
207. plt.plot(xPlot, fitData, color = 'r')
208. plt.show()
209.
210.
211. ##############################################################################################################
212.
213. del yerror
214. del xerror
215.
216. yerror = i1_546nm*0.01+0.005*6
217. xerror = V1_546nm*0.01+0.005*6
218.
219. # Set initial values of fit parameters
220. #          0         1        2
221. #          A         b         c      d         v0
222. pInit =   [0,        0,       10,     0.1,      -1.25]
223. lBounds = [-10,      -10,     0,     0,      -1.5]
224. uBounds = [10,       10,       100,      10,      -1]
225. nPoints = len(V1_546nm)
226. nPars = 5
227.
228. # Run fit
229. output = least_squares(calcChiSq, pInit, args = (V1_546nm, i1_546nm, xerror, yerror),
230.                 bounds = (lBounds, uBounds))
231.
232.
233.
234. # Get least_squares output, stored in array output.x[]
235. A = output.x[0]
236. b = output.x[1]
237. c = output.x[2]
238. d = output.x[3]
239. V0 = output.x[4]
240.
241. # Get errors from our fits using fitStdError(), defined above
242. pErrors = fitStdError(output.jac)
243. d_A = pErrors[0]
244. d_b = pErrors[1]
245. d_c = pErrors[2]
246. d_d = pErrors[3]
247. d_V0 = pErrors[4]
248. # Output fit parameters
249. print("Fitted parameters a: {0:.2f}, b: {1:.2f}, c: {2:.2f}, d: {3:.2f}, V0: {4:.2f}".format(A, b, c, d ,V0))
250. print("Parameter errors: a: {0:.2f}, b: {1:.2f}, c: {2:.2f},d: {3:.2f}, V0: {4:.2f}" .format(d_A, d_b, d_c, d_d, d-V0))
251.
252.
253.
254. # Calculate chis**2 per point, summed chi**2 and chi**2/NDF
255. chiarr = calcChiSq(output.x, V1_546nm, i1_546nm, xerror, yerror)**2
256. chisq = np.sum(chiarr)
257. NDF = nPoints - nPars
258. chisqndf = chisq/NDF
259. print("ChiSq = {:5.2e}, ChiSq/NDF = {:5.2f}.".format(chisq, chisqndf))
260.
261. # Calculate fitted y-values using our fit parameters and the original fit function
262. xPlot = np.linspace(0.9*np.min(V1_546nm), 1.1*np.max(V1_546nm), 100)
263. fitData = fitFunc(output.x, xPlot)
264.
265.
266.
267.
268. fig = plt.figure(figsize = (8, 6))
269. plt.title('546nM')
270. plt.xlabel('V')
271. plt.ylabel('I')
272. plt.grid(color = 'g')
273. plt.errorbar(V1_546nm, i1_546nm, yerror, xerror, marker='.', linestyle='', label='data')
274. plt.plot(xPlot, fitData, color = 'r')
275. plt.show()
276.
277.
278.  ###########################################################################################################
279.
280. del yerror
281. del xerror
282.
283. yerror = i1_578nm*0.01+0.005*6
284. xerror = V1_578nm*0.01+0.005*6
285.
286. # Set initial values of fit parameters
287. #          0         1        2
288. #          A         b         c      d         v0
289. pInit =   [0,        0,       10,     0.1,      -1.25]
290. lBounds = [-10,      -10,     0,     0,      -1.5]
291. uBounds = [10,       10,       100,      10,      -1]
292. nPoints = len(V1_578nm)
293. nPars = 5
294.
295. # Run fit
296. output = least_squares(calcChiSq, pInit, args = (V1_578nm, i1_578nm, xerror, yerror),
297.                 bounds = (lBounds, uBounds))
298.
299.
300.
301. # Get least_squares output, stored in array output.x[]
302. A = output.x[0]
303. b = output.x[1]
304. c = output.x[2]
305. d = output.x[3]
306. V0 = output.x[4]
307.
308. # Get errors from our fits using fitStdError(), defined above
309. pErrors = fitStdError(output.jac)
310. d_A = pErrors[0]
311. d_b = pErrors[1]
312. d_c = pErrors[2]
313. d_d = pErrors[3]
314. d_V0 = pErrors[4]
315. # Output fit parameters
316. print("Fitted parameters a: {0:.2f}, b: {1:.2f}, c: {2:.2f}, d: {3:.2f}, V0: {4:.2f}".format(A, b, c, d ,V0))
317. print("Parameter errors: a: {0:.2f}, b: {1:.2f}, c: {2:.2f},d: {3:.2f}, V0: {4:.2f}" .format(d_A, d_b, d_c, d_d, d-V0))
318.
319.
320.
321. # Calculate chis**2 per point, summed chi**2 and chi**2/NDF
322. chiarr = calcChiSq(output.x, V1_578nm, i1_578nm, xerror, yerror)**2
323. chisq = np.sum(chiarr)
324. NDF = nPoints - nPars
325. chisqndf = chisq/NDF
326. print("ChiSq = {:5.2e}, ChiSq/NDF = {:5.2f}.".format(chisq, chisqndf))
327.
328. # Calculate fitted y-values using our fit parameters and the original fit function
329. xPlot = np.linspace(0.9*np.min(V1_578nm), 1.1*np.max(V1_578nm), 100)
330. fitData = fitFunc(output.x, xPlot)
331.
332.
333.
334.
335. fig = plt.figure(figsize = (8, 6))
336. plt.title('578nM')
337. plt.xlabel('V')
338. plt.ylabel('I')
339. plt.grid(color = 'g')
340. plt.errorbar(V1_578nm, i1_578nm, yerror, xerror, marker='.', linestyle='', label='data')
341. plt.plot(xPlot, fitData, color = 'r')
342. plt.show()
343.
344. ##############################################################################
345.
346. del yerror
347. del xerror
348.
349. yerror = i1_365nm*0.01+0.005*6
350. xerror = V1_365nm*0.01+0.005*6
351.
352. # Set initial values of fit parameters
353. #          0         1        2
354. #          A         b         c      d         v0
355. pInit =   [0,        0,       10,     0.1,      -1.25]
356. lBounds = [-10,      -10,     0,     0,      -1.5]
357. uBounds = [10,       10,       100,      10,      -1]
358. nPoints = len(V1_365nm)
359. nPars = 5
360.
361. # Run fit
362. output = least_squares(calcChiSq, pInit, args = (V1_365nm, i1_365nm, xerror, yerror),
363.                 bounds = (lBounds, uBounds))
364.
365.
366.
367. # Get least_squares output, stored in array output.x[]
368. A = output.x[0]
369. b = output.x[1]
370. c = output.x[2]
371. d = output.x[3]
372. V0 = output.x[4]
373.
374. # Get errors from our fits using fitStdError(), defined above
375. pErrors = fitStdError(output.jac)
376. d_A = pErrors[0]
377. d_b = pErrors[1]
378. d_c = pErrors[2]
379. d_d = pErrors[3]
380. d_V0 = pErrors[4]
381. # Output fit parameters
382. print("Fitted parameters a: {0:.2f}, b: {1:.2f}, c: {2:.2f}, d: {3:.2f}, V0: {4:.2f}".format(A, b, c, d ,V0))
383. print("Parameter errors: a: {0:.2f}, b: {1:.2f}, c: {2:.2f},d: {3:.2f}, V0: {4:.2f}" .format(d_A, d_b, d_c, d_d, d-V0))
384.
385.
386.
387. # Calculate chis**2 per point, summed chi**2 and chi**2/NDF
388. chiarr = calcChiSq(output.x, V1_365nm, i1_365nm, xerror, yerror)**2
389. chisq = np.sum(chiarr)
390. NDF = nPoints - nPars
391. chisqndf = chisq/NDF
392. print("ChiSq = {:5.2e}, ChiSq/NDF = {:5.2f}.".format(chisq, chisqndf))
393.
394. # Calculate fitted y-values using our fit parameters and the original fit function
395. xPlot = np.linspace(0.9*np.min(V1_365nm), 1.1*np.max(V1_365nm), 100)
396. fitData = fitFunc(output.x, xPlot)
397.
398.
399.
400.
401. fig = plt.figure(figsize = (8, 6))
402. plt.title('365nM')
403. plt.xlabel('V')
404. plt.ylabel('I')
405. plt.grid(color = 'g')
406. plt.errorbar(V1_365nm, i1_365nm, yerror, xerror, marker='.', linestyle='', label='data')
407. plt.plot(xPlot, fitData, color = 'r')
408. plt.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.
Top