todo

analisis/plots/20231123_CPTconmicromocion3/CPT_plotter_20231123.py

@@ -1129,16 +1129,20 @@ de que en el punto optimo el ion no esta en el centro de la trampa
 sino que esta a una distancia d
 """
 def hiperbola(x,a,b,c,x0):
-    return a*np.sqrt(((x-x0)**2+c**2))+b
+    #b=0
+    return np.sqrt((a*(x-x0)**2+c**2))+b
 
 medfin = 11
 voltages_dcA = Voltages[0][1:medfin]
 
 
-popthip,pcovhip = curve_fit(hiperbola,voltages_dcA[:10],Betas_vec[:10],p0=(100,0.1,1,-0.15))
+popthip,pcovhip = curve_fit(hiperbola,voltages_dcA[:10],Betas_vec[:10],p0=(20,-0.2,2e-2,-0.13))
+errorhip = np.sqrt(np.diag(pcovhip))
+
 
 xhip = np.linspace(-0.23,0.005,200)
 
+
 plt.figure()
 plt.errorbar(voltages_dcA,Betas_vec[:medfin-1],yerr=ErrorBetas_vec[:medfin-1],fmt='o',capsize=5,markersize=5,color=paleta[1])
 plt.plot(xhip,hiperbola(xhip,*popthip))
@@ -1147,11 +1151,17 @@ plt.ylabel('Modulation factor')
 #plt.yscale('log')
 plt.grid()
 
+betacero = popthip[2]+popthip[1]
+error_betacero = np.sqrt((errorhip[1])**2+(errorhip[2])**2+2*pcovhip[1,2])
+print(f'Beta_0 = {betacero} pm {error_betacero}')
 
 
 
 #%%
 from scipy.special import jv
+import seaborn as sns
+paleta = sns.color_palette("rocket")
+
 
 """
 Por si no quiero correr todo de nuevo:
@@ -1170,19 +1180,6 @@ def expo(x,tau,A,B):
 def cuadratica(x,a,c):
     return a*(x**2)+c
 
-def InverseMicromotionSpectra(beta, A, det, x0, gamma, B):
-    ftrap=22.1
-    #gamma=30
-    P = ((jv(0, beta)**2)/((((det-x0)**2)+(0.5*gamma)**2)**2))*(-2*(det-x0))
-    i = 1
-    #print(P)
-    while i <= 2:
-        P = P + (-2*(det-x0))*((jv(i, beta))**2)/(((((det-x0)+i*ftrap)**2)+(0.5*gamma)**2)**2) + (-2*(det-x0))*(((jv(-i, beta))**2)/((((det-x0)-i*ftrap)**2)+(0.5*gamma)**2)**2)
-        i = i + 1
-        #print(P)
-    #return 1/(A*P+B)
-    return 1/(A*P+B)
-
 
 def MicromotionSpectra(beta,det, gamma):
     ftrap=2*np.pi*22.1e6
@@ -1194,13 +1191,13 @@ def MicromotionSpectra(beta,det, gamma):
         P = P + ((jv(i, beta))**2)/((((det)+i*ftrap)**2)+(0.5*gamma)**2) + ((jv(-i, beta))**2)/((((det)-i*ftrap)**2)+(0.5*gamma)**2) 
         i = i + 1
         #print(P)
-    return P
-
+    return (gamma/2*np.pi)*P
+    #return P
 
-def polynomial(x,a,b,c,d,e):
-    b=0
-    d=0
-    return a+b*x+c*x*x+d*x*x*x+e*x*x*x*x
+    def polynomial(x,a,b,c,d,e):
+        b=0
+        d=0
+        return a+b*x+c*x*x+d*x*x*x+e*x*x*x*x
 
 def InverseDerivMicromotionSpectra(beta, det, gamma):
     ftrap=2*np.pi*22.1e6
@@ -1213,12 +1210,14 @@ def InverseDerivMicromotionSpectra(beta, det, gamma):
         P = P + (-2*(det))*((jv(i, beta))**2)/(((((det)+i*ftrap)**2)+(0.5*gamma)**2)**2) + (-2*(det))*(((jv(-i, beta))**2)/((((det)-i*ftrap)**2)+(0.5*gamma)**2)**2)
         i = i + 1
         #print(P)
-    return 1/P
+    return 1/((gamma/2*np.pi)*P)
+    #return 1/P
 
 
 def FinalTemp(beta,det, C,D):
     """
-    En mK
+    Este deja libre un factor C para corregir la primera parte pero
+    si ajustamos con esto, C da casi 1 entonces es al pedo.
     """
     gamma = 23e6
     hbar = 1.05e-34
