lolo_analisis_superajuste.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Ploteo de datos y ajustes
@author: lolo
"""


import h5py
import matplotlib.pyplot as plt
import numpy as np
# import sys
# import re
# import ast
from scipy.optimize import curve_fit
# import os
# from scipy import interpolate

from numba import jit,njit

from time import time

#%% Importaciones extra

# /home/lolo/Dropbox/marce/LIAF/Trampa_anular/artiq_experiments/analisis/plots/20231123_CPTconmicromocion3/Data/EITfit/MM_eightLevel_2repumps_AnalysisFunctions.py

from Data.EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM


# PARAMETROS = np.load('PARAMETROS.npz',allow_pickle=True)
# for var_name in PARAMETROS.keys():
#     globals()[var_name] = PARAMETROS[var_name]
#     print(f'loaded: {var_name}')


# Funciones auxiliares

from scipy.stats.distributions import  t,chi2

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)) ]
    
    # Cantidad de parámetros
    P = len(parametros)
    
    # Número de datos
    N = len(datos_x)
    
    # Grados de libertas (Degrees Of Freedom)
    dof = N-P-1

    # Cauculamos coordenadas del modelo
    # modelo_x    = datos_x if modelo_x_arr is None else modelo_x_arr
    # modelo_y    = modelo( modelo_x, *parametros )
    
    # Predicción del modelo para los datos_x medidos
    prediccion_modelo = modelo( datos_x, *parametros )
    
    # Calculos de cantidades estadísticas relevantes
    COV       = pcov                                      # Matriz de Covarianza
    SE        = np.sqrt(np.diag( COV  ))                        # Standar Error / Error estandar de los parámetros
    residuos  = datos_y - prediccion_modelo               # diferencia enrte el modelo y los datos
    
    SSE       = sum(( residuos )**2 )                     # Resitual Sum of Squares
    SST       = sum(( datos_y - np.mean(datos_y))**2)        # Total Sum of Squares
    
    # http://en.wikipedia.org/wiki/Coefficient_of_determination
    # Expresa el porcentaje de la varianza que logra explicar el modelos propuesto
    Rsq       =  1 - SSE/SST                               # Coeficiente de determinación
    Rsq_adj   = 1-(1-Rsq) * (N-1)/(N-P-1)                  # Coeficiente de determinación Ajustado   
    
    # https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#In_least_squares_regression_analysis
    # Expresa la correlación que hay entre los datos y la predicción del modelo
    r_pearson = np.corrcoef( datos_y ,  prediccion_modelo )[0,1]
    
    # Reduced chi squared
    # https://en.wikipedia.org/wiki/Reduced_chi-squared_statistic
    chi2_     = sum( residuos**2 )/N
    chi2_red  = sum( residuos**2 )/(N-P)
    
    # Chi squared test
    chi2_test = sum( residuos**2 / abs(prediccion_modelo) )
    # p-value del ajuste
    p_val  = chi2(dof).cdf( chi2_test )
    
    
    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_)
    print('chi2_reduced ',chi2_red)
    print('chi2_test    ',chi2_test)
    print('r-pearson    ',r_pearson)
    print('p-value      ',p_val)
    print('')
    print('Error Estandard (SE):')
    for i in range(P):
        print(f'parametro[{nombres[i]:>5s}]: ' , parametros[i], ' ± ' , SE[i])
    print('')
    print('Intervalo de confianza al '+str((1-alpha)*100)+'%:')
    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)






#%%

"""
Primero tengo mediciones de espectros cpt de un ion variando la tension dc_A
"""

#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20220106_CPT_DosLaseres_v08_TISA_DR\Data
# os.chdir('../20231123_CPTconmicromocion3/Data/')

folder = '../20231123_CPTconmicromocion3/Data/'
CPT_FILES = f"""
{folder}/000016262-IR_Scan_withcal_optimized
{folder}/000016239-IR_Scan_withcal_optimized
{folder}/000016240-IR_Scan_withcal_optimized
{folder}/000016241-IR_Scan_withcal_optimized
{folder}/000016244-IR_Scan_withcal_optimized
{folder}/000016255-IR_Scan_withcal_optimized
{folder}/000016256-IR_Scan_withcal_optimized
{folder}/000016257-IR_Scan_withcal_optimized
"""


def SeeKeys(files):
    for i, fname in enumerate(files.split()):
        data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
        print(fname)
        print(list(data['datasets'].keys()))

print(SeeKeys(CPT_FILES))

#carpeta pc nico labo escritorio:
#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20211101_CPT_DosLaseres_v03\Data

Counts = []
Freqs = []

AmpTisa = []
UVCPTAmp = []
No_measures = []
Voltages = []

for i, fname in enumerate(CPT_FILES.split()):
    print(str(i) + ' - ' + fname)
    #print(fname)
    data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'

    # Aca hago algo repugnante para poder levantar los strings que dejamos
    # que además tenian un error de tipeo al final. Esto no deberá ser necesario
    # cuando se solucione el error este del guardado.
    Freqs.append(np.array(data['datasets']['IR1_Frequencies']))
    Counts.append(np.array(data['datasets']['data_array']))
    #AmpTisa.append(np.array(data['datasets']['TISA_CPT_amp']))
    UVCPTAmp.append(np.array(data['datasets']['UV_CPT_amp']))
    No_measures.append(np.array(data['datasets']['no_measures']))
    Voltages.append(np.array(data['datasets']['scanning_voltages']))

def Split(array,n):
    length=len(array)/n
    splitlist = []
    jj = 0
    while jj<length:
        partial = []
        ii = 0
        while ii < n:
            partial.append(array[jj*n+ii])
            ii = ii + 1
        splitlist.append(partial)
        jj = jj + 1
    return splitlist


CountsSplit = []
CountsSplit.append(Split(Counts[0],len(Freqs[0])))


CountsSplit_2ions = []
CountsSplit_2ions.append(Split(Counts[4],len(Freqs[4])))


#%% Cargo parámetros fiteados de antes


PARAMETROS = np.load('analisis_superajuste_PARAMETROS.npz',allow_pickle=True)
for var_name in PARAMETROS.keys():
    globals()[var_name] = PARAMETROS[var_name]
    print(f'loaded: {var_name}')


if False:
    # Esto es para correr en caso de necesidad de limpiar todos los vectores de parametros
    print('Limpio los vectores de parámetros')
    for var in [ kk for kk in globals().keys() if kk.endswith('_vec') ]:
        print(f'del {var}')
        del(globals()[var])


#%% Definiciones de Numba

@jit
def FitEIT_MM_single_plot(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP):
    #BETA = 1.8
    # SG = 0.6
    # SP = 8.1
    # TEMP = 0.2e-3

    freqs = [2*f*1e-6-offset for f in Freqs]

    Detunings, Fluorescence1 = PerformExperiment_8levels_MM(SG, SP, gPS, gPD, DetDoppler, u, 
                                                            DopplerLaserLinewidth, ProbeLaserLinewidth, 
                                                            TEMP, alpha, phidoppler, titadoppler, 
                                                            phiprobe, titaprobe,  BETA1, drivefreq, 
                                                            min(freqs), max(freqs)+(freqs[1]-freqs[0]), 
                                                            freqs[1]-freqs[0], circularityprobe=CircPr, 
                                                            plot=False, solvemode=1, detpvec=None)

    ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
    return ScaledFluo1, Detunings



@jit
def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP):
    "Esta verison de la función devuelve sólo el eje y, para usar de modelo en un ajuste"
    return FitEIT_MM_single_plot(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP)[0]



param_names = 'offset DetDoppler SG SP SCALE1 OFFSET BETA1 TEMP'.split()

#%%
"""
AHORA INTENTO SUPER AJUSTES O SEA CON OFFSETXPI Y DETDOPPLER INCLUIDOS

