todo

analisis/plots/20240312_RotationalDopplerShift_news/RDS_piezobeamsizes_2024_invariance.py

@@ -48,12 +48,9 @@ def SeeKeys(files):
 
 def Lorentzian( x, A, B, x0, gam,C):
     #C=0
-    return A * gam**2 / ( gam**2 + ( x - x0 )**2) + B - C*(x - x0)
+    return A * gam**1 / ( gam**2 + ( x - x0 )**2) + B - C*(x - x0)
+
 
-def Lorentzian2( x, A, B, x0, gam):
-    C=0
-    A=0
-    return A * gam**2 / ( gam**2 + ( x - x0 )**2) + B - C*(x - x0)
 
 
 Piezo1Counts = []
@@ -124,7 +121,7 @@ palette = sns.color_palette("tab10")
 
 pmlocmedvec1=list(np.arange(0,len(PIEZOS1_FILES)-2,1))
 
-#pmlocmedvec1 = []
+#pmlocmedvec1 = [1]
 
 plt.figure()
 
@@ -148,7 +145,7 @@ for med in pmlocmedvec1:
     Counts = [c for c in Piezo1Counts[med][1:]]
     
     if med==1:
-        popt, pcov = curve_fit(Lorentzian, np.array(Freqs[0:100]), np.array(Counts[0:100]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10)))
+        popt, pcov = curve_fit(Lorentzian, np.array(Freqs[0:100]), np.array(Counts[0:100]), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,1000)))
     
     else:
         popt, pcov = curve_fit(Lorentzian, np.array(Freqs), np.array(Counts), p0=(-2000,2100,435.8,0.05,0.1), bounds=((-100000,0,435.5,0,0),(0,1e5, 436.1, 10,10000)))
@@ -184,7 +181,7 @@ if len(pmlocmedvec1)!=1:
     plt.figure()
     plt.errorbar(np.arange(0,len(Intensityver1),1), [i/np.max(Intensityver1) for i in Intensityver1], yerr=np.sqrt(Intensityver1)/np.max(Intensityver1), fmt='-o', capsize=3,markersize=8)
     plt.errorbar(np.arange(0,len(Intensityver1),1), [p for p in pmdepthsdrver1], yerr= errorpmdepthsdrver1, fmt='o', capsize=3, markersize=8)
-    plt.errorbar(np.arange(0,len(Intensityver1),1), [5*p for p in Gamas1], yerr=ErrorGamas1, fmt='o', capsize=3, markersize=8)
+    plt.errorbar(np.arange(0,len(Intensityver1),1), [5*p for p in Gamas1], yerr=[5*e for e in ErrorGamas1], fmt='o', capsize=3, markersize=8)
     plt.xlabel('Ion position')
     plt.ylabel('Intensity / DR Relative depth')
     #plt.xticks([1,2,3,4,5])
@@ -294,7 +291,7 @@ if len(pmlocmedvec2)!=1:
     plt.figure()
     plt.errorbar(np.arange(0,len(Intensityver2),1), [i/np.max(Intensityver2) for i in Intensityver2], yerr=np.sqrt(Intensityver2)/np.max(Intensityver2), fmt='-o', capsize=3,markersize=8)
     plt.errorbar(np.arange(0,len(Intensityver2),1), [p for p in pmdepthsdrver2], yerr= errorpmdepthsdrver2, fmt='o', capsize=3, markersize=8)
-    plt.errorbar(np.arange(0,len(Intensityver2),1), [5*p for p in Gamas2], yerr=ErrorGamas2, fmt='o', capsize=3, markersize=8)
+    plt.errorbar(np.arange(0,len(Intensityver2),1), [5*p for p in Gamas2], yerr=[5*e for e in ErrorGamas2], fmt='o', capsize=3, markersize=8)
     plt.xlabel('Ion position (steps)')
     plt.ylabel('DR Relative depth')
     #plt.xticks([1,2,3,4,5])
@@ -359,6 +356,10 @@ plt.ylabel('DR depth / Beam intensity')
 Intento ajustar con un modelo
 
 F(r)=A*( r²/(r²+b²)
+        
+y con un modelo tipo
+
+F(r) = A*(J_0 (b/r))²
 
 """
 from scipy.optimize import curve_fit
@@ -391,29 +392,24 @@ plt.legend()
 
 plt.grid()
 
-#%%
-"""
-Va otro modelo con funciones de bessel
 
-F(r) = 1-(J_0 (b*r))²
-"""
 
-from scipy.special import jv
 
+from scipy.special import jv
 
 
-def modelo2(r,a,b):
-    return a*(1-(jv(0, b*(r))**2))  
+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(modelo2,rfit,yfit,p0=(0.7,0.005))
+popt,pcov = curve_fit(modelo3,rfit,yfit,p0=(0.7,1))
 
 print(popt)
-plt.figure()
+#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)
@@ -421,20 +417,17 @@ plt.plot([x*20 for x in xchicofinal],[i/np.max(IntensityChico) for i in Intensit
 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.plot(rlong,modelo3(rlong,*popt),zorder=2,color='green',linewidth=3,linestyle='dashed',label=r'$I(r)=a\,(J_0(\frac{b}{r}))^2$')
 
 plt.xlim(-20,600+x2+10)
-
-plt.legend(loc='upper left')
+plt.legend()
 plt.grid()
 
-
 #%%
 """
-Va otro modelo con funciones de bessel
-
+Veo el modelo 
 
-F(r) = A*(J_0 (b/r))²
+F(r)=1/((a²+b/(r²)))
 
 """
 
@@ -442,7 +435,7 @@ from scipy.special import jv
 
 
 def modelo3(r,a,b):
-    return a*((jv(0, b*(1/r)))**2)  
+    return 1/((a**2+b/(r**2)))
 
 rfit = np.arange(x2,len(Intensityver2)*20+x2,20)
 yfit = [p for p in pmdepthsdrver2]
@@ -469,8 +462,49 @@ plt.grid()
 
 
 
+#%%
+"""
+Intento ajustar con un modelo con bessels
+Uso los primeros tres puntos de la medicion con haz chico y el resto con el haz grande
+"""
+from scipy.optimize import curve_fit
+
 
+def modelo3(r,a,b,d):
+    return (a*jv(0, b*(1/(r-d)))*jv(0,2*b*(1/(r-d))))**2
 
+rfit = np.arange(x1,len(Intensityver2)*20+x2,20)
+yfit = [pmdepthsdrver1[0],pmdepthsdrver1[1],pmdepthsdrver1[2]] + [p for p in pmdepthsdrver2]
+
+popt,pcov = curve_fit(modelo3,rfit,yfit)
+
+
+rfit = np.arange(x2,len(Intensityver2)*20+x2,20)
+yfit = [p for p in pmdepthsdrver2]
+
+rlong = np.arange(1e-10,1000,1)
+
+popt,pcov = curve_fit(modelo3,rfit,yfit,p0=(10,10,10))
+
+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='green',linewidth=3,label=r'$I(r)=a\,(J_0(\frac{b}{r}))^2$')
+
+plt.xlim(-20,600+x2+10)
+plt.legend()
+plt.grid()
+
+# plt.figure()
+# plt.plot(rlong,popt[1]/rlong)
+# plt.xlim(0,200)
+# plt.ylim(-10,30)