Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download
Views: 50
Kernel: Python 3 (Anaconda)
import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as plt from matplotlib.pyplot import *
#data1 tesla_1=open("1 tesla.txt","r") #import data data_1=[] for line in tesla_1: x,y=line.split() #split lines x,y=float(x),float(y) #convert to numbers data_1.append([x,y]) #store data to array td=[] pcd=[] for i,j in data_1: td.append(i) pcd.append(j) time=np.asarray(td) currrent=np.asarray(pcd) def intensity(t,a,tor,w,phi,c): return ((a*np.exp(-t/tor))*np.sin(w*t+phi)+c) g1=0.80 g2=0.00000000015 g3=-4*10**10 g4=np.pi/3 g5=0 xdata=time*10**(-12) ydata=intensity(xdata,g1,g2,g3,g4,g5) popt,pcov=curve_fit(intensity,xdata,pcd,p0=[g1,g2,g3,g4,g5]) #plot(xdata, intensity(xdata, popt[0], popt[1], popt[2],popt[3],popt[4])) #print(popt[0], popt[1], popt[2],popt[3],popt[4]) #plot(xdata, intensity(xdata, g1, g2, g3,g4,g5)) figure(figsize=[25,15]) plot(xdata,intensity(xdata, popt[0], popt[1], popt[2],popt[3],popt[4]),'k-.',label='Optimised Data') plot(xdata,pcd,'c',label='Original Data') plot(xdata,pcd,'r.',label='Original Data Dotted') plot(xdata,intensity(xdata,g1,g2,g3,g4,g5),'r',label='guess') ylabel("Intensity", fontsize=20) xlabel("time", fontsize=20) plt.axvline(x=0.5e-10,color='k') plt.axvline(x=1.75e-10,color='k') plt.axvline(x=3e-10,color='k') plt.axvline(x=4.25e-10,color='k') legend(loc=1) title("1 Tesla") print("Omega = (%.2f +/- %.2f)x10rad/s" % (popt[2]*10**(-10), np.sqrt(pcov[2, 2])*10**(-10)))
Omega = (-4.99 +/- 0.04)x10rad/s
Image in a Jupyter notebook
#data1 tesla_1=open("1 tesla.txt","r") #import data data_1=[] for line in tesla_1: x,y=line.split() #split lines x,y=float(x),float(y) #convert to numbers data_1.append([x,y]) #store data to array td=[] pcd=[] for i,j in data_1: td.append(i) pcd.append(j) time=np.asarray(td) currrent=np.asarray(pcd) def intensity(t,a,tor,w,phi,c): return ((a*np.exp(-t/tor))*np.sin(w*t+phi)+c) g1=0.80 g2=0.00000000015 g3=-4*10**10 g4=np.pi/3 g5=0 xdata=time*10**(-12)
#data2 tesla_2=open("2 tesla.txt","r") #import data data_2=[] for line in tesla_2: x2,y2=line.split() #split lines x2,y2=float(x2),float(y2) #convert to numbers data_2.append([x2,y2]) #store data to array td2=[] pcd2=[] for i,j in data_2: td2.append(i) pcd2.append(j) time2=np.asarray(td2) currrent2=np.asarray(pcd2) g1_2=1.2 g2_2=5.2e-10 g3_2=-10*10**10 g4_2=np.pi/4 g5_2=0 xdata2=time2*10**(-12)
#data3 tesla_3=open("3 tesla.txt","r") #import data data_3=[] for line in tesla_3: x3,y3=line.split() #split lines x3,y3=float(x3),float(y3) #convert to numbers data_3.append([x3,y3]) #store data to array td3=[] pcd3=[] for i,j in data_3: td3.append(i) pcd3.append(j) time3=np.asarray(td3) currrent3=np.asarray(pcd3) g1_3=-0.60 g2_3=4e-10 g3_3=-15*10**10 g4_3=np.pi g5_3=0 xdata3=time3*10**(-12)
#data4 tesla_4=open("4 tesla.txt","r") #import data data_4=[] for line in tesla_4: x4,y4=line.split() #split lines x4,y4=float(x4),float(y4) #convert to numbers data_4.append([x4,y4]) #store data to array td4=[] pcd4=[] for i,j in data_4: td4.append(i) pcd4.append(j) time4=np.asarray(td4) currrent4=np.asarray(pcd4) g1_4=-0.5 g2_4=3e-10 g3_4=-20*10**10 g4_4=-np.pi/2.5 g5_4=0 xdata4=time4*10**(-12)
def omega(xdata,pd,a,b,c,d,e): intensity1=intensity(xdata,a,b,c,d,e) popt,pcov=curve_fit(intensity,xdata,pd,p0=[(a,b,c,d,e)]) intensity2=intensity(xdata, popt[0], popt[1], popt[2],popt[3],popt[4]) figure(figsize=[20,10]) plot(xdata,intensity1,'r',label='guess') plot(xdata,intensity2,'k-.',label='Optimised Data') plot(xdata,pd,'c',label='Original Data') plot(xdata,pd,'r.',label='Original Data Dotted') ylabel("Intensity (A)", fontsize=20) xlabel("Time delay (s)", fontsize=20) print(("Omega = (%.3f +/- %.3f)x10^(10)rad/s" % (popt[2]*10**(-10), np.sqrt(pcov[2, 2])*10**(-10)))) legend(loc=1)
omega1=omega(xdata,pcd,g1,g2,g3,g4,g5) omega2=omega(xdata2,pcd2,g1_2,g2_2,g3_2,g4_2,g5_2) omega3=omega(xdata3,pcd3,g1_3,g2_3,g3_3,g4_3,g5_3) omega4=omega(xdata4[8:181],pcd4[8:181],g1_4,g2_4,g3_4,g4_4,g5_4)
Omega = (-4.991 +/- 0.038)x10^(10)rad/s Omega = (-9.942 +/- 0.052)x10^(10)rad/s Omega = (-14.815 +/- 0.119)x10^(10)rad/s Omega = (-20.119 +/- 0.068)x10^(10)rad/s
Image in a Jupyter notebookImage in a Jupyter notebookImage in a Jupyter notebookImage in a Jupyter notebook
plt.figure(figsize=[25,15]) plt.plot(xdata4,pcd4) plt.axvline(x=0.45e-10,color='k') plt.axvline(x=0.75e-10,color='k') plt.axvline(x=1.05e-10,color='k') plt.axvline(x=1.35e-10,color='k')
<matplotlib.lines.Line2D at 0x7ffbb3f51cc0>
Image in a Jupyter notebook
#for curve fit values y=np.asarray([4.99e10,9.94e10,14.82e10,20.12e10]) x=[8.782e10,1.756e11,2.634e11,3.513e11] plt.figure(figsize=[25,15]) plt.plot(x,y) ylabel("Omega (rad/s)", fontsize=20) xlabel(" eB/2m (rad/s)", fontsize=20) slope,intercept=np.polyfit(x,y,1) errorbar(x,y,yerr=[0.038e10,0.052e10,0.119e10,0.068e10],color='r', linewidth=(1)) print(slope)
0.572397962283
Image in a Jupyter notebook
#for other method y1=np.asarray([5.03e10,9.97e10,15.7e10,20.9e10]) x1=[8.782e10,1.756e11,2.634e11,3.513e11] plt.figure(figsize=[25,15]) plt.plot(x1,y1, color='c') ylabel("Omega (rad/s)", fontsize=20) xlabel(" eB/2m (rad/s)", fontsize=20) slope1,intercept1=np.polyfit(x1,y1,1) print(slope1) errorbar(x1,y1,yerr=[0.402e10,1.58e10,3.93e10,6.97e10],color='r', linewidth=(1))
0.607351657471
<Container object of 3 artists>
Image in a Jupyter notebook
np.fft?
from scipy.fftpack import fft xaxisvalues=xdata funcfft=fft(pcd) funcfft=np.abs(funcfft) #print(funcfft) plt.plot(1/xaxisvalues,funcfft) #plt.axvline(x=0.7e-10,color='k') plt.axvline(x=4.1e10,color='k') plt.show()
/ext/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:6: RuntimeWarning: divide by zero encountered in true_divide
Image in a Jupyter notebook
from scipy.fftpack import fft xaxisvalues=xdata3 funcfft=fft(pcd3) funcfft=np.abs(funcfft) #print(funcfft) plt.plot(xaxisvalues,funcfft) #plt.axvline(x=0.7e-10,color='k') #plt.axvline(x=11e10,color='k') plt.show()
Image in a Jupyter notebook