La 0 no ajusta bien incluso con todos los parametros libres
De la 1 a la 11 ajustan bien
"""


#from EITfit.MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels
from scipy.optimize import curve_fit

"""
SUPER AJUSTE (SA)
"""

phidoppler, titadoppler = 0, 90
phirepump, titarepump = 0,  0
phiprobe = 0
titaprobe = 90

Temp = 0.5e-3

sg = 0.544
sp = 4.5
sr = 0
DetRepump = 0


lw = 0.1
DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth = lw, lw, lw #ancho de linea de los laseres


u = 32.5e6

#B = (u/(2*np.pi))/c

correccion = 13

#DetDoppler = -11.5-correccion

gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
alpha = 0


drivefreq = 2*np.pi*22.135*1e6


SelectedCurveVec = [1,2,3,4,5,6,7,8,9,10,11]
#SelectedCurveVec = [10]

CircPr = 1
alpha = 0

t0 = time()

if not 'popt_SA_vec' in globals().keys() or len(popt_SA_vec)==0:

    popt_SA_vec = []
    pcov_SA_vec = []
    Detuningsshort_vec = []
    Counts_vec = []
    Detuningslong_vec = []
    FittedCounts_vec = []

    Betas_vec = []
    ErrorBetas_vec = []
    Temp_vec = []
    ErrorTemp_vec = []

    DetuningsUV_vec = []
    ErrorDetuningsUV_vec = []
    
    Estadistica_vec = []

    for selectedcurve in SelectedCurveVec:
        print(f"{round(time()-t0,1):6.1f}: Procesando la curva {selectedcurve}")
        #selectedcurve = 2 #IMPORTANTE: SELECCIONA LA MEDICION

        FreqsDR = Freqs[0]
        CountsDR = CountsSplit[0][selectedcurve]


        if selectedcurve==1:
            CountsDR[100]=0.5*(CountsDR[99]+CountsDR[101])
            CountsDR[105]=0.5*(CountsDR[104]+CountsDR[106])
        if selectedcurve==2:
            CountsDR[67]=0.5*(CountsDR[66]+CountsDR[68])
            CountsDR[71]=0.5*(CountsDR[70]+CountsDR[72])
        if selectedcurve==6:
            CountsDR[1]=0.5*(CountsDR[0]+CountsDR[2])
            CountsDR[76]=0.5*(CountsDR[75]+CountsDR[77])
        if selectedcurve==7:
            CountsDR[117]=0.5*(CountsDR[116]+CountsDR[118])

        freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))


        if True:
            popt_3_SA, pcov_3_SA = curve_fit(FitEIT_MM_single, FreqsDR, CountsDR, p0=[430, -25, 0.9, 6.2, 3e4, 1.34e3, 2, (np.pi**2)*1e-3], bounds=((0, -50, 0, 0, 0, 0, 0, 0), (1000, 0, 2, 20, 5e4, 5e4, 10, (np.pi**2)*10e-3)))

            popt_SA_vec.append(popt_3_SA)
            pcov_SA_vec.append(pcov_3_SA)

            FittedEITpi_3_SA_short, Detunings_3_SA_short = FitEIT_MM_single_plot(FreqsDR, *popt_3_SA)
            freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
            FittedEITpi_3_SA_long, Detunings_3_SA_long = FitEIT_MM_single_plot(freqslong, *popt_3_SA)
            
            
            # estadistica(datos_x,datos_y,modelo,pcov,parametros,nombres=None,alpha=0.05)
            est_tmp = estadistica(FreqsDR,CountsDR,FitEIT_MM_single,pcov_3_SA,popt_3_SA,
                                  nombres=param_names,alpha=0.05)
            Estadistica_vec.append(est_tmp)

        DetuningsUV_vec.append(popt_3_SA[1])
        ErrorDetuningsUV_vec.append(np.sqrt(pcov_3_SA[1,1]))

        Betas_vec.append(popt_3_SA[6])
        ErrorBetas_vec.append(np.sqrt(pcov_3_SA[6,6]))
        Temp_vec.append(popt_3_SA[7])
        ErrorTemp_vec.append(np.sqrt(pcov_3_SA[7,7]))

        Detuningsshort_vec.append(Detunings_3_SA_short)
        Counts_vec.append(CountsDR)
        Detuningslong_vec.append(Detunings_3_SA_long)
        FittedCounts_vec.append(FittedEITpi_3_SA_long)


#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#%% Graficamos todos los fiteos

# tmp_datos=(Detuningsshort_vec,Counts_vec,Detuningslong_vec,FittedCounts_vec,SelectedCurveVec)

# for Detunings_3_SA_short,CountsDR,Detunings_3_SA_long,FittedEITpi_3_SA_long,selectedcurve in zip(*tmp_datos):
#     plt.figure()
#     plt.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
#     plt.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='darkolivegreen', linewidth=3, label=f'med {selectedcurve}')
#     #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
#     plt.xlabel('Detuning (MHz)')
#     plt.ylabel('Counts')
#     plt.legend(loc='upper left', fontsize=20)
#     plt.grid()

#     print(f'listo med {selectedcurve}')
#     print(popt_3_SA)


fig, axx = plt.subplots( 3,4, figsize=(13,8) ,  constrained_layout=True, sharex=True , sharey=True )
fig.set_constrained_layout_pads(w_pad=2/72, h_pad=2/72, hspace=0, wspace=0)

tmp_datos=(Detuningsshort_vec,Counts_vec,Detuningslong_vec,FittedCounts_vec,SelectedCurveVec,axx.flatten())

for Detunings_3_SA_short,CountsDR,Detunings_3_SA_long,FittedEITpi_3_SA_long,selectedcurve,ax in zip(*tmp_datos):
    ax.errorbar(Detunings_3_SA_short, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.3, capsize=2, markersize=2)
    ax.plot(Detunings_3_SA_long, FittedEITpi_3_SA_long, color='black', linewidth=2, label=f'med {selectedcurve}', alpha=0.7)
    #plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {round(popt[2]*1e3, 2)} mK, detDop: {DetDoppler} MHz')
    # ax.set_xlabel('Detuning (MHz)')
    # ax.set_ylabel('Counts')
    ax.legend(loc='upper left', fontsize=12)
    ax.grid(True, ls=":")

    print(f'listo med {selectedcurve}')
    print(popt_3_SA)


for ax in axx[:,0]:
    ax.set_ylabel('Counts')

for ax in axx[-1,:]:
    ax.set_xlabel('Detuning (MHz)')




#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#%% Inspección de parámetros

param_names = 'offset DetDoppler SG SP SCALE1 OFFSET BETA1 TEMP'.split()

err_vecs = np.array([ np.sqrt(np.diag(el)) for el in pcov_SA_vec ])
num_med  = np.arange(len(pcov_SA_vec)) +1

r2_values = np.array([ el['R2_adj'] for el in Estadistica_vec ])


fig, axx = plt.subplots( len(popt_SA_vec[0])+1,1, figsize=(13,8) ,  constrained_layout=True, sharex=True , sharey=False )
fig.set_constrained_layout_pads(w_pad=2/72, h_pad=2/72, hspace=0, wspace=0)

for ax,param_vec,err_vec,par_name in zip(axx,popt_SA_vec.T,err_vecs.T,param_names) :
    ax.plot(num_med, param_vec, '.-')
    ax.errorbar( num_med, param_vec, yerr=err_vec, 
                 fmt='s', mfc='none', elinewidth = 1, capsize=3, ms=1)

    ax.grid(True, ls=":", color='lightgray')
    ax.set_ylabel(par_name)


ax=axx[-1]
ax.plot( num_med , r2_values, '.-')
ax.set_ylabel(r'$R^2$')
ax.grid(True, ls=":", color='lightgray')



fig.align_ylabels()
ax.set_xticks(num_med)
ax.set_xlabel('Num. de medición')


#%%
"""
Grafico distintas variables que salieron del SUper ajuste
"""

import seaborn as sns
paleta = sns.color_palette("rocket")

voltages_dcA = Voltages[0][1:10]

def lineal(x,a,b):
    return a*x+b

def hiperbola(x,a,b,c,x0):
    return a*np.sqrt(((x-x0)**2+c**2))+b

hiperbola_or_linear = True

if hiperbola_or_linear:
    popthip,pcovhip = curve_fit(hiperbola,voltages_dcA,Betas_vec,p0=(100,0.1,1,-0.15))

    xhip = np.linspace(-0.23,0.005,200)

    plt.figure()
    plt.errorbar(voltages_dcA,Betas_vec,yerr=ErrorBetas_vec,fmt='o',capsize=5,markersize=5,color=paleta[1])
    plt.plot(xhip,hiperbola(xhip,*popthip))
    plt.xlabel('Endcap voltage (V)')
    plt.ylabel('Modulation factor')
    plt.grid()

else:
    poptini,pcovini = curve_fit(lineal,voltages_dcA[0:3],Betas_vec[0:3])
    poptfin,pcovfin = curve_fit(lineal,voltages_dcA[4:],Betas_vec[4:])

    minimum_voltage = -(poptini[1]-poptfin[1])/(poptini[0]-poptfin[0]) #voltaje donde se intersectan las rectas, es decir, donde deberia estar el minimo de micromocion
    minimum_modulationfactor = lineal(minimum_voltage,*poptini) #es lo mismo si pongo *poptfin

    xini = np.linspace(-0.23,-0.13,100)
    xfin = np.linspace(-0.15,0.005,100)

    plt.figure()
    plt.errorbar(voltages_dcA,Betas_vec,yerr=ErrorBetas_vec,fmt='o',capsize=5,markersize=5,color=paleta[1])
    plt.plot(xini,lineal(xini,*poptini))
    plt.plot(xfin,lineal(xfin,*poptfin))
    plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
    plt.xlabel('Endcap voltage (V)')
    plt.ylabel('Modulation factor')
    plt.grid()


print([t*1e3 for t in Temp_vec])

plt.figure()
plt.errorbar(voltages_dcA,[t*1e3 for t in Temp_vec],yerr=[t*1e3 for t in ErrorTemp_vec],fmt='o',capsize=5,markersize=5,color=paleta[3])
# plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
print(f'\n\nTE FALTA DEFINIR LA VARIABLE minimum_voltage\n\n')
plt.axhline(0.538)
plt.xlabel('Endcap voltage (V)')
plt.ylabel('Temperature (mK)')
plt.grid()
#plt.ylim(0,2)

#%%
"""
Ahora hago un ajuste con una hiperbola porque tiene mas sentido, por el hecho
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

