Shareddatachallenge1.ipynbOpen in CoCalc
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)))

#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)



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>
#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)
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
#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')
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>
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
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()