fix hiperbola offset

analisis/plots/20231123_CPTconmicromocion3/lolo_graficos_pub.py

@@ -36,12 +36,14 @@ from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
 
 # Funciones auxiliares
 
+
 from scipy.stats.distributions import  t,chi2
+import inspect
 
 def estadistica(datos_x,datos_y,modelo,pcov,parametros,nombres=None,alpha=0.05):
     
     if nombres is None:
-        nombres = [ f'{j}' for j in range(len(parametros)) ]
+        nombres = inspect.getfullargspec(modelo).args[1:]
     
     # Cantidad de parámetros
     P = len(parametros)
@@ -90,13 +92,27 @@ def estadistica(datos_x,datos_y,modelo,pcov,parametros,nombres=None,alpha=0.05):
     sT = t.ppf(1.0 - alpha/2.0, N - P ) # student T multiplier
     CI = sT * SE                        # Confidence Interval
     
-    print('R-squared    ',Rsq)
-    print('R-sq_adjusted',Rsq_adj)
-    print('chi2         ',chi2_)
+
+    return dict(R2=Rsq,R2_adj=Rsq_adj,chi2=chi2_,chi2_red=chi2_red,
+                chi2_test=chi2_test,r=r_pearson,pvalue=p_val,
+                SE=SE,CI=CI,parametros=parametros,nombres=nombres,alpha=alpha,N=N)
+
+
+def print_estadistica(dic):
+    
+    pars = 'N,pvalue,R2,R2_adj,chi2,chi2_red,chi2_test,r,SE,CI,parametros,nombres,alpha'.split(',') 
+    
+    N,pvalue,R2,R2_adj,chi2,chi2_red,chi2_test,r,SE,CI,parametros,nombres,alpha = [ dic[k] for k in pars  ]
+    P = len(parametros)
+    
+    
+    print('R-squared    ',R2)
+    print('R-sq_adjusted',R2_adj)
+    print('chi2         ',chi2)
     print('chi2_reduced ',chi2_red)
     print('chi2_test    ',chi2_test)
-    print('r-pearson    ',r_pearson)
-    print('p-value      ',p_val)
+    print('r-pearson    ',r)
+    print('p-value      ',pvalue)
     print('')
     print('Error Estandard (SE):')
     for i in range(P):
@@ -106,10 +122,6 @@ def estadistica(datos_x,datos_y,modelo,pcov,parametros,nombres=None,alpha=0.05):
     for i in range(P):
         print(f'parametro[{nombres[i]:>5s}]: ' , parametros[i], ' ± ' , CI[i])
 
-    return dict(R2=Rsq,R2_adj=Rsq_adj,chi2=chi2_,chi2_red=chi2_red,
-                chi2_test=chi2_test,r=r_pearson,pvalue=p_val,
-                SE=SE,CI=CI)
-
 
 
 