popthip,pcovhip = curve_fit(hiperbola,voltages_dcA,Betas_vec,p0=(100,0.1,1,-0.15))

xhip = np.linspace(-0.23,0.005,200)

plt.figure()
plt.errorbar(voltages_dcA,Betas_vec,yerr=ErrorBetas_vec,fmt='o',capsize=5,markersize=5,color=paleta[1])
plt.plot(xhip,hiperbola(xhip,*popthip))
plt.xlabel('Endcap voltage (V)')
plt.ylabel('Modulation factor')
plt.grid()




#%%
from scipy.special import jv


def expo(x,tau,A,B):
    return A*np.exp(x/tau)+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 <= 5:
        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 InverseMicromotionSpectra_raw(beta, A, det, B):
    ftrap=22.1
    gamma=21
    P = ((jv(0, beta)**2)/((((det)**2)+(0.5*gamma)**2)**2))*(-2*(det))
    i = 1
    #print(P)
    while i <= 3:
        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 A/P+B


"""
Temperatura vs beta con un ajuste exponencial 
"""

popt_exp, pcov_exp = curve_fit(expo,Betas_vec[:10],[t*1e3 for t in Temp_vec[:10]])
popt_quad, pcov_quad = curve_fit(cuadratica,Betas_vec[:10],[t*1e3 for t in Temp_vec[:10]],p0=(1,10))
#popt_rho22, pcov_rho22 = curve_fit(InverseMicromotionSpectra,Betas_vec,[t*1e3 for t in Temp_vec],p0=(10,10,-10,1,20)) #esto ajusta muy bien
#popt_rho22, pcov_rho22 = curve_fit(InverseMicromotionSpectra,Betas_vec, [t*1e3 for t in Temp_vec],p0=(-10,-10,10,1,20)) #esto ajusta muy bien
popt_rho22_raw, pcov_rho22_raw = curve_fit(InverseMicromotionSpectra_raw,Betas_vec[:10], [t*1e3 for t in Temp_vec[:10]],p0=(-10, -10, 1)) #esto ajusta muy bien


