todo

analisis/plots/20240312_RotationalDopplerShift_news/RDS_piezobeamsizes_2024_invariance.py

@@ -122,9 +122,9 @@ power=1
 
 palette = sns.color_palette("tab10")
 
-pmlocmedvec1=list(np.arange(0,len(PIEZOS1_FILES),1))
+pmlocmedvec1=list(np.arange(0,len(PIEZOS1_FILES)-1,1))
 
-#pmlocmedvec1 = [6]
+#pmlocmedvec1 = []
 
 plt.figure()
 
@@ -219,7 +219,7 @@ palette = sns.color_palette("tab10")
 
 pmlocmedvec2=list(np.arange(0,len(PIEZOS2_FILES),1))
 
-#pmlocmedvec2 = [0,1,2]
+#pmlocmedvec2 = [4]
 
 plt.figure()
 
@@ -243,6 +243,8 @@ for med in pmlocmedvec2:
     if med==3:    
         popt, pcov = curve_fit(Lorentzian, np.array(Freqs[0:70]), np.array(Counts[0:70]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10000)))
 
+    if med==4:
+        popt, pcov = curve_fit(Lorentzian, np.array(Freqs[0:-4]), np.array(Counts[0:-4]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10000)))
     # elif med==7:
     #     k=60
     #     popt, pcov = curve_fit(Lorentzian, np.array(Freqs[0:10]+Freqs[k:]), np.array(Counts[0:10]+Counts[k:]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10000)))
@@ -250,8 +252,11 @@ for med in pmlocmedvec2:
     # elif med==10:
     #     popt, pcov = curve_fit(Lorentzian, np.array(Freqs[6:80]), np.array(Counts[6:80]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10000)))
 
-    # elif med==12:
-    #     popt, pcov = curve_fit(Lorentzian, np.array(Freqs[:75]), np.array(Counts[:75]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10000)))
+    elif med==25:
+        print('che')
+        k1,k2=0,40
+        k3,k4=-20,-1
+        popt, pcov = curve_fit(Lorentzian, np.array(Freqs[k1:k2]+Freqs[k3:k4]), np.array(Counts[k1:k2]+Counts[k3:k4]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10000)))
 
         
     else:
@@ -347,6 +352,127 @@ plt.grid()
 
 plt.xlabel('Distance to center of beam (steps)')
 plt.ylabel('DR depth / Beam intensity')