@@ -1228,45 +1227,44 @@ def FinalTemp(beta,det, C,D):
     #det = -12e6
     #print(MicromotionSpectra(beta,det,gamma))
     #D=0
+    #C=1
     return ((hbar/(kb))*C*MicromotionSpectra(beta,2*np.pi*det,2*np.pi*gamma)+D*(beta**2))*InverseDerivMicromotionSpectra(beta, 2*np.pi*det, 2*np.pi*gamma)
     #return (C*MicromotionSpectra(beta,det,gamma))*InverseDerivMicromotionSpectra(beta, det, gamma)
 
 
-
-def FinalTemp_fixedall(beta, det, D):
+def FinalTemp_withoutC(beta,det, D):
     """
-    En mK
+    Este anda barbaro y es el que hay que usar
     """
     gamma = 23e6
-    #det=-11
-    #D=-0.8
-    #C = 1.68656122e-03  
-    #D = 6.64227010e-02
-    #D=1e-29
     hbar = 1.05e-34
     kb = 1.38e-23
-    D = 2e-9
     C = 1
-    det = -11e6
+    #D=0.5e-27
+    #det = -12e6
     #print(MicromotionSpectra(beta,det,gamma))
-    return 1e3*((hbar/(kb))*C*MicromotionSpectra(beta,2*np.pi*det,2*np.pi*gamma)+D*beta**2)*InverseDerivMicromotionSpectra(beta, 2*np.pi*det, 2*np.pi*gamma)
+    #D=0
+    #C=1
+    return ((hbar/(kb))*C*MicromotionSpectra(beta,2*np.pi*det,2*np.pi*gamma)+D*(beta**2))*InverseDerivMicromotionSpectra(beta, 2*np.pi*det, 2*np.pi*gamma)
     #return (C*MicromotionSpectra(beta,det,gamma))*InverseDerivMicromotionSpectra(beta, det, gamma)
 
 
+
 """
 Temperatura vs beta con un ajuste exponencial
 """
 
-#popt_rho22_balance, pcov_rho22_balance = curve_fit(FinalTemp,list(np.arange(0,0.3,0.001))+list(Betas_vec[:9]), list(0.5e-3*np.ones(len(np.arange(0,0.3,0.001))))+[t for t in Temp_vec[:9]],p0=(-10e6, 1,1e-27), bounds=((-30e6,0,0),(0e6,2,1e-25))) #esto ajusta muy bien
-popt_rho22_balance, pcov_rho22_balance = curve_fit(FinalTemp,list(Betas_vec[:9]), [t for t in Temp_vec[:9]],p0=(-11e6,0.8,4e-28)) #esto ajusta muy bien
-#popt_rho22_balance, pcov_rho22_balance = curve_fit(FinalTemp_fixedall,list(Betas_vec[:9]), [t*1e3 for t in Temp_vec[:9]],p0=(-10,1e-20), bounds=((-30,0),(0,1e6))) #esto ajusta muy bien
+#popt_rho22_balance, pcov_rho22_balance = curve_fit(FinalTemp,list(Betas_vec[:9]), [t for t in Temp_vec[:9]],p0=(-11e6,0.8,4e-28)) #esto ajusta muy bien
+#popt_rho22_balance, pcov_rho22_balance = curve_fit(FinalTemp,list(Betas_vec[:10]), [t for t in Temp_vec[:10]],sigma=ErrorTemp_vec[:10],absolute_sigma=False,p0=(-10e6,0.8,5e-20)) #esto ajusta muy bien
+
+popt_rho22_balance, pcov_rho22_balance = curve_fit(FinalTemp_withoutC,list(Betas_vec[:10]), [t for t in Temp_vec[:10]],sigma=ErrorTemp_vec[:10],absolute_sigma=False,p0=(-10e6,5e-20)) #esto ajusta muy bien
 
 
 
 betaslong = np.arange(0,2.8,0.01)
 
-#print(f'Min temp predicted: {FinalTemp_fixedall(betaslong,*popt_rho22_balance)[0]}')
-print(f'Min temp predicted: {1e3*FinalTemp(betaslong,*popt_rho22_balance)[0]} mK')
+print(f'Min temp predicted: {FinalTemp_withoutC(betaslong,*popt_rho22_balance)[0]}')
+#print(f'Min temp predicted: {1e3*FinalTemp(betaslong,*popt_rho22_balance)[0]} mK')
 
 print(f'Detuning: {1e-6*popt_rho22_balance[0]} MHz')
 