print(popt_rho22_raw)

betaslong = np.arange(0,2*2.7,0.01)

print(f'Min temp predicted: {InverseMicromotionSpectra_raw(betaslong,*popt_rho22_raw)[100]}')

plt.figure()
plt.errorbar(Betas_vec[:10],[t*1e3 for t in Temp_vec[:10]],xerr=ErrorBetas_vec[:10], yerr=[t*1e3 for t in ErrorTemp_vec[:10]],fmt='o',capsize=5,markersize=5,color=paleta[3])
#plt.plot(betaslong,expo(betaslong,*popt_exp),label='Ajuste exponencial')
#plt.plot(betaslong,cuadratica(betaslong,*popt_quad),label='Ajuste cuadratico')
#plt.plot(betaslong,InverseMicromotionSpectra(betaslong,*popt_rho22),label='Ajuste cuadratico')
plt.plot(betaslong,InverseMicromotionSpectra_raw(betaslong,*popt_rho22_raw),label='Ajuste cuadratico')

#plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
#plt.axhline(0.538)
plt.xlabel('Modulation factor')
plt.ylabel('Temperature (mK)')
plt.grid()


#%%
"""
Esto no es del super ajuste sino de los ajustes anteriores en donde DetDoppler y offset son puestos a mano
Aca grafico los betas con su error en funcion de la tension variada.
Ademas, hago ajuste lineal para primeros y ultimos puntos, ya que espero que
si la tension hace que la posicion del ion varie linealmente, el beta varia proporcional a dicha posicion.
"""