-plt.savefig('/home/nico/Documents/artiq_experiments/analisis/plots/20240312_RotationalDopplerShift_news/scaleinvariance.jpg')
+#plt.savefig('/home/nico/Documents/artiq_experiments/analisis/plots/20240312_RotationalDopplerShift_news/scaleinvariance.jpg')
+
+#%%
+"""
+Intento ajustar con un modelo
+
+F(r)=A*( r²/(r²+b²)
+
+"""
+from scipy.optimize import curve_fit
+
+
+def modelo1(r,a,b):
+    k = 2
+    return a*(r**k/(b**k+r**k))
+
+rfit = np.arange(x2,len(Intensityver2)*20+x2,20)
+yfit = [p for p in pmdepthsdrver2]
+
+rlong = np.arange(0,1000,1)
+
+popt,pcov = curve_fit(modelo1,rfit,yfit)
+
+
+plt.figure()
+
+plt.errorbar(np.arange(x1,len(Intensityver1)*20+x1,20), [p for p in pmdepthsdrver1], yerr= errorpmdepthsdrver1, fmt='o',color='red', capsize=3, markersize=8,zorder=1)
+plt.plot([x*20 for x in xchicofinal],[i/np.max(IntensityChico) for i in IntensityChico],'o',color='red',alpha=0.3)
+
+plt.errorbar(np.arange(x2,len(Intensityver2)*20+x2,20), [p for p in pmdepthsdrver2], yerr= errorpmdepthsdrver2, fmt='o',color='blue', capsize=3, markersize=8,zorder=1)
+plt.plot([x*20 for x in xgrandefinal],[i/np.max(IntensityGrande) for i in IntensityGrande],'o',color='blue',alpha=0.3)
+
+plt.plot(rlong,modelo1(rlong,*popt),zorder=2,color='black',linewidth=3,linestyle='dashed',label=r'$I(r)=A\,\frac{r^2}{r^2+b^2}$')
+
+plt.xlim(-20,600+x2+10)
+plt.legend()
+
+plt.grid()
+
+#%%
+"""
+Va otro modelo con funciones de bessel
+
+F(r) = 1-(J_0 (b*r))²
+"""
+
+from scipy.special import jv
+
+
+
+def modelo2(r,a,b):
+    return a*(1-(jv(0, b*(r))**2))  
+
+rfit = np.arange(x2,len(Intensityver2)*20+x2,20)
+yfit = [p for p in pmdepthsdrver2]
+
+rlong = np.arange(0,1000,1)
+
+popt,pcov = curve_fit(modelo2,rfit,yfit,p0=(0.7,0.005))
+
+print(popt)
+plt.figure()
+
+plt.errorbar(np.arange(x1,len(Intensityver1)*20+x1,20), [p for p in pmdepthsdrver1], yerr= errorpmdepthsdrver1, fmt='o',color='red', capsize=3, markersize=8,zorder=1)
+plt.plot([x*20 for x in xchicofinal],[i/np.max(IntensityChico) for i in IntensityChico],'o',color='red',alpha=0.3)
+
+plt.errorbar(np.arange(x2,len(Intensityver2)*20+x2,20), [p for p in pmdepthsdrver2], yerr= errorpmdepthsdrver2, fmt='o',color='blue', capsize=3, markersize=8,zorder=1)
+plt.plot([x*20 for x in xgrandefinal],[i/np.max(IntensityGrande) for i in IntensityGrande],'o',color='blue',alpha=0.3)
+
+plt.plot(rlong,modelo2(rlong,*popt),zorder=2,color='black',linewidth=3,linestyle='dashed',label=r'$I(r)=a\,(1-(J_0(br)))$')
+
+plt.xlim(-20,600+x2+10)
+
+plt.legend(loc='upper left')
+plt.grid()
+
+
+#%%
+"""
+Va otro modelo con funciones de bessel
+
+
+F(r) = A*(J_0 (b/r))²
+
+"""
+
+from scipy.special import jv
+
+
+def modelo3(r,a,b):
+    return a*((jv(0, b*(1/r)))**2)  
+
+rfit = np.arange(x2,len(Intensityver2)*20+x2,20)
+yfit = [p for p in pmdepthsdrver2]
+
+rlong = np.arange(0,1000,1)
+
+popt,pcov = curve_fit(modelo3,rfit,yfit,p0=(0.7,1))
+
+print(popt)
+plt.figure()
+
+plt.errorbar(np.arange(x1,len(Intensityver1)*20+x1,20), [p for p in pmdepthsdrver1], yerr= errorpmdepthsdrver1, fmt='o',color='red', capsize=3, markersize=8,zorder=1)
+plt.plot([x*20 for x in xchicofinal],[i/np.max(IntensityChico) for i in IntensityChico],'o',color='red',alpha=0.3)
+
+plt.errorbar(np.arange(x2,len(Intensityver2)*20+x2,20), [p for p in pmdepthsdrver2], yerr= errorpmdepthsdrver2, fmt='o',color='blue', capsize=3, markersize=8,zorder=1)
+plt.plot([x*20 for x in xgrandefinal],[i/np.max(IntensityGrande) for i in IntensityGrande],'o',color='blue',alpha=0.3)
+
+plt.plot(rlong,modelo3(rlong,*popt),zorder=2,color='black',linewidth=3,linestyle='dashed',label=r'$I(r)=a\,(J_0(\frac{b}{r}))^2$')
+
+plt.xlim(-20,600+x2+10)
+plt.legend()
+plt.grid()
+
+
+
+
+
+
+
+
+