-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfit_thermal_state.py
116 lines (89 loc) · 3.68 KB
/
fit_thermal_state.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 6 14:38:26 2022
@author: ions
"""
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from statsmodels.stats.proportion import proportion_confint
import math
class Fit_thermal_state:
def __init__(self,t,pop,bounds = [] ,p0 = [] , num_shots = 100, sideband = 'BSB'):
self.t = t
self.pop = pop
self.bounds = bounds
self.p0 = p0
self.num_shots = num_shots
self.fit = None
self.fit_cov = None
self.sideband = sideband
self.get_std(self.pop, self.num_shots)
def thermal_state_bsb(self,t,nbar,Ω,η,γ):
Pb = 0
for n in range(800):
Pn = (1/(nbar+1))*(nbar/(1+nbar))**n
Ωnnp1 = np.sqrt(n+1)*η*Ω
γn = γ*np.sqrt(n+1)
Pb+= Pn*np.cos(Ωnnp1*t)*np.exp(-γn*t)
return 1-0.5*(1+Pb)
def thermal_state_rsb(self,t,nbar,Ω,η,γ):
Pb = 0
for n in range(800):
Pn = (1/(nbar+1))*(nbar/(1+nbar))**n
Ωnnp1 = np.sqrt(n+1)*η*Ω
γn = γ*np.sqrt(n+1)
Pb+= Pn*np.cos(Ωnnp1*t)*np.exp(-γn*t)
return 1-0.5*(1+Pb)
def fit_thermal_state(self):
bounds = self.bounds
p0 = self.p0
t = self.t
pop = self.pop
if self.sideband == 'BSB':
if len(bounds) >0 and len(p0)>0:
fit_bsb, fit_bsb_cov = curve_fit(self.thermal_state_bsb,t,pop, bounds = bounds, p0 = p0,sigma = self.std)
elif len(bounds)>0 and len(p0)==0:
fit_bsb, fit_bsb_cov = curve_fit(self.thermal_state_bsb,t,pop, bounds = bounds,sigma = self.std)
else:
fit_bsb, fit_bsb_cov = curve_fit(self.thermal_state_bsb,t,pop,sigma = self.std)
else:
if len(bounds) >0 and len(p0)>0:
fit_bsb, fit_bsb_cov = curve_fit(self.thermal_state_rsb,t,pop, bounds = bounds, p0 = p0,sigma = self.std)
elif len(bounds)>0 and len(p0)==0:
fit_bsb, fit_bsb_cov = curve_fit(self.thermal_state_rsb,t,pop, bounds = bounds,sigma = self.std)
else:
fit_bsb, fit_bsb_cov = curve_fit(self.thermal_state_rsb,t,pop,sigma = self.std)
self.fit = fit_bsb
self.fit_cov = fit_bsb_cov
return fit_bsb, fit_bsb_cov
def plot_fit(self,color):
if self.fit:
fit = self.fit
fit_cov = self.fit_cov
else:
self.fit_thermal_state()
fit = self.fit
fit_cov = self.fit_cov
print(self.fit)
t = self.t
pop = self.pop
std = self.std
t_fit = np.linspace(min(t),max(t),1000)
plt.errorbar(t*1e6,pop,std,marker = '.', ls = 'none',c= color)
fit_std = np.sqrt(np.diag(fit_cov))[0]
if self.sideband =='BSB':
plt.plot(t_fit*1e6,self.thermal_state_bsb(t_fit,*self.fit),label = r'nbar:{:.1f}$({:.0f})$'.format(self.fit[0],10*fit_std),c = color)
else:
plt.plot(t_fit*1e6,self.thermal_state_rsb(t_fit,*self.fit),label = r'nbar:{:.1f}$({:.0f})$'.format(self.fit[0],10*fit_std),c = color)
plt.xlabel('BSB pulse duration (μs)')
plt.ylabel(r'$P_{\uparrow}$',fontsize = 14)
plt.legend();
def get_std(self,prob,N):
std = []
for p in prob:
k = p*N
confint = proportion_confint(k, N, alpha=0.3173, method='beta')
uncertainty = (confint[1]-confint[0])/2
std.append(uncertainty)
self.std = std