import seaborn as sns

def lineal(x,a,b):
    return a*x+b

paleta = sns.color_palette("rocket")

betavector = [beta1,beta2,beta3,beta4,beta5,beta6,beta7,beta8,beta9]
errorbetavector = [errorbeta1,errorbeta2,errorbeta3,errorbeta4,errorbeta5,errorbeta6,errorbeta7,errorbeta8,errorbeta9]
voltages_dcA = Voltages[0][1:10]

poptini,pcovini = curve_fit(lineal,voltages_dcA[0:3],betavector[0:3])
poptfin,pcovfin = curve_fit(lineal,voltages_dcA[4:],betavector[4:])

minimum_voltage = -(poptini[1]-poptfin[1])/(poptini[0]-poptfin[0]) #voltaje donde se intersectan las rectas, es decir, donde deberia estar el minimo de micromocion
minimum_modulationfactor = lineal(minimum_voltage,*poptini) #es lo mismo si pongo *poptfin

xini = np.linspace(-0.23,-0.13,100)
xfin = np.linspace(-0.15,0.005,100)

plt.figure()
plt.errorbar(voltages_dcA,betavector,yerr=errorbetavector,fmt='o',capsize=5,markersize=5,color=paleta[1])
plt.plot(xini,lineal(xini,*poptini))
plt.plot(xfin,lineal(xfin,*poptfin))
plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
plt.xlabel('Endcap voltage (V)')
plt.ylabel('Modulation factor')
plt.grid()


