API tools faq
paste
Advertisement
Guest User

Planks constant

a guest
Feb 14th, 2020
115
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 12.14 KB | None | 0 0
raw download clone embed print report
  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()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!