@@ -312,6 +324,15 @@ def hiperbola(x,a,y0,b,x0):
     """
     return a*np.sqrt(((x-x0)**2+b**2))+y0
 
+
+def hiperbola2(x,a,b,x0):
+    """
+    Hiperbola de ecuación:
+        1 =(y)²/a² - (x-x0)²/b² 
+    """
+    return a*np.sqrt(((x-x0)**2+b**2))
+
+
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #%% Graficamos todos los fiteos
 
@@ -440,11 +461,31 @@ par_inicial = (12,0.1,1,-0.13)
 param,pcov = curve_fit(hiperbola,voltages_dcA[I],Betas_vec[I],p0=par_inicial)
 x_hip = np.linspace(-0.23,0.005,200)
 
+EST=estadistica(voltages_dcA[I],Betas_vec[I],hiperbola,pcov,param,nombres=None,alpha=0.05)
+print_estadistica(EST)
+
+
+I = slice(None,9)
+par_inicial = (12,1,-0.13)
+param2,pcov2 = curve_fit(hiperbola2,voltages_dcA[I],Betas_vec[I],p0=par_inicial)
+x_hip = np.linspace(-0.23,0.005,200)
+
+EST2=estadistica(voltages_dcA[I],Betas_vec[I],hiperbola2,pcov2,param2,nombres=None,alpha=0.05)
+print_estadistica(EST2)
+
+
+
+
+
+
+
 
 ax = ax_central
 ax.errorbar(voltages_dcA[I],Betas_vec[I],yerr=ErrorBetas_vec[I],fmt='o',
             capsize=4,markersize=3,color='C3', label=r'fitted $\beta$')
 ax.plot(x_hip,hiperbola(x_hip,*param),color='C0', label=r'hyperbola model')
+ax.plot(x_hip,hiperbola2(x_hip,*param2),'--',color='C4', label=r'hyperbola model')
+
 ax.set_ylabel(r'Modulation factor $\beta$', labelpad=-5)
 ax.set_ylim(-0.05,3)
 ax.set_xlim(-0.22,0)
@@ -453,6 +494,8 @@ ax.set_title(f'(g)', x=0.95, y=0.006, color='gray')
 ax = ax_res
 ax.errorbar(voltages_dcA[I],Betas_vec[I]-hiperbola(voltages_dcA[I],*param),
             yerr=ErrorBetas_vec[I],fmt='o',capsize=4,markersize=3,color='C3')
+ax.errorbar(voltages_dcA[I]+0.005,Betas_vec[I]-hiperbola2(voltages_dcA[I],*param2),
+            yerr=ErrorBetas_vec[I],fmt='o',capsize=4,markersize=3,color='C4')
 ax.axhline( 0 , color='C0')
 ax.set_ylabel('Res.', labelpad=-5)
 ax.set_xlabel('Endcap voltage [V]')
@@ -536,7 +579,7 @@ gs2 = gridspec.GridSpec(4,3, width_ratios=width_ratios,left=0.16,right=0.9,botto
 
 ax_central = fig.add_subplot(gs2[1:-1, 1])
 
-ax_dib = fig.add_subplot(gs2[0, 1], sharex=ax_res)
+ax_dib = fig.add_subplot(gs2[0, 1], sharex=ax_res, zorder=-20)
 ax_dib.spines[:].set_visible(False)
 ax_dib.xaxis.set_tick_params(labelbottom=False)
 ax_dib.yaxis.set_tick_params(labelleft=False)
@@ -609,27 +652,30 @@ axx[1,2].set_ylabel('Counts')
 
 ## Grafico central ###############################
 I = slice(None,9)
-par_inicial = (12,0.1,1,-0.13)
-param,pcov = curve_fit(hiperbola,voltages_dcA[I],Betas_vec[I],p0=par_inicial)
+par_inicial = (12,1,-0.13)
+param,pcov = curve_fit(hiperbola2,voltages_dcA[I],Betas_vec[I],p0=par_inicial)
 x_hip = np.linspace(-0.23,0.005,200)
 
+EST=estadistica(voltages_dcA[I],Betas_vec[I],hiperbola2,pcov,param,nombres=None,alpha=0.05)
+print_estadistica(EST)
+
 
 ax = ax_central
 ax.errorbar(voltages_dcA[I],Betas_vec[I],yerr=ErrorBetas_vec[I],fmt='o',
             capsize=4,markersize=3,color='C3', label=r'fitted $\beta$')
-ax.plot(x_hip,hiperbola(x_hip,*param),color='C0', label=r'hyperbola model')
+ax.plot(x_hip,hiperbola2(x_hip,*param),color='C0', label=r'hyperbola model')
 ax.set_ylabel(r'Modulation factor $\beta$', labelpad=-5)
 ax.set_ylim(-0.05,3)
 ax.set_xlim(-0.22,0)
-ax.set_title(f'(g)', x=0.95, y=0.006, color='gray')
+ax.set_title(f'(h)', x=0.95, y=0.006, color='gray')
 
 ax = ax_res
-ax.errorbar(voltages_dcA[I],Betas_vec[I]-hiperbola(voltages_dcA[I],*param),
+ax.errorbar(voltages_dcA[I],Betas_vec[I]-hiperbola2(voltages_dcA[I],*param),
             yerr=ErrorBetas_vec[I],fmt='o',capsize=4,markersize=3,color='C3')
 ax.axhline( 0 , color='C0')
 ax.set_ylabel('Res.', labelpad=-5)
 ax.set_xlabel('Endcap voltage [V]')
-ax.set_title(f'(h)', x=0.95, y=0.72, color='gray')
+ax.set_title(f'(i)', x=0.95, y=0.72, color='gray')
 
 
 ax_res.get_xticklabels()[-1].set_visible(False)
@@ -644,19 +690,25 @@ for ax in [ax_central,ax_res]:
 # Anotaciones ##########################
 for jj,ax in zip(selection,axes_vec):
     
+    rad = 0.2 if jj==8 or jj==3 else -0.2
     Axes_x = 1 if axx.flatten().tolist().index(ax)%3==0 else 0
     ax_central.annotate("",
                          xy=(ax_central.get_lines()[0].get_xdata()[jj], ax_central.get_lines()[0].get_ydata()[jj]), 
                          xycoords=ax_central.transData,
                          xytext=(Axes_x, 0.5), textcoords=ax.transAxes,
-                         arrowprops=dict(arrowstyle="<-",connectionstyle="arc3,rad=-0.2", color='gray', alpha=0.5),
+                         arrowprops=dict(arrowstyle="<-",connectionstyle=f"arc3,rad={rad}", color='gray', alpha=0.5),
                          zorder=-2)
     
 
 
 # dibujos ##########################
 
-a, y0, b, x0 = param
+# a, y0, b, x0 = param
+a, b, x0 = param
+
+
+y0 = -0.212663964284234
+
 
 ax_dib.plot( x_hip , x_hip*0+y0 , color='gray')
 ax_dib.set_ylim(-y0*0.1,y0*1.2)
@@ -670,6 +722,8 @@ ax_dib.plot( voltages_dcA[I], voltages_dcA[I]*0+y0, '.' )
 for Vx in voltages_dcA[I]:
     ax_dib.plot( [Vx,x0], [y0,0], ':', color='lightgray', lw=1)
 
+
+
 import matplotlib.patches as mpatches
 
 arr = mpatches.FancyArrowPatch((x0, y0*1.1), (x0/2, y0*1.1),
@@ -678,6 +732,15 @@ ax_dib.add_patch(arr)
 ax_dib.annotate("ion position", (.5, .5), xycoords=arr, ha='center', va='bottom', color='C1')
 
 
+# r
+arr = mpatches.FancyArrowPatch((x0,0), ( voltages_dcA[5], y0),
+                               arrowstyle='->,head_width=.25', mutation_scale=10, color='black')
+ax_dib.add_patch(arr)
+ax_dib.annotate(r"$\vec r$", (.5, .5), xycoords=arr, ha='center', va='bottom', color='black')
+
+
+
+
 bbox = ax_dib.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
 aspect_ratio = bbox.width / bbox.height
 
@@ -693,6 +756,10 @@ ax_dib.text(0.1, 0.13, "x", transform=ax_dib.transAxes,
 ax_dib.text(0.03, 0.13*aspect_ratio, "y", transform=ax_dib.transAxes,
             fontsize=8,  va='top')
 
+
+ax_dib.set_title(f'(g)', x=0.95, y=0.5, color='gray')
+
+
 # Leyenda ##############################
 # h1, l1 = ax_central.get_legend_handles_labels()
 # h2, l2 = axx[0,0].get_legend_handles_labels()
@@ -922,7 +989,7 @@ ax.plot(betaslong,[t*1e3 for t in FinalTemp_withoutC(betaslong,*popt_rho22_balan
         '-', color='C0',  label='model w/ RF heating')
 
 ax.plot(betaslong,[t*1e3 for t in FinalTemp_withoutmicromotion(betaslong,*popt_rho22_balance_noMM)],
-        '--', color='C1',label='model w/ RF heating')
+        '--', color='C1',label='model w/o RF heating')
 
 
 ax.fill_between( [0, 2.8] , np.ones(2)*Doppler_limit_area[0],   np.ones(2)*Doppler_limit_area[1],