#%%
"""
Aca veo la temperatura del ion en funcion del voltaje del endcap, ya que
al cambiar la cantidad de micromocion, cambia la calidad del enfriado
"""
tempvector = np.array([temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8,temp9])*1e3
errortempvector = np.array([errortemp1,errortemp2,errortemp3,errortemp4,errortemp5,errortemp6,errortemp7,errortemp8,errortemp9])*1e3

voltages_dcA = Voltages[0][1:10]

plt.figure()
plt.errorbar(voltages_dcA,tempvector,yerr=errortempvector,fmt='o',capsize=5,markersize=5,color=paleta[3])
plt.axvline(minimum_voltage,linestyle='dashed',color='grey')
plt.xlabel('Endcap voltage (V)')
plt.ylabel('Temperature (mK)')
plt.grid()
plt.ylim(0,2)


#%%
"""
Por las dudas, temperatura en funcion de beta
"""
plt.figure()
plt.errorbar(betavector,tempvector,yerr=errortempvector,xerr=errorbetavector,fmt='o',capsize=5,markersize=5)
plt.xlabel('Modulation factor')
plt.ylabel('Temperature (mK)')
plt.grid()



#%%
"""
Si quiero ver algun parametro del ajuste puntual. el orden es: 0:SG, 1:SP, 2:SCALE1, 3:OFFSET
"""
ki=2
plt.errorbar(np.arange(0,9,1),[popt_1[ki],popt_2[ki],popt_3[ki],popt_4[ki],popt_5[ki],popt_6[ki],popt_7[ki],popt_8[ki],popt_9[ki]],yerr=[np.sqrt(pcov_1[ki,ki]),np.sqrt(pcov_2[ki,ki]),np.sqrt(pcov_3[ki,ki]),np.sqrt(pcov_4[ki,ki]),np.sqrt(pcov_5[ki,ki]),np.sqrt(pcov_6[ki,ki]),np.sqrt(pcov_7[ki,ki]),np.sqrt(pcov_8[ki,ki]),np.sqrt(pcov_9[ki,ki])], fmt='o',capsize=3,markersize=3)


#%%

if False:
    GUARDAR = {}
    # for var in [ kk for kk in globals().keys() if kk.startswith('pop') ]:
    #     print(var)
    #     GUARDAR[var] = globals()[var]
    # print('')
    # for var in [ kk for kk in globals().keys() if kk.startswith('pcov') ]:
    #     print(var)
    #     GUARDAR[var] = globals()[var]

    # print('')
    # for var in [ kk for kk in globals().keys() if kk.startswith('Fitted') ]:
    #     print(var)
    #     GUARDAR[var] = globals()[var]
    # print('')
    for var in [ kk for kk in globals().keys() if kk.endswith('_vec') ]:
        print(var)
        GUARDAR[var] = globals()[var]


    np.savez('analisis_superajuste_PARAMETROS.npz', **GUARDAR )