@@ -1278,11 +1276,11 @@ k_plot = 9
 
 plt.figure()
 plt.errorbar(Betas_vec[:k_plot],[t*1e3 for t in Temp_vec[:k_plot]],xerr=ErrorBetas_vec[:k_plot], yerr=[t*1e3 for t in ErrorTemp_vec[:k_plot]],fmt='o',capsize=5,markersize=5,color=paleta[3])
-plt.plot(betaslong,[t*1e3 for t in FinalTemp(betaslong,*popt_rho22_balance)],label='Ajuste con espectro modulado')
+#plt.plot(betaslong,[t*1e3 for t in FinalTemp(betaslong,*popt_rho22_balance)],label='Ajuste con espectro modulado')
+plt.plot(betaslong,[t*1e3 for t in FinalTemp_withoutC(betaslong,*popt_rho22_balance)],label='Ajuste con espectro modulado')
 #plt.plot(betaslong,FinalTemp_fixedall(betaslong,*popt_rho22_balance),label='Ajuste con espectro modulado')
 # plt.xlim(-0.1,1.1)
 #plt.ylim(0,1)
-
 #plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
 #plt.axhline(0.538)
 plt.xlabel('Modulation factor')
@@ -1290,13 +1288,69 @@ plt.ylabel('Temperature (mK)')
 plt.title('With RF heating')
 plt.grid()
 
+#%%
+"""
+Convierto D a RF heating rate.
+Analizo como es la derivada de n (fonones) respecto a la temperatura
+y veo en que regimenes puedo convertir n a T facilmente y cuando no.
+"""
+
+
+hbar = 1.05e-34
+kb = 1.38e-23
+
+bigfactor = hbar*((2*np.pi/0.397e-6)**2)*2*np.pi*23e6*kb/6.6e-26
+
+Tdot = 2*(popt_rho22_balance[1]*bigfactor)/kb #en unidades de k/s
+ErrorTdot = 2*(np.sqrt(pcov_rho22_balance[1,1])*bigfactor)/kb
+
+
+print()
+
+wtrap = 2*np.pi*1e6
+
+
+Temp = 1e-6
+
+ndot_factor_vec = []
+Tempvec = np.logspace(-6,-2,1000)
+
+for Temp in Tempvec:
+    ndot_factor = (hbar*wtrap/kb)*np.exp(hbar*wtrap/(kb*Temp))/((Temp**2)*(np.exp(hbar*wtrap/(kb*Temp))-1)**2)
+    ndot_factor_vec.append(ndot_factor)
+
+plt.figure()
+plt.semilogx(Tempvec,ndot_factor_vec)
+plt.xlabel('Temperature (K)')
+plt.ylabel('dn/dT')
+plt.grid()
+
+print(ndot_factor*Tdot) 
+
+"""
+se ve que para temperaturas de mas de 0.2 mK es recontra constante esa funcion y no depende de T
+"""
+
+#%%
+"""
+Calculo el heating rate en unidades de fonones/s, para beta = 1. Para otro beta hay que multiplicar por beta**2
+"""
+
+Temp = 0.5e-3
+wtrap = 2*np.pi*1e6
+ndot_factor_dopplerlimit = (hbar*wtrap/kb)*np.exp(hbar*wtrap/(kb*Temp))/((Temp**2)*(np.exp(hbar*wtrap/(kb*Temp))-1)**2)
+
+ndot = ndot_factor_dopplerlimit*Tdot
+errorndot = ndot_factor_dopplerlimit*ErrorTdot
+
+print(ndot)
 
 #%%
 
 hbar = 1.05e-34
 kb = 1.38e-23
 
-beta_vector = np.arange(0,1.5,0.01)
+beta_vector = np.arange(0,2,0.01)
 
 comparison = []
 
@@ -1306,10 +1360,11 @@ for beta_target in beta_vector:
     comparison.append(factor2/factor1)
 
 plt.figure()
-plt.plot(beta_vector,comparison)
-plt.xlabel('Beta')
-plt.ylabel('Ratio RF heat/Normal heat')
-plt.axhline(1)
+plt.axhline(1, color=paleta[4],linestyle='dashed')
+plt.plot(beta_vector,comparison,color=paleta[1],linewidth=3)
+plt.xlabel('Modulation factor beta')
+plt.ylabel('RF heating/Recoil heating')
+plt.grid()
 
 
 #%%