# -*- coding: utf-8 -*- """ Created on Sun Mar 24 11:40:42 2019 @author: Marnik Metting van Rijn """ import numpy as np import uncertainties as unc import uncertainties.unumpy as unp import matplotlib.pyplot as plt import scipy.stats from scipy.stats import chi2 def fit_polynomial(x, y, degree, uncertainties ): dof = x.size - degree fit, cov = np.polyfit(x, y, degree, cov=True) residuals = np.sum((y - np.polyval(fit, x))**2 / uncertainties**2) chisq = residuals / (dof) return fit, cov, residuals, chisq, dof def evaluate_chisq(chisq, dof): return 1-chi2.cdf(chisq, dof) if __name__ =='__main__': edef = 1.60217656535 * 10 **-19 #Determination of Ctot omega1 = unc.ufloat(290, 3) * 10**3 #wavetek shows 289, fehler sehr gross, variert c1 = unc.ufloat(5365, 7) * 10**-12 omega2 = unc.ufloat(127, 3) * 10**3 #wavetek shows 125 Ctot = omega2**2*c1/(omega1**2-omega2**2) L = 1./(omega1**2*Ctot) ##################################################### #Determination of R at 60dB rho = unc.ufloat(33, 0.5) * 10**3 A1in = unc.ufloat(1.78, 0.02) * 10**-3 A1out = unc.ufloat(1.41, 0.02) A1 = A1out/A1in U1 = unp.uarray(6.76, .02) * 10**-3 #signal of function generator U2 = unp.uarray(np.array([2.1, 2.08, 2.06, 2.05, 2.04, 2.04, 2.02, 2.01, 2, 1.98, 1.97, 1.96, 1.95, 1.94, 1.93, 1.92, 1.91, 1.89, 1.86, 1.86, 1.84, 1.83, 1.82]), .01) /A1 Im = unp.uarray(np.arange(8, 31, 1), 1) * 10**-3 #diode current R = U2/(U1-U2) * rho Rerror = unp.std_devs(R) fitr, covr, residualr, chisqr, dofr = fit_polynomial(unp.nominal_values(Im), unp.nominal_values(R), 1, unp.std_devs(R)) pr = evaluate_chisq(chisqr, dofr) #Linear Regression for R mr, qr = np.polyfit(unp.nominal_values(Im), unp.nominal_values(R), 1) Yr = unp.nominal_values(Im)*mr + qr ######################################################## #Determination of m, G_I constant, rectifier amplifier calibration Vt = unp.uarray([.06, .1, .12, .14, .16, .18, .2, .22, .24, .26, .28, .3, .32, .34, .36, .38], .02) * 10**-3 It = unp.uarray(np.array([1, 4, 5, 7, 9, 11, 13, 16, 19, 22, 26, 30, 35, 40, 44, 49]) , 1) * 10**-6 mabs, q = np.polyfit(unp.nominal_values(Vt)**2, unp.nominal_values(It), 1) fitm, covm, residualsm, chisqm, dofm = fit_polynomial(unp.nominal_values(Vt)**2, unp.nominal_values(It), 1, unp.std_devs(It)) pm = evaluate_chisq(chisqm, dofm) Y = unp.nominal_values(Vt)**2*mabs + q m = unc.ufloat(fitm[0], np.abs(fitm[1])) ######################################################### I = unp.uarray(np.arange(10, 31, 1), 1) * 10**-3 #input current left Iz = unp.uarray(np.array([6, 6, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 12])-1, 1) * 10**-6 #output current right with correction of 1, since the rectifier shows 1 microamp when the amplifier is connected and no current in noise generator (also outpout voltage less than 0.1 microVolt out of the noise generator) def Rf(I): return I * mr + qr def df(I): return 1./(4. * Rf(I) * Ctot) e = Iz / (m * Rf(unp.nominal_values(I))**2 * 2. *I * df(unp.nominal_values(I))) def plotr(): plt.figure() plt.plot(unp.nominal_values(Im)*1000, Yr/1000, '--',label ='linear Regression') plt.errorbar(unp.nominal_values(Im)*1000, unp.nominal_values(R)/1000, Rerror/1000, fmt = '+', capsize = 2, label = 'measurement') plt.rc('text', usetex=True) plt.rc('font', family = 'serif') plt.xlabel(r'Diode current $\langle I(t) \rangle$ in [mA]', fontsize = 13) plt.ylabel('Resistance $R$ in [K$\Omega$]', fontsize = 13) plt.legend(prop={'size': 11}) plt.grid() plt.savefig('ResMes') plt.show() def plotrh(): pullr = (unp.nominal_values(R) - Yr)/unp.std_devs(R) plt.hist(pullr, 6, density=1, histtype='stepfilled', facecolor = '#99bbff', alpha = 0.75) x = np.linspace(-3.0, 3.0, 50) # always label the axes, also for histograms plt.xlabel(r'pull', fontsize = 13) plt.ylabel(r'count (normalized)', fontsize= 13) plt.savefig('ResMesH') plt.show() def plotg(): plt.figure() plt.plot(unp.nominal_values(Vt)**2 * 10**6, Y*10**6, '--', label ='linear Regression') plt.errorbar(unp.nominal_values(Vt)**2 * 10**6, unp.nominal_values(It)*10**6, unp.std_devs(It)*10**6, unp.std_devs(Vt)**2*10**6, fmt = '+', capsize = 3, label = 'measurement') plt.legend(prop = {'size': 11}) plt.xlabel(r'Voltage $\langle U_{amp}^{2}(t) \rangle$ in $[nV]^{2}$', fontsize = 13) plt.ylabel('Current $I_z$ shown on the rectifier in [$\mu$A] ', fontsize = 13) plt.grid() plt.savefig('GaMes') plt.show() def plotgh(): pullg = (unp.nominal_values(It)-Y)/unp.std_devs(It) plt.hist(pullg, 5, density=1, histtype='stepfilled', facecolor='#99bbff', alpha=0.75) #by using sturges rule k = 1 + 3.322 log(N), N = size(It) x = np.linspace(-3.0, 3.0, 50) # always label the axes, also for histograms plt.xlabel(r'pull', fontsize = 13) plt.ylabel(r'count (normalized)', fontsize= 13) plt.savefig('GaMesH') plt.show() def plote(): plt.figure() plt.plot(unp.nominal_values(I), np.abs((1 - unp.nominal_values(e)/edef)*100), 'x', label = 'Deviation of the calculated elementary charge to literature value') plt.legend() plt.xlabel('Measurement, diode current [mA]') plt.ylabel('Deviation in $\%$') plt.axis([0.01, 0.03, 0, 100]) plt.grid() plt.show() def ploteh(): pulle = (unp.nominal_values(e)- edef)/unp.std_devs(e) plt.hist(pulle, 6, density = 1, histtype = 'stepfilled', facecolor = '#99bbff', alpha = 0.75) x = np.linspace(-3.0, 3.0, 50) # always label the axes, also for histograms plt.xlabel(r'pull', fontsize = 13) plt.ylabel(r'count (normalized)', fontsize= 13) plt.savefig('EMesH') plt.show() def plotex(): pulle = (unp.nominal_values(e)- edef)/unp.std_devs(e) plt.hist(pulle, 6, density = 1, histtype = 'stepfilled') x = np.linspace(0, 3.0, 100) plt.plot(x,scipy.stats.lognorm.pdf(x, 1)) plt.show()