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. OK, I Understand
Top