compu caro

analisis/plots/20250516_CPT_bladetrap/CPT_bladetrap.py

@@ -15,8 +15,8 @@ Primero tengo mediciones de espectros cpt de un ion variando la tension dc_A
 
 #os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20231123_CPTconmicromocion3/Data/')
 #os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20250516_CPT_bladetrap/Data')
-os.chdir(r'C:\Users\Carolina\Documents\artiq_experiments\analisis\plots\20250516_CPT_bladetrap\Data')
-
+#os.chdir(r'C:\Users\Carolina\Documents\artiq_experiments\analisis\plots\20250516_CPT_bladetrap\Data')
+os.chdir("/home/carolina/Documents/artiq_experiments/analisis/plots/20250516_CPT_bladetrap/Data2")
 
 CPT_FILES = """000021278-IR_Scan_simple
 000021360-IR_Scan_simple

analisis/plots/20250516_CPT_bladetrap/CPT_bladetrap2.py

@@ -22,8 +22,8 @@ Primero tengo mediciones de espectros cpt de un ion variando la tension dc_A
 
 #os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20231123_CPTconmicromocion3/Data/')
 #os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20250516_CPT_bladetrap/Data')
-os.chdir(r'C:\Users\Carolina\Documents\artiq_experiments\analisis\plots\20250516_CPT_bladetrap\Data2')
-
+# os.chdir(r'C:\Users\Carolina\Documents\artiq_experiments\analisis\plots\20250516_CPT_bladetrap\Data2')
+os.chdir("/home/carolina/Documents/artiq_experiments/analisis/plots/20250516_CPT_bladetrap/Data2")
 
 VDown = [-177, -175, -173, -171, -170, -170, -168, -166]
 meds = [21368, 21367, 21369, 21370, 21366,21390, 21391, 21392]
@@ -39,12 +39,12 @@ def SeeKeys(files):
 
 Counts = []
 Freqs = []
-
 AmpTisa = []
 UVCPTAmp = []
 No_measures = []
 Voltages = []
 
+
 for i, fname in enumerate(CPT_FILES):
     print(str(i) + ' - ' + fname)
     #print(fname)
@@ -57,8 +57,190 @@ for i, fname in enumerate(CPT_FILES):
     Counts.append(np.array(data['datasets']['counts_spectrum']))
     #AmpTisa.append(np.array(data['datasets']['TISA_CPT_amp']))
     print(np.array(data['datasets']['no_measures']))
+    
+
+#%%
+i = 4
+
+x, y = Freqs[i]*2*1e-6, Counts[i]
+plt.figure(1), plt.clf()
+plt.errorbar(x, y, yerr=np.sqrt(y), capsize=3, fmt=".", label = "Medicion " + str(meds[i]))
+plt.xlabel("Frequency [MHz]")
+plt.ylabel("Counts")
+plt.legend()
+#%%
+
+fig, axes = plt.subplots(4, 2, figsize=(12, 10))  # 4 filas, 2 columnas
+axes = axes.transpose().flatten()  # Para poder iterar fácilmente
+
+for i in range(8):
+    ax = axes[i]
+    x, y = Freqs[i]*2*1e-6, Counts[i]
+    ax.errorbar(x, y, yerr=2*np.sqrt(y), fmt=".", capsize=3)
+    ax.set_title(f'VDown = {VDown[i]} V.  Medición {meds[i]}')
+    ax.set_xlabel('Detuning [MHz]')
+    ax.set_ylabel('Counts')
+    ax.grid(True)
+
+plt.tight_layout()
+plt.show()
+
+
+#%%
+""" Ajustes """
+
+from EITfit.lolo_modelo_full_8niveles import PerformExperiment_8levels_MM
+from scipy.optimize import curve_fit
+import time
+import pandas as pd
+"""
+SUPER AJUSTE (SA)
+"""
+SelectedCurveVec = [i]
+
+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
+
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+alpha = 0
+
+
+drivefreq = 2*np.pi*7.262*1e6
+
+
+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 = []
+
+for selectedcurve in SelectedCurveVec:
+
+    FreqsDR = Freqs[selectedcurve]
+    CountsDR = Counts[selectedcurve]
+
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    CircPr = 1
+    alpha = 0
+
+
+    def FitEIT_MM_single(Freqs, offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP, U, plot=False):
+
+        freqs = [2*f*1e-6-offset for f in Freqs]
+        # BETA1=0
+        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])
+        if plot:
+            return ScaledFluo1, Detunings
+        else:
+            return ScaledFluo1
+        #return ScaledFluo1
+
+    t1 = time.time()
+
+    print(f'Beginning {selectedcurve}')
+    popt_3_SA, pcov_3_SA = curve_fit(FitEIT_MM_single, FreqsDR[1:], CountsDR, p0=[440, -15, 0.3, 10, 1e4, 500, 2, (np.pi**2)*1e-3, 49.5e6], bounds=((0, -50, 0, 0, 0, 0, 0, 0, 48e6), (1000, 0, 2, 20, 5e4, 2e3, 3, (np.pi**2)*10e-3, 50e6)))
+    print(f'Ended {selectedcurve}, time elapsed: {round(time.time()-t1)} s')
+
+    #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, 32e6], bounds=((0, -50, 0, 0, 0, 0, 0, 0, 25e6), (1000, 0, 2, 20, 5e4, 5e4, 10, (np.pi**2)*10e-3, 40e6)))
+
+    popt_SA_vec.append(popt_3_SA)
+    pcov_SA_vec.append(pcov_3_SA)
+
+
+#offset, DetDoppler, SG, SP, SCALE1, OFFSET, BETA1, TEMP, U
 
+    FittedEITpi_3_SA_short, Detunings_3_SA_short = FitEIT_MM_single(FreqsDR, *np.array([popt_3_SA[0],popt_3_SA[1],popt_3_SA[2],popt_3_SA[3],popt_3_SA[4],popt_3_SA[5],popt_3_SA[6],popt_3_SA[7],popt_3_SA[8]]), plot=True)
+    freqslong = np.arange(min(FreqsDR), max(FreqsDR)+FreqsDR[1]-FreqsDR[0], 0.1*(FreqsDR[1]-FreqsDR[0]))
+    
+    # popt_3_SA[4]=1
+    # popt_3_SA[5]=0
+    # popt_3_SA[6]=0
+    FittedEITpi_3_SA_long, Detunings_3_SA_long = FitEIT_MM_single(freqslong, *np.array([popt_3_SA[0],popt_3_SA[1],popt_3_SA[2],popt_3_SA[3],popt_3_SA[4],popt_3_SA[5],popt_3_SA[6],popt_3_SA[7],popt_3_SA[8]]), plot=True)
+    corri = 0
+    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)
 
+    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([c-corri for c in 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)
+    
+
+    # Guardar datos en .csv
+
+    df = pd.DataFrame({
+        'x_data': Detunings_3_SA_short[:],
+        'y_data': CountsDR,
+        'y_error': 2*np.sqrt(CountsDR),
+        })
+    
+    fit_x = [c - corri for c in Detunings_3_SA_long] 
+    fit_df = pd.DataFrame({
+        'fit_x': fit_x,
+        'fit_y': FittedEITpi_3_SA_long
+        })
+    
+    max_len = max(len(df), len(fit_df))
+    df = df.reindex(range(max_len))
+    fit_df = fit_df.reindex(range(max_len))
+
+    # Concatenar los dos DataFrames
+    df_combined = pd.concat([df, fit_df], axis=1)
+    # Añadir los parámetros del fit al final como una nueva fila (puedes usar una fila separadora)
+    fit_param_labels = ['offset', 'DetDoppler', 'SG', 'SP', 'SCALE1', 'OFFSET', 'BETA1', 'TEMP', 'U']
+    fit_params_row = pd.Series([''] * df_combined.shape[1], index=df_combined.columns)
+    param_strs = [f"{label}={val:.5g}" for label, val in zip(fit_param_labels, popt_3_SA)]
+    fit_params_row.iloc[0] = 'Fit parameters:'
+    fit_params_row.iloc[1] = '; '.join(param_strs)
+    
+    # Añadir la fila al final
+    df_combined = pd.concat([df_combined, pd.DataFrame([fit_params_row])], ignore_index=True)
+    
+    # Guardar en CSV
+    df_combined.to_csv(f'DataAndFit_IR_Scan_simple_{meds[i]}.csv', index=False)
 

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/MM_eightLevel_2repumps_AnalysisFunctions.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jul  2 16:30:09 2020
+
+@author: oem
+"""
+
+"""
+ESTE ES EL CODIGO QUE PLOTEA CPT CON MICROMOCION BIEN
+"""
+
+import os
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+from scipy.signal import argrelextrema
+#from EITfit.MM_eightLevel_2repumps_python_scripts import CPTspectrum8levels_MM
+import random
+from scipy.signal import savgol_filter as sf
+
+
+def PerformExperiment_8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None):
+    """
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    #tinicial = time.time()
+    ProbeDetuningVectorL, Fluovector = CPTspectrum8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe, beta, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
+    #tfinal = time.time()
+    
+    #print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
+          
+    return ProbeDetuningVectorL, Fluovector
+
+
+def GenerateNoisyCPT_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, kg, kr, v0, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, kg, kr, v0, drivefreq, freqMin, freqMax, freqStep, circularityprobe, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+def GenerateNoisyCPT_MM_fit(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, beta, drivefreq, freqs, circularityprobe=1, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    
+    Frequencyvector, Fluovector = PerformExperiment_8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, beta, drivefreq, freqs[0], freqs[-1], freqs[1]-freqs[0], circularityprobe, plot=False, solvemode=1, detpvec=None)
+    #NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, Fluovector
+
+
+def SmoothNoisyCPT(Fluo, window=11, poly=3):
+    SmoothenFluo = sf(Fluo, window, poly)
+    return SmoothenFluo
+
+
+
+
+
+

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/MM_eightLevel_2repumps_CPTPlotter.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Sep  1 17:58:39 2020
+
+@author: nico
+"""
+
+
+
+import os
+import numpy as np
+#os.chdir('/home/oem/Nextcloud/G_liaf/liaf-TrampaAnular/Código General/EIT-CPT/Buenos Aires/Experiment Simulations/CPT scripts/Eight Level 2 repumps')
+#from MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels, GenerateNoisyCPT, SmoothNoisyCPT
+import matplotlib.pyplot as plt
+import time
+#from threeLevel_2repumps_AnalysisFunctions import MeasureRelativeFluorescenceFromCPT, IdentifyPolarizationCoincidences, RetrieveAbsoluteCoincidencesBetweenMaps, GetClosestIndex
+import seaborn as sns
+
+#C:\Users\Usuario\Nextcloud\G_liaf\liaf-TrampaAnular\Código General\EIT-CPT\Buenos Aires\Experiment Simulations\CPT scripts\Eight Level 2 repumps
+
+ub = 9.27e-24 #magneton de bohr
+h = 6.63e-34 #cte de planck
+c = (ub/h)*1e-4 #en unidades de MHz/G
+
+u = 2e6 #proportional to the magnetic field of around 5 G
+
+B = (u/(2*np.pi))/c
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+
+lw = 0. #linewidth of the lasers, 0.1 MHz are the actual linewidths of both lasers
+DopplerLaserLinewidth, ProbeLaserLinewidth = lw, lw #ancho de linea de los laseres
+
+
+TempVec = [0e-3] #Temperature vector
+alpha = 0 #angle between lasers, which is zero
+
+#Polarization angles (we can keep it fixed in 90)
+phidoppler, titadoppler = 0, 90
+titaprobe = 90
+phiprobe = 0
+
+#este es el desfasaje exp(i.phi) de la componente de la polarizacion y respecto a la x. Con 1 la polarizacion es lineal
+CircPr = 1  #this has to do with the circularity of the polarizations and since both are linear it is one
+
+#Simulation parameters
+center = -10
+span = 200
+freqMin = center-span*0.5
+freqMax = center+span*0.5
+freqStep = 2e-1
+
+noiseamplitude = 0 #i dont know what it is
+
+#parametros de saturacion de los laseres. g: doppler. p: probe (un rebombeo que scanea), r: repump (otro rebombeo fijo)
+
+"""
+Good case: sg=0.6, sp=9, DetDoppler=-15
+"""
+
+DetDoppler = -25 #nice range: -30 to 0 
+sgvec = [0.6] #nice range: 0.1 to 10 #g is for green but is the doppler
+sp = 8 #nice range: 0.1 to 20 #p is for probe but is the repump
+
+
+
+drivefreq=2*np.pi*22.135*1e6 #ignore it
+#betavec = np.arange(0,1.1,0.1) #ignore it
+betavec=[0] #ignore it
+alphavec = [0] #ignore it
+
+
+
+fig1, ax1 = plt.subplots()
+
+FrequenciesVec = []
+FluorescencesVec = []
+
+for sg in sgvec:
+    for T in TempVec:
+        for alpha in alphavec:
+            for beta in betavec:
+                Frequencies, Fluorescence = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+                
+                FrequenciesVec.append(Frequencies)
+                FluorescencesVec.append(Fluorescence)
+                
+                ax1.plot(Frequencies, [100*f for f in Fluorescence], label=fr'$\alpha={int(alpha*180/np.pi)}°$')  
+                ax1.set_xlabel('Detuning Rebombeo (MHz)')
+                ax1.set_ylabel('Fluorescencia (AU)')
+                ax1.set_title(f'Sdop: {sg}, Spr: {sp}, Temp: {int(T*1e3)} mK')
+                #ax1.legend()
+ax1.grid()
+
+#%%
+import seaborn as sns
+
+paleta=sns.color_palette('mako')
+
+plt.figure()
+plt.plot(Frequencies, [100*f for f in Fluorescence], color=paleta[1], linewidth=3)  
+plt.grid()
+plt.axvline(-25,color=paleta[2], linestyle='dashed')
+plt.xlabel(r'$\Delta_2$ (MHz)', fontsize=25, fontname='STIXgeneral')
+plt.ylabel('Fluorescence', fontsize=18, fontname='STIXgeneral')
+#%%
+#Este bloque ajusta a las curvas con un beta de micromocion de 0
+
+from scipy.optimize import curve_fit
+
+def FitEIT_MM(freqs, Temp):
+    BETA = 0
+    scale=1
+    offset=0
+    Detunings, Fluorescence = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA, drivefreq, freqMin, freqMax, freqStep, circularityprobe=CircPr, plot=False, solvemode=1, detpvec=None)
+    ScaledFluo = [f*scale + offset for f in Fluorescence]
+    return ScaledFluo
+
+TempMedidas = []
+FittedEIT_fluosVec = []
+
+for j in range(len(betavec)):
+    SelectedFluo = FluorescencesVec[j]
+    SelectedFreqs = FrequenciesVec[j]
+    
+    popt_mm, pcov_mm = curve_fit(FitEIT_MM, SelectedFreqs, SelectedFluo, p0=[1e-3], bounds=((0), (10e-3)))
+    
+    TempMedidas.append(1e3*popt_mm[2])
+    
+    print(popt_mm)
+    
+    FittedEIT_fluo = FitEIT_MM(SelectedFreqs, *popt_mm)
+    
+    FittedEIT_fluosVec.append(FittedEIT_fluo)
+    
+    plt.figure()
+    plt.plot(SelectedFreqs, SelectedFluo, 'o')
+    plt.plot(SelectedFreqs, FittedEIT_fluo)
+
+plt.figure()
+for i in range(len(FluorescencesVec)):
+    plt.plot(SelectedFreqs, FluorescencesVec[i], 'o', markersize=3)
+    plt.plot(SelectedFreqs, FittedEIT_fluosVec[i])    
+
+
+plt.figure()
+plt.plot(betavec, TempMedidas, 'o', markersize=10)
+plt.xlabel('Beta')
+plt.ylabel('Temperatura medida (mK)')
+plt.axhline(T*1e3, label='Temperatura real', linestyle='--', color='red')
+plt.legend()
+plt.grid()
\ No newline at end of file

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/MM_eightLevel_2repumps_python_scripts.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr  7 22:30:01 2020
+
+@author: nico
+"""
+
+"""
+ESTE ES EL CODIGO QUE PLOTEA CPT CON MICROMOCION BIEN
+"""
+
+#ESTE CODIGO ES EL PRINCIPAL PARA PLOTEAR CPT TEORICOS
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+from scipy.signal import argrelextrema
+
+
+"""
+Scripts para el calculo de la curva CPT
+"""
+
+
+def H0matrix(Detg, Detp, u): 
+    """
+    Calcula la matriz H0 en donde dr es el detuning del doppler, dp es el retuning del repump y u es el campo magnético en Hz/Gauss.
+    Para esto se toma la energía del nivel P como 0
+    """
+    
+    eigenEnergies = (Detg-u, Detg+u, -u/3, u/3, Detp-6*u/5, Detp-2*u/5, Detp+2*u/5, Detp+6*u/5) #pagina 26 de Oberst. los lande del calcio son iguales a Bario.
+    H0 = np.diag(eigenEnergies)
+    
+    return H0
+
+
+def HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe=1):
+    """
+    Calcula la matriz de interacción Hsp + Hpd, en donde rabR es la frecuencia de rabi de la transición Doppler SP, 
+    rabP es la frecuencia de rabi de la transición repump DP, y las componentes ei_r y ei_p son las componentes de la polarización
+    del campo eléctrico incidente de doppler y repump respectivamente. Deben estar normalizadas a 1
+    """
+    HI = np.zeros((8, 8), dtype=np.complex_)
+    
+    
+    i, j = 1, 3
+    HI[i-1, j-1] = (rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1]) 
+    
+    i, j = 1, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 2, 3
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(-1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+     
+    i, j = 2, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])   
+    
+    i, j = 3, 5
+    HI[i-1, j-1] = -(rabP/2) * np.sin(titaprobe)*(np.cos(phiprobe)-1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 3, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 3, 7
+    HI[i-1, j-1] = rabP/np.sqrt(12) * np.sin(titaprobe)*(np.cos(phiprobe)+1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    
+    i, j = 4, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(12)) * np.sin(titaprobe)*(np.cos(phiprobe)-1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 4, 7
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe) 
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 4, 8
+    HI[i-1, j-1] = (rabP/2) * np.sin(titaprobe)*(np.cos(phiprobe)+1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    
+    return HI
+
+
+def LtempCalculus(beta, drivefreq, forma=1):
+    Hint = np.zeros((8, 8), dtype=np.complex_)
+    
+    ampg=beta*drivefreq
+    ampr=beta*drivefreq*(397/866)        
+    #ampr=beta*drivefreq
+    Hint[0,0] = ampg
+    Hint[1,1] = ampg
+    Hint[4,4] = ampr
+    Hint[5,5] = ampr
+    Hint[6,6] = ampr
+    Hint[7,7] = ampr
+    
+    if forma==1:
+        Ltemp = np.zeros((64, 64), dtype=np.complex_)
+        """
+        for r in range(8):
+            for q in range(8):
+                for k in range(8):
+                    for j in range(8):
+                        Ltemp[r*8+q][k*8+j] = (-1j)*(Hint[r,k]*int(j==q) - Hint[j,q]*int(r==k))
+        """
+        """
+        for r in range(8):
+            for q in range(8):
+                for k in range(8):
+                    for j in range(8):
+                        if r==k and j==q:
+                            Ltemp[r*8+q][k*8+j] = (-1j)*(Hint[r,k] - Hint[j,q])
+        """
+
+        
+        for r in range(8):
+            for q in range(8):
+                if r!=q:
+                    Ltemp[r*8+q][r*8+q] = (-1j)*(Hint[r,r] - Hint[q,q])
+        
+        
+    if forma==2:
+        deltaKro = np.diag([1, 1, 1, 1, 1, 1, 1, 1])
+        Ltemp = (-1j)*(np.kron(Hint, deltaKro) - np.kron(deltaKro, Hint))
+
+    Omega = np.zeros((64, 64), dtype=np.complex_)
+    
+    for i in range(64):
+        Omega[i, i] = (1j)*drivefreq
+    
+
+    return np.matrix(Ltemp), np.matrix(Omega)
+
+
+def GetL1(Ltemp, L0, Omega, nmax):
+    """
+    Devuelve Splus0 y Sminus0
+    """
+    
+    Sp = (-1)*(np.matrix(np.linalg.inv(L0 - (nmax+1)*Omega))*0.5*np.matrix(Ltemp))
+    Sm = (-1)*(np.matrix(np.linalg.inv(L0 + (nmax+1)*Omega))*0.5*np.matrix(Ltemp))
+    
+    for n in list(range(nmax+1))[(nmax+1)::-1][0:len(list(range(nmax+1))[(nmax+1)::-1])-1]: #jaja esto solo es para que vaya de nmax a 1 bajando. debe haber algo mas facil pero kcio
+        Sp = (-1)*(np.matrix(np.linalg.inv(L0 - n*Omega + (0.5*Ltemp*np.matrix(Sp))))*0.5*np.matrix(Ltemp))
+        Sm = (-1)*(np.matrix(np.linalg.inv(L0 + n*Omega + (0.5*Ltemp*np.matrix(Sm))))*0.5*np.matrix(Ltemp))
+    L1 = 0.5*np.matrix(Ltemp)*(np.matrix(Sp) + np.matrix(Sm))
+    
+    return L1
+
+
+def EffectiveL(gPS, gPD, lwg, lwp):
+
+    """
+    Siendo Heff = H + EffectiveL, calcula dicho EffectiveL que es (-0.5j)*sumatoria(CmDaga*Cm) que luego sirve para calcular el Liouvilliano
+    """
+    Leff = np.zeros((8, 8), dtype=np.complex_)
+    
+    Leff[0, 0] = 2*lwg
+    Leff[1, 1] = 2*lwg
+    Leff[2, 2] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[3, 3] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[4, 4] = 2*lwp
+    Leff[5, 5] = 2*lwp
+    Leff[6, 6] = 2*lwp
+    Leff[7, 7] = 2*lwp
+    
+    return (-0.5j)*Leff
+
+
+def CalculateSingleMmatrix(gPS, gPD, lwg, lwp):
+    """
+    Si tomamos el Liuvilliano como L = (-j)*(Heff*deltak - Heffdaga*deltak) + sum(Mm),
+    esta funcion calcula dichos Mm, que tienen dimensión 64x64 ya que esa es la dimensión del L. Estas componentes
+    salen de hacer la cuenta a mano conociendo los Cm y considerando que Mm[8*(r-1)+s, 8*(k-1)+j] = Cm[r,l] + Cmdaga[j,s] = Cm[r,l] + Cm[s,j]
+    ya que los componentes de Cm son reales.
+    Esta M es la suma de las 8 matrices M.
+    
+    """
+    
+    M = np.matrix(np.zeros((64, 64), dtype=np.complex_))
+    
+    M[0,27] = (2/3)*gPS
+    
+    M[9,18] = (2/3)*gPS
+    
+    M[0,18] = (1/3)*gPS
+    M[1,19] = -(1/3)*gPS
+    M[8,26] = -(1/3)*gPS
+    M[9,27] = (1/3)*gPS
+    
+    M[36,18] = (1/2)*gPD
+    M[37,19] = (1/np.sqrt(12))*gPD
+    M[44,26] = (1/np.sqrt(12))*gPD
+    M[45,27] = (1/6)*gPD
+    
+    M[54,18] = (1/6)*gPD
+    M[55,19] = (1/np.sqrt(12))*gPD
+    M[62,26] = (1/np.sqrt(12))*gPD
+    M[63,27] = (1/2)*gPD
+    
+    M[45,18] = (1/3)*gPD
+    M[46,19] = (1/3)*gPD
+    M[53,26] = (1/3)*gPD
+    M[54,27] = (1/3)*gPD
+    
+    M[0,0] = 2*lwg
+    M[1,1] = 2*lwg
+    M[8,8] = 2*lwg
+    M[9,9] = 2*lwg
+
+    
+    #M[36, 45] = lwp
+    
+    
+    for k in [36, 37, 38, 39, 44, 45, 46, 47, 52, 53, 54, 55, 60, 61, 62, 63]:
+        M[k,k]=2*lwp
+     
+    return M
+
+def dopplerBroadening(wlg, wlp, alpha, T, mcalcio = 6.655e-23*1e-3):
+    """
+    Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
+    wlg es la longitud de onda doppler, wlp la longitud de onda repump, T la temperatura del ion en kelvin, y alpha (en rads) el ángulo
+    que forman ambos láseres.
+    """
+    
+    kboltzmann = 1.38e-23 #J/K
+    
+    gammaD = (2*np.pi)*np.sqrt((1/(wlg*wlg)) + (1/(wlp*wlp)) - 2*(1/(wlg*wlp))*np.cos(alpha))*np.sqrt(kboltzmann*T/(2*mcalcio))
+
+    return gammaD
+    
+
+def FullL_MM(rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp = 0, 
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, beta=0, drivefreq=2*np.pi*22.135*1e6, T = 0, alpha = 0, circularityprobe=1):
+
+    """
+    Calcula el Liouvilliano total de manera explícita índice a índice. Suma aparte las componentes de las matrices M.
+    Es la más eficiente hasta ahora.
+    """
+    
+    
+    db = dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)
+    
+    lwg = np.sqrt(lwg**2 + db**2)
+    lwp = np.sqrt(lwp**2 + db**2)    
+    
+    CC = EffectiveL(gPS, gPD, lwg, lwp)
+    
+    Heff = H0matrix(Detg, Detp, u) + HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe) + CC
+    
+    Heffdaga = np.matrix(Heff).getH()
+    
+    Lfullpartial = np.zeros((64, 64), dtype=np.complex_)
+    
+    for r  in range(8):
+        for q in range(8):
+            for k in range(8):
+                for j in range(8):
+                    if j!=q and r!=k:
+                        pass
+                    elif j==q and r!=k:
+                        if (r < 2 and k > 3) or (k < 2 and r > 3) or (r > 3 and k > 3) or (r==0 and k==1) or (r==1 and k==0) or (r==2 and k==3) or (r==3 and k==2): #todo esto sale de analizar explicitamente la matriz y tratar de no calcular cosas de más que dan cero
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k])
+                    elif j!=q and r==k:
+                        if (j < 2 and q > 3) or (q < 2 and j > 3) or (j > 3 and q > 3) or (j==0 and q==1) or (j==1 and q==0) or (j==2 and q==3) or (j==3 and q==2):
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(-Heffdaga[j,q])
+                    else:
+                        if Heff[r,k] == Heffdaga[j,q]:
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k]-Heffdaga[j,q])
+
+    M = CalculateSingleMmatrix(gPS, gPD, lwg, lwp) 
+    L0 = np.array(np.matrix(Lfullpartial) + M)
+    
+    #ESTA PARTE ES CUANDO AGREGAS MICROMOCION
+    nmax = 3
+    
+    #print(nmax)
+    Ltemp, Omega = LtempCalculus(beta, drivefreq)
+    #print(factor)
+    L1 = GetL1(Ltemp, L0, Omega, nmax)
+    Lfull = L0 + L1 #ESA CORRECCION ESTA EN L1
+    #HASTA ACA
+    
+    #NORMALIZACION DE RHO
+    i = 0 
+    while i < 64:
+        if i%9 == 0:
+            Lfull[0, i] = 1
+        else:
+            Lfull[0, i] = 0
+        i = i + 1    
+    
+    return Lfull 
+
+
+"""
+Scripts para correr un experimento y hacer el análisis de los datos
+"""
+
+
+def CPTspectrum8levels_MM(sg, sp, gPS, gPD, Detg, u, lwg, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, Circularityprobe, beta, drivefreq, freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+
+    """
+    ESTA ES LA FUNCION QUE ESTAMOS USANDO
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+    
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180) 
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6+0*freqStep*1e6, freqStep*1e6)
+    Detg = 2*np.pi*Detg*1e6
+    #lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+    lwg, lwp = lwg*1e6, lwp*1e6
+
+    rabG = sg*gPS
+    rabP = sp*gPD
+    #u = 2*np.pi*u*1e6
+    
+    Fluovector = []
+    
+    tinicial = time.time()
+    
+    for Detp in DetProbeVector:
+        L = FullL_MM(rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, Temp, alpha, Circularityprobe)
+        if solvemode == 1:
+            rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)]))
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27]))
+            Fluovector.append(Fluo)                
+
+        if solvemode == 2:
+            Linv = np.linalg.inv(L) 
+            rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)
+        
+    tfinal = time.time()
+        
+    print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
+    
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+    
+    if plot:
+        plt.xlabel('Probe detuning (MHz)')
+        plt.ylabel('Fluorescence (A.U.)')
+        plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+        plt.legend()
+        
+    return DetProbeVectorMHz, Fluovector
+
+#%%
+if __name__ == "__main__":
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+    
+    B = 25 #campo magnetico en gauss
+    
+    u = c*B
+    
+    
+    sg, sr, sp = 0.5, 1.5, 4 #parámetros de saturación del doppler y repump
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+    rabG, rabR, rabP = sg*gPS, sr*gPD, sp*gPD #frecuencias de rabi
+    
+    lwg, lwr, lwp = 0.3, 0.3, 0.3 #ancho de linea de los laseres
+    
+    Detg = -25
+    Detr = 20 #detuning del doppler y repump
+    
+    Temp = 0.0e-3 #temperatura en K
+    alpha = 0*(np.pi/180) #angulo entre los láseres
+    
+    phidoppler, titadoppler = 0, 90
+    phirepump, titarepump = 0, 90
+    
+    phiprobe, titaprobe = 0, 90
+    
+    
+    plotCPT = False
+    
+    freqMin = -50
+    freqMax = 50
+    freqStep = 5e-2
+    
+    Frequencyvector, Fluovector = CPTspectrum8levels_MM(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=plotCPT, solvemode=1)
+    
+    plt.plot(Frequencyvector, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    plt.xlabel('Probe detuning (MHz)')
+    plt.ylabel('Fluorescence (A.U.)')
+    
+

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/lolo_modelo_full_3niveles.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+4 de ene 2024
+
+@author: lolo
+
+reingenieria del código que anda
+
+Surge de fusionar
+Data/EITfit/threeLevel_2repumps_AnalysisFunctions.py
+Data/EITfit/threeLevel_2repumps_linealpol_python_scripts.py
+
+MAPA de FUNCIONES
+
+
+GenerateNoisyCPT_fit
+ |--> PerformExperiment_8levels_fixedRabi
+   |--> CPTspectrum8levels_fixedRabi --> ndarray,ndarray
+     |--> FullL_efficient --> ndarray
+       |--> dopplerBroadening --> float
+       |--> EffectiveL --> ndarray
+       |--> H0matrix --> ndarray
+       |--> HImatrix --> ndarray
+       |--> CalculateSingleMmatrix --> ndarray
+       |--> Lplusminus --> ndarray,ndarray,ndarray
+       |--> GetL1 --> ndarray
+
+"""
+
+# pylint: disable=C0301,R0913,R0914,W0621
+
+
+
+import os
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+from scipy.signal import argrelextrema
+# from EITfit.threeLevel_2repumps_linealpol_python_scripts import CPTspectrum8levels, CPTspectrum8levels_fixedRabi
+import random
+from scipy.signal import savgol_filter as sf
+
+from  numba import jit,njit
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+@njit
+def make_diag(vec):
+    "Construye matris diagonal desde una lista o vector"
+    return np.eye(len(vec))*np.array(vec).reshape(-1,1)
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+"""
+Scripts para el calculo de la curva CPT
+"""
+
+@njit
+def H0matrix(Detg, Detp, u):
+    """
+    Calcula la matriz H0 en donde dr es el detuning del doppler, dp es el retuning
+    del repump y u es el campo magnético en Hz/Gauss.
+    Para esto se toma la energía del nivel P como 0
+    """
+
+    eigenEnergies = (Detg-u, Detg+u, -u/3, u/3, Detp-6*u/5,
+                    Detp-2*u/5, Detp+2*u/5, Detp+6*u/5)
+    #pagina 26 de Oberst. los lande del calcio son iguales a Bario.
+    # H0 = np.diag(eigenEnergies)
+    # H0 = np.eye(len(eigenEnergies))*np.array(eigenEnergies).reshape(-1,1)
+    H0 = make_diag(eigenEnergies)
+    return H0
+
+@njit
+def HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe):
+    """
+    Calcula la matriz de interacción Hsp + Hpd, en donde rabR es la frecuencia de rabi de la transición Doppler SP,
+    rabP es la frecuencia de rabi de la transición repump DP, y las componentes ei_r y ei_p son las componentes de la polarización
+    del campo eléctrico incidente de doppler y repump respectivamente. Deben estar normalizadas a 1
+    """
+    HI = np.zeros((8, 8), dtype=np.complex_)
+
+    i, j = 1, 3
+    HI[i-1, j-1] = (rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 1, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 2, 3
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(-1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 2, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+
+    i, j = 3, 5
+    HI[i-1, j-1] = -(rabP/2) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 3, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 3, 7
+    HI[i-1, j-1] = rabP/np.sqrt(12) * np.sin(titaprobe)*np.exp(1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+
+    i, j = 4, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(12)) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 4, 7
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 4, 8
+    HI[i-1, j-1] = (rabP/2) * np.sin(titaprobe)*np.exp(1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    return HI
+
+
+@njit
+def Lplusminus(detr, detp, phirepump, titarepump, forma=1):
+    Hintplus = np.zeros((8, 8), dtype=np.complex_)
+    Hintminus = np.zeros((8, 8), dtype=np.complex_)
+
+
+    Hintplus[4, 2] = (-1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
+    Hintplus[5, 2] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintplus[6, 2] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
+
+    Hintplus[5, 3] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
+    Hintplus[6, 3] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintplus[7, 3] = (1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
+
+
+    Hintminus[2, 4] = (-1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
+    Hintminus[2, 5] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintminus[2, 6] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
+
+    Hintminus[3, 5] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
+    Hintminus[3, 6] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintminus[3, 7] = (1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
+
+    # if forma==1:
+    if True:
+        Lplus = np.zeros((64, 64), dtype=np.complex_)
+        Lminus = np.zeros((64, 64), dtype=np.complex_)
+        DeltaBar = np.zeros((64, 64), dtype=np.complex_)
+
+        for r  in range(8):
+            for q in range(8):
+                for k in range(8):
+                    for j in range(8):
+                        if j==q:
+                            if (k==2 or k==3) and r > 3:
+                                Lplus[r*8+q][k*8+j] = (-1j)*(Hintplus[r,k])
+                            if (r==2 or r==3) and k > 3:
+                                Lminus[r*8+q][k*8+j] = (-1j)*(Hintminus[r,k])
+                        elif r==k:
+                            if (q==2 or q==3) and j > 3:
+                                Lplus[r*8+q][k*8+j] = (-1j)*(- Hintplus[j,q])
+                            if (j==2 or j==3) and q > 3:
+                                Lminus[r*8+q][k*8+j] = (-1j)*(- Hintminus[j,q])
+    # if forma==2:
+    #     deltaKro = np.diag([1, 1, 1, 1, 1, 1, 1, 1])
+    #     Lplus = (-1j)*(np.kron(Hintplus, deltaKro) - np.kron(deltaKro, Hintplus))
+    #     Lminus = (-1j)*(np.kron(Hintminus, deltaKro) - np.kron(deltaKro, Hintminus))
+
+    DeltaBar = np.zeros((64, 64), dtype=np.complex_)
+    for i in range(64):
+
+        DeltaBar[i, i] = (1j)*(detr - detp)
+
+    # return np.matrix(Lminus), np.matrix(Lplus), np.matrix(DeltaBar)
+    return Lminus, Lplus, DeltaBar
+
+@njit
+def GetL1(Lplus, Lminus, DeltaBar, L0, rabR, nmax):
+    """
+    Devuelve Splus0 y Sminus0
+    """
+
+    # Sp = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 - (nmax+1)*DeltaBar))*np.matrix(Lplus))
+    # Sm = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 + (nmax+1)*DeltaBar))*np.matrix(Lminus))
+
+    Sp = (-1)*(0.5*rabR)*(np.linalg.inv(L0 - (nmax+1)*DeltaBar)).dot(Lplus)
+    Sm = (-1)*(0.5*rabR)*(np.linalg.inv(L0 + (nmax+1)*DeltaBar)).dot(Lminus)
+
+    for n in list(range(nmax+1))[(nmax+1)::-1][0:len(list(range(nmax+1))[(nmax+1)::-1])-1]: #jaja esto solo es para que vaya de nmax a 1 bajando. debe haber algo mas facil pero kcio
+        # Sp = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 - n*DeltaBar + rabR*(Lminus*np.matrix(Sp))))*np.matrix(Lplus))
+        # Sm = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 + n*DeltaBar + rabR*(Lplus*np.matrix(Sm))))*np.matrix(Lminus))
+
+        Sp = (-1)*(rabR)*((np.linalg.inv(L0 - n*DeltaBar + rabR*(Lminus.dot(Sp)))).dot(Lplus))
+        Sm = (-1)*(rabR)*((np.linalg.inv(L0 + n*DeltaBar + rabR*(Lplus.dot(Sm)))).dot(Lminus))
+
+    L1 = 0.5*rabR*(Lminus.dot(Sp) + Lplus.dot(Sm))
+
+    return L1
+
+@njit
+def EffectiveL(gPS, gPD, lwg, lwr, lwp):
+
+    """
+    Siendo Heff = H + EffectiveL, calcula dicho EffectiveL que es (-0.5j)*sumatoria(CmDaga*Cm) que luego sirve para calcular el Liouvilliano
+    """
+    Leff = np.zeros((8, 8), dtype=np.complex_)
+
+    Leff[0, 0] = 2*lwg
+    Leff[1, 1] = 2*lwg
+    Leff[2, 2] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[3, 3] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[4, 4] = 2*(lwr + lwp)
+    Leff[5, 5] = 2*(lwr + lwp)
+    Leff[6, 6] = 2*(lwr + lwp)
+    Leff[7, 7] = 2*(lwr + lwp)
+
+    return (-0.5j)*Leff
+
+@njit
+def CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp):
+    """
+    Si tomamos el Liuvilliano como L = (-j)*(Heff*deltak - Heffdaga*deltak) + sum(Mm),
+    esta funcion calcula dichos Mm, que tienen dimensión 64x64 ya que esa es la dimensión del L. Estas componentes
+    salen de hacer la cuenta a mano conociendo los Cm y considerando que Mm[8*(r-1)+s, 8*(k-1)+j] = Cm[r,l] + Cmdaga[j,s] = Cm[r,l] + Cm[s,j]
+    ya que los componentes de Cm son reales.
+    Esta M es la suma de las 8 matrices M.
+
+    """
+
+    M =np.zeros((64, 64), dtype=np.complex_)
+
+    M[0,27] = (2/3)*gPS
+
+    M[9,18] = (2/3)*gPS
+
+    M[0,18] = (1/3)*gPS
+    M[1,19] = -(1/3)*gPS
+    M[8,26] = -(1/3)*gPS
+    M[9,27] = (1/3)*gPS
+
+    M[36,18] = (1/2)*gPD
+    M[37,19] = (1/np.sqrt(12))*gPD
+    M[44,26] = (1/np.sqrt(12))*gPD
+    M[45,27] = (1/6)*gPD
+
+    M[54,18] = (1/6)*gPD
+    M[55,19] = (1/np.sqrt(12))*gPD
+    M[62,26] = (1/np.sqrt(12))*gPD
+    M[63,27] = (1/2)*gPD
+
+    M[45,18] = (1/3)*gPD
+    M[46,19] = (1/3)*gPD
+    M[53,26] = (1/3)*gPD
+    M[54,27] = (1/3)*gPD
+
+    M[0,0] = 2*lwg
+    M[1,1] = 2*lwg
+    M[8,8] = 2*lwg
+    M[9,9] = 2*lwg
+
+    factor1 = 1
+    factor2 = 1
+    factor3 = 1
+    factor4 = 1
+
+
+    #M[36, 45] = lwp
+
+
+    M[36,36] = 2*(lwr + factor1*lwp)
+    M[37,37] = 2*(lwr + factor1*lwp)
+    M[38,38] = 2*(lwr + factor1*lwp)
+    M[39,39] = 2*(lwr + factor1*lwp)
+    M[44,44] = 2*(lwr + factor2*lwp)
+    M[45,45] = 2*(lwr + factor2*lwp)
+    M[46,46] = 2*(lwr + factor2*lwp)
+    M[47,47] = 2*(lwr + factor2*lwp)
+    M[52,52] = 2*(lwr + factor3*lwp)
+    M[53,53] = 2*(lwr + factor3*lwp)
+    M[54,54] = 2*(lwr + factor3*lwp)
+    M[55,55] = 2*(lwr + factor3*lwp)
+    M[60,60] = 2*(lwr + factor4*lwp)
+    M[61,61] = 2*(lwr + factor4*lwp)
+    M[62,62] = 2*(lwr + factor4*lwp)
+    M[63,63] = 2*(lwr + factor4*lwp)
+
+    return M
+
+
+@njit
+def dopplerBroadening(wlg, wlp, alpha, T, mcalcio = 6.655e-23*1e-3):
+    """
+    Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
+    wlg es la longitud de onda doppler, wlp la longitud de onda repump, T la temperatura del ion en kelvin, y alpha (en rads) el ángulo
+    que forman ambos láseres.
+    """
+
+    kboltzmann = 1.38e-23 #J/K
+
+    gammaD = (2*np.pi)*np.sqrt((1/(wlg*wlg)) + (1/(wlp*wlp)) - 2*(1/(wlg*wlp))*np.cos(alpha))*np.sqrt(kboltzmann*T/(2*mcalcio))
+
+    return gammaD
+
+@njit
+def FullL_efficient(rabG, rabR, rabP, gPS = 0, gPD = 0, Detg = 0, Detr = 0, Detp = 0, u = 0, lwg = 0, lwr=0, lwp = 0,
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, phirepump=0, titarepump=0, T = 0, alpha = 0):
+
+    """
+    Calcula el Liouvilliano total de manera explícita índice a índice. Suma aparte las componentes de las matrices M.
+    Es la más eficiente hasta ahora.
+    """
+
+
+    db = dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)
+
+    #lwr = np.sqrt(lwr**2 + dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)**2)
+    lwg = np.sqrt(lwg**2 + db**2)
+    lwr = np.sqrt(lwr**2 + db**2)
+
+    CC = EffectiveL(gPS, gPD, lwg, lwr, lwp)
+
+    Heff = H0matrix(Detg, Detp, u) + HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe) + CC
+
+    # Heffdaga = np.matrix(Heff).getH()
+    Heffdaga = Heff.transpose().conj()
+
+    Lfullpartial = np.zeros((64, 64), dtype=np.complex_)
+
+    for r  in range(8):
+        for q in range(8):
+            for k in range(8):
+                for j in range(8):
+                    if j!=q and r!=k:
+                        pass
+                    elif j==q and r!=k:
+                        if (r < 2 and k > 3) or (k < 2 and r > 3) or (r > 3 and k > 3) or (r==0 and k==1) or (r==1 and k==0) or (r==2 and k==3) or (r==3 and k==2): #todo esto sale de analizar explicitamente la matriz y tratar de no calcular cosas de más que dan cero
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k])
+                    elif j!=q and r==k:
+                        if (j < 2 and q > 3) or (q < 2 and j > 3) or (j > 3 and q > 3) or (j==0 and q==1) or (j==1 and q==0) or (j==2 and q==3) or (j==3 and q==2):
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(-Heffdaga[j,q])
+                    else:
+                        if Heff[r,k] == Heffdaga[j,q]:
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k]-Heffdaga[j,q])
+
+    M = CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp)
+    L0 = Lfullpartial + M
+
+    nmax = 1
+    Lminus, Lplus, DeltaBar = Lplusminus(Detr, Detp, phirepump, titarepump)
+    factor1 = np.exp(1j*0.2*np.pi)
+    factor2 = np.exp(-1j*0.2*np.pi)
+    #print(factor)
+    L1 = GetL1(factor1*Lplus, factor2*Lminus, DeltaBar, L0, rabR, nmax)
+    Lfull = L0 + L1
+
+
+    #NORMALIZACION DE RHO
+    i = 0
+    while i < 64:
+        if i%9 == 0:
+            Lfull[0, i] = 1
+        else:
+            Lfull[0, i] = 0
+        i = i + 1
+
+    return Lfull
+
+
+"""
+Scripts para correr un experimento y hacer el análisis de los datos
+"""
+
+def CalculoTeoricoDarkResonances(u, titadoppler):
+    if titadoppler==0:
+        NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
+    elif titadoppler==90:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+
+    PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
+
+    return NegativeDR, PositiveDR
+
+
+def CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump,
+                      freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
+
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
+    Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
+    lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+
+    #u = 2*np.pi*u*1e6
+
+    Fluovector = []
+
+    tinicial = time.time()
+
+    for Detp in DetProbeVector:
+        L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
+
+        # if solvemode == 1:
+        if True:
+            rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)],dtype=np.complex_))
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)
+
+        # if solvemode == 2:
+        #     Linv = np.linalg.inv(L)
+        #     rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+        #     Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+        #     Fluovector.append(Fluo)
+
+    tfinal = time.time()
+
+    print('Done, Total time: ', round((tfinal-tinicial), 2), "s")
+
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+
+    if plot:
+        plt.xlabel('Probe detuning (MHz)')
+        plt.ylabel('Fluorescence (A.U.)')
+        plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+        plt.legend()
+
+    return DetProbeVectorMHz, Fluovector
+
+@njit
+def CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump,
+                      freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+
+
+
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
+
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
+    Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
+    #lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+    lwg, lwr, lwp = lwg*1e6, lwr*1e6, lwp*1e6
+
+    rabG = sg*gPS
+    rabR = sr*gPD
+    rabP = sp*gPD
+    #u = 2*np.pi*u*1e6
+
+    Fluovector = []
+
+    # tinicial = time.time()
+
+    for Detp in DetProbeVector:
+
+        L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
+
+        # if solvemode == 1:
+        if True:
+            coh = 5
+            rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)],dtype=np.complex_))
+            #Fluo = np.abs(rhovectorized[coh])
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)
+
+        # if solvemode == 2:
+        #     Linv = np.linalg.inv(L)
+        #     rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+        #     Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+        #     Fluovector.append(Fluo)
+
+    # tfinal = time.time()
+
+    # print('Done, Total time: ', round((tfinal-tinicial), 2), "s")
+
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+
+    # if plot:
+    #     plt.xlabel('Probe detuning (MHz)')
+    #     plt.ylabel('Fluorescence (A.U.)')
+    #     plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    #     plt.legend()
+
+    return DetProbeVectorMHz, Fluovector
+#%%
+
+if __name__ == "__main__":
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+
+    B = 25 #campo magnetico en gauss
+
+    u = c*B
+
+
+    sg, sr, sp = 0.5, 1.5, 4 #parámetros de saturación del doppler y repump
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+    rabG, rabR, rabP = sg*gPS, sr*gPD, sp*gPD #frecuencias de rabi
+
+    lwg, lwr, lwp = 0.3, 0.3, 0.3 #ancho de linea de los laseres
+
+    Detg = -25
+    Detr = 20 #detuning del doppler y repump
+
+    Temp = 0.0e-3 #temperatura en K
+    alpha = 0*(np.pi/180) #angulo entre los láseres
+
+    phidoppler, titadoppler = 0, 90
+    phirepump, titarepump = 0, 90
+
+    phiprobe, titaprobe = 0, 90
+
+
+    plotCPT = False
+
+    freqMin = -50
+    freqMax = 50
+    freqStep = 5e-2
+
+    Frequencyvector, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=plotCPT, solvemode=1)
+
+    NegativeDR, PositiveDR = CalculoTeoricoDarkResonances(u/(2*np.pi*1e6), titadoppler)
+
+    plt.plot(Frequencyvector, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    plt.xlabel('Probe detuning (MHz)')
+    plt.ylabel('Fluorescence (A.U.)')
+
+    for PDR in PositiveDR:
+        plt.axvline(Detr+PDR, linestyle='--', linewidth=0.5, color='red')
+    for NDR in NegativeDR:
+        plt.axvline(Detg+NDR, linestyle='--', linewidth=0.5, color='blue')
+
+
+
+
+    #parametros que andan piola:
+
+
+    """
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+
+    B = 17 #campo magnetico en gauss
+
+    u = c*B
+    #u = 80e6
+
+    sr, sp = 0.53, 4.2
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+    rabR, rabP = sr*gPS, sp*gPD
+    lw = 2*np.pi * 0.33e6
+    lwr, lwp = lw, lw #ancho de linea de los laseres
+    dr_spec = - 2*np.pi* 26e6
+
+
+    freqSteps = 500
+    freqMin = -100e6
+    freqMax = 100e6
+
+    dps = 2*np.pi*np.linspace(freqMin, freqMax, freqSteps)
+
+    #dps = [-30e6]
+
+    alfar = 90*(np.pi/180)
+    ex_r, ey_r, ez_r = np.sin(alfar)*np.cos(0), np.sin(alfar)*np.sin(0), np.cos(alfar)
+
+
+    alfap = 90*(np.pi/180)
+    ex_p, ey_p, ez_p = np.sin(alfap)*np.cos(0), np.sin(alfap)*np.sin(0), np.cos(alfap)
+
+    """
+
+
+
+
+
+
+
+
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+
+
+
+def CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump):
+
+    if titadoppler==0:
+        NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
+    elif titadoppler==90:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+    else:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+
+    PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
+
+    return [detuningdoppler + dr for dr in NegativeDR], [detuningrepump + dr for dr in PositiveDR]
+
+
+def GetClosestIndex(Vector, value, tolerance=1e-3):
+
+    i = 0
+    while i<len(Vector):
+        if abs(Vector[i] - value) < tolerance:
+            return i
+        else:
+            i = i + 1
+    return GetClosestIndex(Vector, value, tolerance=2*tolerance)
+
+
+def FindDRFrequencies(Freq, Fluo, TeoDR, entorno=3):
+    """
+    Busca los indices y la frecuencia de los minimos en un entorno cercano al de la DR.
+    Si no encuentra, devuelve el valor teórico.
+    """
+
+    IndiceDRteo1, IndiceEntornoinicialDRteo1, IndiceEntornofinalDRteo1 = GetClosestIndex(Freq, TeoDR[0]), GetClosestIndex(Freq, TeoDR[0]-entorno), GetClosestIndex(Freq, TeoDR[0]+entorno)
+    IndiceDRteo2, IndiceEntornoinicialDRteo2, IndiceEntornofinalDRteo2 = GetClosestIndex(Freq, TeoDR[1]), GetClosestIndex(Freq, TeoDR[1]-entorno), GetClosestIndex(Freq, TeoDR[1]+entorno)
+    IndiceDRteo3, IndiceEntornoinicialDRteo3, IndiceEntornofinalDRteo3 = GetClosestIndex(Freq, TeoDR[2]), GetClosestIndex(Freq, TeoDR[2]-entorno), GetClosestIndex(Freq, TeoDR[2]+entorno)
+    IndiceDRteo4, IndiceEntornoinicialDRteo4, IndiceEntornofinalDRteo4 = GetClosestIndex(Freq, TeoDR[3]), GetClosestIndex(Freq, TeoDR[3]-entorno), GetClosestIndex(Freq, TeoDR[3]+entorno)
+    IndiceDRteo5, IndiceEntornoinicialDRteo5, IndiceEntornofinalDRteo5 = GetClosestIndex(Freq, TeoDR[4]), GetClosestIndex(Freq, TeoDR[4]-entorno), GetClosestIndex(Freq, TeoDR[4]+entorno)
+    IndiceDRteo6, IndiceEntornoinicialDRteo6, IndiceEntornofinalDRteo6 = GetClosestIndex(Freq, TeoDR[5]), GetClosestIndex(Freq, TeoDR[5]-entorno), GetClosestIndex(Freq, TeoDR[5]+entorno)
+
+    EntornoFreqDR1, EntornoFreqDR2 = Freq[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Freq[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
+    EntornoFreqDR3, EntornoFreqDR4 = Freq[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Freq[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
+    EntornoFreqDR5, EntornoFreqDR6 = Freq[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Freq[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
+
+    EntornoFluoDR1, EntornoFluoDR2 = Fluo[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Fluo[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
+    EntornoFluoDR3, EntornoFluoDR4 = Fluo[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Fluo[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
+    EntornoFluoDR5, EntornoFluoDR6 = Fluo[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Fluo[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
+
+    IndiceFluoMinimaEntorno1, IndiceFluoMinimaEntorno2 = argrelextrema(np.array(EntornoFluoDR1), np.less)[0], argrelextrema(np.array(EntornoFluoDR2), np.less)[0]
+    IndiceFluoMinimaEntorno3, IndiceFluoMinimaEntorno4 = argrelextrema(np.array(EntornoFluoDR3), np.less)[0], argrelextrema(np.array(EntornoFluoDR4), np.less)[0]
+    IndiceFluoMinimaEntorno5, IndiceFluoMinimaEntorno6 = argrelextrema(np.array(EntornoFluoDR5), np.less)[0], argrelextrema(np.array(EntornoFluoDR6), np.less)[0]
+
+    try:
+        FreqDR1 = EntornoFreqDR1[int(IndiceFluoMinimaEntorno1)]
+        IndiceDR1 = GetClosestIndex(Freq, FreqDR1)
+    except:
+        FreqDR1 = TeoDR[0]
+        IndiceDR1 = IndiceDRteo1
+    try:
+        FreqDR2 = EntornoFreqDR2[int(IndiceFluoMinimaEntorno2)]
+        IndiceDR2 = GetClosestIndex(Freq, FreqDR2)
+    except:
+        FreqDR2 = TeoDR[1]
+        IndiceDR2 = IndiceDRteo2
+    try:
+        FreqDR3 = EntornoFreqDR3[int(IndiceFluoMinimaEntorno3)]
+        IndiceDR3 = GetClosestIndex(Freq, FreqDR3)
+    except:
+        FreqDR3 = TeoDR[2]
+        IndiceDR3 = IndiceDRteo3
+    try:
+        FreqDR4 = EntornoFreqDR4[int(IndiceFluoMinimaEntorno4)]
+        IndiceDR4 = GetClosestIndex(Freq, FreqDR4)
+    except:
+        FreqDR4 = TeoDR[3]
+        IndiceDR4 = IndiceDRteo4
+    try:
+        FreqDR5 = EntornoFreqDR5[int(IndiceFluoMinimaEntorno5)]
+        IndiceDR5 = GetClosestIndex(Freq, FreqDR5)
+    except:
+        FreqDR5 = TeoDR[4]
+        IndiceDR5 = IndiceDRteo5
+    try:
+        FreqDR6 = EntornoFreqDR6[int(IndiceFluoMinimaEntorno6)]
+        IndiceDR6 = GetClosestIndex(Freq, FreqDR6)
+    except:
+        FreqDR6 = TeoDR[5]
+        IndiceDR6 = IndiceDRteo6
+
+    return [IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6], [FreqDR1, FreqDR2, FreqDR3, FreqDR4, FreqDR5, FreqDR6]
+
+
+def FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=-100, getindices=False):
+    """
+    Toma los indices donde estan las DR y evalua su fluorescencia. Esos indices son minimos locales en un entorno
+    cercano a las DR teoricas y, si no hay ningun minimo, toma la teorica.
+    Luego, hace el cociente de esa fluorescencia y un factor de normalización segun NormalizationCriterium:
+        1: Devuelve la fluorescencia absoluta de los minimos
+        2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no
+        hubiera una resonancia oscura y hubiera una recta. Ese valor esta a  DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
+        3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
+        o el valor limite por izquierda de la curva
+        4: Deuelve el cociente entre la fluorescencia del minimo y el valor de fluorescencia a detuning 0 MHz
+
+    """
+    IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6 = IndicesDR[0], IndicesDR[1], IndicesDR[2], IndicesDR[3], IndicesDR[4], IndicesDR[5]
+    FluorescenceOfMinimums = [Fluo[IndiceDR1], Fluo[IndiceDR2], Fluo[IndiceDR3], Fluo[IndiceDR4], Fluo[IndiceDR5], Fluo[IndiceDR6]]
+    FrequencyOfMinimums = [Freq[IndiceDR1], Freq[IndiceDR2], Freq[IndiceDR3], Freq[IndiceDR4], Freq[IndiceDR5], Freq[IndiceDR6]]
+    DistanciaFrecuenciaCociente = 25
+
+    if NormalizationCriterium==0:
+        print('che')
+        return FrequencyOfMinimums, FluorescenceOfMinimums
+
+    if NormalizationCriterium==1:
+        Fluorescenciacerodetuning = Fluo[GetClosestIndex(Freq, 0)]
+        Fluorescenciaasintotica = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
+        return FrequencyOfMinimums, np.array([Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica])
+
+    if NormalizationCriterium==2:
+        k = 0
+        while k < len(Freq):
+            if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
+                FluoIzquierda = Fluo[k]
+                indiceizquierda = k
+                print('Izq:', Freq[k])
+                break
+            else:
+                k = k + 1
+
+        l = 0
+        while l < len(Freq):
+            if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
+                FluoDerecha = Fluo[l]
+                indicederecha = l
+                print('Der: ', Freq[l])
+                break
+            else:
+                l = l + 1
+        FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
+        print(FluoNormDivisor)
+
+    if NormalizationCriterium==3:
+        #asintotico
+        FluoNormDivisor = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
+
+    if NormalizationCriterium==4:
+        #este te tira la fluorescencia de detuning 0
+        FluoNormDivisor = Fluo[GetClosestIndex(Freq, 0)]
+
+
+    RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
+    print('Esto: ', RelativeFluorescenceOfMinimums)
+
+    if NormalizationCriterium==2 and getindices==True:
+        return FrequencyOfMinimums, RelativeFluorescenceOfMinimums, indiceizquierda, indicederecha
+
+    return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
+
+
+def GetFinalMaps(MapasDR1, MapasDR2, MapasDR3, MapasDR4, MapasDR5, MapasDR6):
+    """
+    Nota: esto vale para polarizacion del 397 sigma+ + sigma-. Sino hay que cambiar los coeficientes.
+    La estructura es:
+        MapasDRi = [MapaMedido_criterio1_DRi, MapaMedido_criterio2_DRi, MapaMedido_criterio3_DRi, MapaMedido_criterio4_DRi]
+    """
+    Mapa1 = MapasDR1[0]
+
+    Mapa2pi = np.sqrt(3)*(MapasDR2[1] + MapasDR5[1])
+    Mapa2smas = np.sqrt(12/2)*MapasDR3[1] + (2/np.sqrt(2))*MapasDR6[1]
+    Mapa2smenos = (2/np.sqrt(2))*MapasDR1[1] + np.sqrt(12/2)*MapasDR4[1]
+
+    Mapa3pi = np.sqrt(3)*(MapasDR2[2] + MapasDR5[2])
+    Mapa3smas = np.sqrt(12/2)*MapasDR3[2] + (2/np.sqrt(2))*MapasDR6[2]
+    Mapa3smenos = (2/np.sqrt(2))*MapasDR1[2] + np.sqrt(12/2)*MapasDR4[2]
+
+    return Mapa1, [Mapa2pi, Mapa2smas, Mapa2smenos], [Mapa3pi, Mapa3smas, Mapa3smenos]
+
+
+def CombinateDRwithCG(RelMinMedido1, RelMinMedido2, RelMinMedido3, RelMinMedido4):
+    Fluo1 = RelMinMedido1[0]
+
+    Fluo2pi = np.sqrt(3)*(RelMinMedido2[1] + RelMinMedido2[4])
+    Fluo2smas = np.sqrt(12/2)*RelMinMedido2[2] + (2/np.sqrt(2))*RelMinMedido2[5]
+    Fluo2smenos = (2/np.sqrt(2))*RelMinMedido2[0] + np.sqrt(12/2)*RelMinMedido2[3]
+
+    Fluo3pi = np.sqrt(3)*(RelMinMedido3[1] + RelMinMedido3[4])
+    Fluo3smas = np.sqrt(12/2)*RelMinMedido3[2] + (2/np.sqrt(2))*RelMinMedido3[5]
+    Fluo3smenos = (2/np.sqrt(2))*RelMinMedido3[0] + np.sqrt(12/2)*RelMinMedido3[3]
+
+    return Fluo1, [Fluo2pi, Fluo2smas, Fluo2smenos], [Fluo3pi, Fluo3smas, Fluo3smenos]
+
+
+
+def IdentifyPolarizationCoincidences(theoricalmap, target, tolerance=1e-1):
+    """
+    Busca en un mapa 2D la presencia de un valor target (medido) con tolerancia tolerance.
+    Si lo encuentra, pone un 1. Sino, un 0. Al plotear con pcolor se verá
+    en blanco la zona donde el valor medido se puede hallar.
+    """
+    CoincidenceMatrix = np.zeros((len(theoricalmap), len(theoricalmap[0])))
+
+    i = 0
+    while i<len(theoricalmap):
+        j = 0
+        while j<len(theoricalmap[0]):
+            if abs(theoricalmap[i][j]-target) < tolerance:
+                CoincidenceMatrix[i][j] = 1
+            j=j+1
+        i=i+1
+    return CoincidenceMatrix
+
+def RetrieveAbsoluteCoincidencesBetweenMaps(MapsVectors):
+
+    MatrixSum = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
+    AbsoluteCoincidencesMatrix = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
+
+
+    MatrixMapsVectors = []
+
+    for i in range(len(MapsVectors)):
+        MatrixMapsVectors.append(np.matrix(MapsVectors[i]))
+
+    for i in range(len(MatrixMapsVectors)):
+        MatrixSum = MatrixSum + MatrixMapsVectors[i]
+
+    MaxNumberOfCoincidences = np.max(MatrixSum)
+
+    ListMatrixSum = [list(i) for i in list(np.array(MatrixSum))]
+
+    for i in range(len(ListMatrixSum)):
+        for j in range(len(ListMatrixSum[0])):
+            if ListMatrixSum[i][j] == MaxNumberOfCoincidences:
+                AbsoluteCoincidencesMatrix[i][j] = 1
+
+    return AbsoluteCoincidencesMatrix, MaxNumberOfCoincidences
+
+
+def MeasureMeanValueOfEstimatedArea(AbsoluteCoincidencesMap, X, Y):
+    NonZeroIndices = np.nonzero(AbsoluteCoincidencesMap)
+
+    Xsum = 0
+    Xvec = []
+    Ysum = 0
+    Yvec = []
+    N = len(NonZeroIndices[0])
+    for i in range(N):
+        Xsum = Xsum + X[NonZeroIndices[1][i]]
+        Xvec.append(X[NonZeroIndices[1][i]])
+        Ysum = Ysum + Y[NonZeroIndices[0][i]]
+        Yvec.append(Y[NonZeroIndices[0][i]])
+    Xaverage = Xsum/N
+    Yaverage = Ysum/N
+    Xspread = np.std(Xvec)
+    Yspread = np.std(Yvec)
+    return Xaverage, Yaverage, N, Xspread, Yspread
+
+
+def MeasureRelativeFluorescenceFromCPT(Freq, Fluo, u, titadoppler, detuningrepump, detuningdoppler, frefasint=-100, entorno=3):
+    ResonanciasTeoricas, ResonanciasPositivas = CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)
+    IndicesDR, FreqsDR = FindDRFrequencies(Freq, Fluo, ResonanciasTeoricas, entorno=entorno)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums0 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=0, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums1 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums2, indiceizquierda, indicederecha = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=2, frecuenciareferenciacriterioasintotico=frefasint, getindices=True)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums3 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=3, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums4 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=4, frecuenciareferenciacriterioasintotico=frefasint)
+    print('hola')
+    print(RelativeFluorescenceOfMinimums0)
+    return RelativeFluorescenceOfMinimums0, RelativeFluorescenceOfMinimums1, RelativeFluorescenceOfMinimums2, RelativeFluorescenceOfMinimums3, RelativeFluorescenceOfMinimums4, IndicesDR, [indiceizquierda, indicederecha]
+
+
+def GenerateNoisyCPT(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+
+def GenerateNoisyCPT_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+def GenerateNoisyCPT_fit(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqs, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, min(freqs), max(freqs) + freqs[1]-freqs[0], freqs[1]-freqs[0], plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+
+def AddNoiseToCPT(Fluo, noisefactor):
+    return [f+noisefactor*(2*random.random()-1) for f in Fluo]
+
+
+def SmoothNoisyCPT(Fluo, window=11, poly=3):
+    SmoothenFluo = sf(Fluo, window, poly)
+    return SmoothenFluo
+
+
+def GetMinimaInfo(Freq, Fluo, u, titadoppler, detuningdoppler, detuningrepump, MinimumCriterium=2, NormalizationCriterium=1):
+    """
+    FUNCION VIEJA
+
+    Esta funcion devuelve valores de frecuencias y fluorescencia relativa de los minimos.
+    Minimumcriterion:
+        1: Saca los minimos con funcion argelextrema
+        2: Directamente con las frecuencias teoricas busca las fluorescencias
+
+    Normalizationcriterium:
+        1: Devuelve la fluorescencia absoluta de los minimos
+        2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no
+        hubiera una resonancia oscura y hubiera una recta. Ese valor esta a  DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
+        3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
+        o el valor limite por izquierda de la curva
+
+    """
+    FluorescenceOfMaximum = max(Fluo)
+    FrequencyOfMaximum = Freq[Fluo.index(FluorescenceOfMaximum)]
+
+
+    #criterio para encontrar los minimos
+    #criterio usando minimos de la fluorescencia calculados con la curva
+    if MinimumCriterium == 1:
+        LocationOfMinimums = argrelextrema(np.array(Fluo), np.less)[0]
+        FluorescenceOfMinimums = np.array([Fluo[i] for i in LocationOfMinimums])
+        FrequencyOfMinimums = np.array([Freq[j] for j in LocationOfMinimums])
+
+
+    #criterio con las DR teoricas
+    if MinimumCriterium == 2:
+        FrecuenciasDRTeoricas, FrecuenciasDRTeoricasPositivas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
+        FrequencyOfMinimums = []
+        FluorescenceOfMinimums =[]
+        print(FrecuenciasDRTeoricas)
+        k=0
+
+        ventanita = 0.001
+        while k < len(Freq):
+            if Freq[k] < FrecuenciasDRTeoricas[0] + ventanita and Freq[k] > FrecuenciasDRTeoricas[0] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[1] + ventanita and Freq[k] > FrecuenciasDRTeoricas[1] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[2] + ventanita and Freq[k] > FrecuenciasDRTeoricas[2] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[3] + ventanita and Freq[k] > FrecuenciasDRTeoricas[3] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[4] + ventanita and Freq[k] > FrecuenciasDRTeoricas[4] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[5] + ventanita and Freq[k] > FrecuenciasDRTeoricas[5] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            k = k + 1
+        print(FrequencyOfMinimums)
+        if len(FrequencyOfMinimums) !=  len(FrecuenciasDRTeoricas):
+            print('NO ANDA BIEN ESTO PAPI, revisalo')
+
+
+    #esto es para establecer un criterio para la fluorescencia relativa
+    DistanciaFrecuenciaCociente = 15
+
+    if NormalizationCriterium==1:
+        FluoNormDivisor = 1
+
+    if NormalizationCriterium==2:
+        k = 0
+        while k < len(Freq):
+            if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
+                FluoIzquierda = Fluo[k]
+                print('Izq:', Freq[k])
+                break
+            else:
+                k = k + 1
+
+        l = 0
+        while l < len(Freq):
+            if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
+                FluoDerecha = Fluo[l]
+                print('Der: ', Freq[l])
+                break
+            else:
+                l = l + 1
+        FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
+        print(FluoNormDivisor)
+
+    if NormalizationCriterium==3:
+        FluoNormDivisor = Fluo[0]
+
+    RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
+
+    return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
+
+
+def GetPlotsofFluovsAngle_8levels(FrequencyOfMinimumsVector, RelativeFluorescenceOfMinimumsVector, u, titadoppler, detuningdoppler, detuningrepump, ventana=0.25, taketheoricalDR=False):
+
+    #primero buscamos las frecuencias referencia que se parezcan a las 6:
+    i = 0
+    FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[0]
+
+    FrecuenciasDRTeoricas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
+
+    while i < len(FrequencyOfMinimumsVector):
+        if len(FrequencyOfMinimumsVector[i])==len(FrecuenciasDRTeoricas):
+            FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[i]
+            print('Cool! Taking the DR identified with any curve')
+            break
+        else:
+            i = i + 1
+            if i==len(FrequencyOfMinimumsVector):
+                print('No hay ningun plot con 5 resonancias oscuras. Tomo las teóricas')
+                FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
+
+    if taketheoricalDR:
+        FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
+
+
+
+    Ventana = abs(ventana*(FrecuenciasReferenciaBase[1] - FrecuenciasReferenciaBase[0])) #ventana separadora de resonancias
+
+    print('Ventana = ', Ventana)
+    DarkResonance1Frequency = []
+    DarkResonance1Fluorescence = []
+
+    DarkResonance2Frequency = []
+    DarkResonance2Fluorescence = []
+
+    DarkResonance3Frequency = []
+    DarkResonance3Fluorescence = []
+
+    DarkResonance4Frequency = []
+    DarkResonance4Fluorescence = []
+
+    DarkResonance5Frequency = []
+    DarkResonance5Fluorescence = []
+
+    DarkResonance6Frequency = []
+    DarkResonance6Fluorescence = []
+
+    i = 0
+
+
+    while i < len(FrequencyOfMinimumsVector):
+        j = 0
+        FrecuenciasReferencia = [i for i in FrecuenciasReferenciaBase]
+
+
+        while j < len(FrequencyOfMinimumsVector[i]):
+
+
+            if abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[0])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[0])-Ventana):
+                DarkResonance1Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance1Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[0] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[1])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[1])-Ventana):
+                DarkResonance2Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance2Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[1] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[2])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[2])-Ventana):
+                DarkResonance3Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance3Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[2] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[3])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[3])-Ventana):
+                DarkResonance4Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance4Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[3] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[4])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[4])-Ventana):
+                DarkResonance5Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance5Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[4] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[5])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[5])-Ventana):
+                DarkResonance6Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance6Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[5] = 0
+
+            else:
+                #print('Algo anduvo mal, por ahi tenes que cambiar la ventana che')
+                pass
+            j = j + 1
+
+        if np.count_nonzero(FrecuenciasReferencia) > 0:
+            if FrecuenciasReferencia[0] != 0:
+                DarkResonance1Frequency.append(FrecuenciasReferencia[0])
+                DarkResonance1Fluorescence.append()
+
+            if FrecuenciasReferencia[1] != 0:
+                DarkResonance2Frequency.append(FrecuenciasReferencia[1])
+                DarkResonance2Fluorescence.append(0)
+
+            if FrecuenciasReferencia[2] != 0:
+                DarkResonance3Frequency.append(FrecuenciasReferencia[2])
+                DarkResonance3Fluorescence.append(0)
+
+            if FrecuenciasReferencia[3] != 0:
+                DarkResonance4Frequency.append(FrecuenciasReferencia[3])
+                DarkResonance4Fluorescence.append(0)
+
+            if FrecuenciasReferencia[4] != 0:
+                DarkResonance5Frequency.append(FrecuenciasReferencia[4])
+                DarkResonance5Fluorescence.append(0)
+
+            if FrecuenciasReferencia[5] != 0:
+                DarkResonance6Frequency.append(FrecuenciasReferencia[5])
+                DarkResonance6Fluorescence.append(0)
+
+        i = i + 1
+
+    return DarkResonance1Frequency, DarkResonance1Fluorescence, DarkResonance2Frequency, DarkResonance2Fluorescence, DarkResonance3Frequency, DarkResonance3Fluorescence, DarkResonance4Frequency, DarkResonance4Fluorescence, DarkResonance5Frequency, DarkResonance5Fluorescence, DarkResonance6Frequency, DarkResonance6Fluorescence, FrecuenciasReferenciaBase
+
+
+
+def PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
+
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+    Fluovectors = []
+
+    for titaprobe in titaprobeVec:
+
+        tinicial = time.time()
+
+        ProbeDetuningVectorL, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
+
+        tfinal = time.time()
+
+        print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")
+
+        if plot:
+            plt.figure()
+            plt.xlabel('Repump detuning (MHz')
+            plt.ylabel('Fluorescence (A.U.)')
+            plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
+            plt.legend()
+
+        Fluovectors.append(Fluovector)
+
+        if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
+            Fluovectors = Fluovector
+
+    return ProbeDetuningVectorL, Fluovectors
+
+
+
+def PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
+
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+    Fluovectors = []
+
+    for titaprobe in titaprobeVec:
+
+        tinicial = time.time()
+
+        ProbeDetuningVectorL, Fluovector = CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
+
+        tfinal = time.time()
+
+        print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")
+
+        if plot:
+            plt.figure()
+            plt.xlabel('Repump detuning (MHz')
+            plt.ylabel('Fluorescence (A.U.)')
+            plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
+            plt.legend()
+
+        Fluovectors.append(Fluovector)
+
+        if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
+            Fluovectors = Fluovector
+
+    return ProbeDetuningVectorL, Fluovectors

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/lolo_modelo_full_8niveles.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+20 de dic 2023
+
+@author: lolo
+
+reingenieria del código que anda
+
+Surge de fusionar
+Data/EITfit/MM_eightLevel_2repumps_AnalysisFunctions.py
+Data/EITfit/MM_eightLevel_2repumps_python_scripts.py
+
+
+MAPA de FUNCIONES
+
+CPTspectrum8levels_MM
+ |--> FullL_MM => ndarray(64,64,np.complex_)
+       |--> dopplerBroadening => float
+       |--> EffectiveL => ndarray(8,8,np.complex_)
+       |--> H0matrix => ndarray(8,8,np.complex_)
+       |--> HImatrix => ndarray(8,8,np.complex_)
+       |--> CalculateSingleMmatrix => ndarray(64,64,np.complex_)
+       |--> LtempCalculus => ndarray,ndarray
+       |--> GetL1 => ndarray
+
+"""
+
+# pylint: disable=C0301,R0913,R0914,W0621
+
+import time
+import random
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.signal import savgol_filter as sf
+
+from  numba import jit,njit
+
+@njit
+def PerformExperiment_8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None):
+    """
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    rta = CPTspectrum8levels_MM(sg, sp, gPS, gPD, DetDoppler, u,
+                                DopplerLaserLinewidth, ProbeLaserLinewidth,
+                                T, alpha, phidoppler, titadoppler, phiprobe,
+                                titaprobe, circularityprobe, beta, drivefreq,
+                                freqMin=freqMin, freqMax=freqMax, freqStep=freqStep,
+                                plot=False, solvemode=1)
+    ProbeDetuningVectorL, Fluovector = rta
+    return ProbeDetuningVectorL, Fluovector
+
+
+def GenerateNoisyCPT_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, kg, kr, v0, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    """
+    Genera un resultado de PerformExperiment_8levels_MM con ruido normal agregado
+    """
+    nFrequencyvector, Fluovector = PerformExperiment_8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, kg, kr, v0, drivefreq, freqMin, freqMax, freqStep, circularityprobe, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return nFrequencyvector, NoisyFluovector
+
+def GenerateNoisyCPT_MM_fit(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, beta, drivefreq, freqs, circularityprobe=1, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    """
+    Este no se qué hace
+    """
+    Frequencyvector, Fluovector = PerformExperiment_8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, beta, drivefreq, freqs[0], freqs[-1], freqs[1]-freqs[0], circularityprobe, plot=False, solvemode=1, detpvec=None)
+    #NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, Fluovector
+
+
+def SmoothNoisyCPT(Fluo, window=11, poly=3):
+    SmoothenFluo = sf(Fluo, window, poly)
+    return SmoothenFluo
+
+
+
+#%% Estos son los auxiliares ###################################################
+
+"""
+Esta parte es la del modelo
+"""
+
+@njit
+def make_diag(vec):
+    "Construye matris diagonal desde una lista o vector"
+    return np.eye(len(vec))*np.array(vec).reshape(-1,1)
+
+# @njit
+# def kron(a,b):
+#     "Hago el producto de Kronecker a mano"
+#     return np.vstack( [  np.hstack( [ a[k,j]*b for j in range(a.shape[1]) ] )  for k in range(a.shape[0])])
+
+# @njit
+# def kron(A, B):
+#     cola = A.shape[1]
+#     rowa = A.shape[0]
+#     colb = B.shape[1]
+#     rowb = B.shape[0]
+#
+#     C = [[0] * (cola * colb)  for _ in range(rowa * rowb) ]
+#
+#     for i in range(rowa):
+#         for k in range(cola):
+#             for j in range(rowb):
+#                 for l in range(colb):
+#                     C[i * rowb + k][j * colb + l]  = A[i][j] * B[k][l]
+#     return np.array(C)
+import numba
+
+@jit
+def kron(A,B):
+    out=np.empty((A.shape[0],B.shape[0],A.shape[1],B.shape[1]),dtype=A.dtype)
+    for i in numba.prange(A.shape[0]):
+        for j in range(B.shape[0]):
+            for k in range(A.shape[1]):
+                for l in range(B.shape[1]):
+                    out[i,j,k,l]=A[i,k]*B[j,l]
+    return out
+
+
+@njit
+def H0matrix(Detg, Detp, u):
+    """
+    Calcula la matriz H0 en donde dr es el detuning del doppler, dp es el retuning
+    del repump y u es el campo magnético en Hz/Gauss.
+    Para esto se toma la energía del nivel P como 0
+    """
+
+    eigenEnergies = (Detg-u, Detg+u, -u/3, u/3, Detp-6*u/5,
+                    Detp-2*u/5, Detp+2*u/5, Detp+6*u/5)
+    #pagina 26 de Oberst. los lande del calcio son iguales a Bario.
+    # H0 = np.diag(eigenEnergies)
+    # H0 = np.eye(len(eigenEnergies))*np.array(eigenEnergies).reshape(-1,1)
+    H0 = make_diag(eigenEnergies)
+    return H0
+
+@njit
+def HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe=1):
+    """
+    Calcula la matriz de interacción Hsp + Hpd, en donde rabR es la frecuencia de rabi de
+    la transición Doppler SP, rabP es la frecuencia de rabi de la transición repump DP,
+    y las componentes ei_r y ei_p son las componentes de la polarización
+    del campo eléctrico incidente de doppler y repump respectivamente. Deben
+    estar normalizadas a 1
+    """
+    HI = np.zeros((8, 8), dtype=np.complex_)
+
+
+    i, j = 1, 3
+    HI[i-1, j-1] = (rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 1, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 2, 3
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(-1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 2, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 3, 5
+    HI[i-1, j-1] = -(rabP/2) * np.sin(titaprobe)*(np.cos(phiprobe)-1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 3, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 3, 7
+    HI[i-1, j-1] = rabP/np.sqrt(12) * np.sin(titaprobe)*(np.cos(phiprobe)+1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+
+    i, j = 4, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(12)) * np.sin(titaprobe)*(np.cos(phiprobe)-1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 4, 7
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+    i, j = 4, 8
+    HI[i-1, j-1] = (rabP/2) * np.sin(titaprobe)*(np.cos(phiprobe)+1j*np.sin(phiprobe)*circularityprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+
+
+    return HI
+
+# @jit
+# def LtempCalculus(beta:float, drivefreq:float, forma=1):
+#     Hint = np.zeros((8, 8), dtype=np.complex_)
+
+#     ampg=beta*drivefreq
+#     ampr=beta*drivefreq*(397/866)
+#     Hint[0,0] = ampg
+#     Hint[1,1] = ampg
+#     Hint[4,4] = ampr
+#     Hint[5,5] = ampr
+#     Hint[6,6] = ampr
+#     Hint[7,7] = ampr
+
+#     if forma==1:
+#         Ltemp = np.zeros((64, 64), dtype=np.complex_)
+
+#         for r in range(8):
+#             for q in range(8):
+#                 if r!=q:
+#                     Ltemp[r*8+q][r*8+q] = (-1j)*(Hint[r,r] - Hint[q,q])
+
+#     if forma==2:
+#         # deltaKro = np.diag([1, 1, 1, 1, 1, 1, 1, 1])
+#         deltaKro = make_diag([1., 1., 1., 1., 1., 1., 1., 1.]).astype(np.complex_)
+#         # Ltemp = (-1j)*(np.kron(Hint, deltaKro) - np.kron(deltaKro, Hint))
+#         Ltemp = (-1j)*(kron(Hint, deltaKro) - kron(deltaKro, Hint))
+
+#     Omega = np.zeros((64, 64), dtype=np.complex_)
+
+#     for i in range(64):
+#         Omega[i, i] = (1j)*drivefreq
+
+#     return Ltemp, Omega
+
+@njit
+def LtempCalculus(beta:float, drivefreq:float, forma=1):
+    Hint = np.zeros((8, 8), dtype=np.complex_)
+
+    ampg=beta*drivefreq
+    ampr=beta*drivefreq*(397/866)
+    Hint[0,0] = ampg
+    Hint[1,1] = ampg
+    Hint[4,4] = ampr
+    Hint[5,5] = ampr
+    Hint[6,6] = ampr
+    Hint[7,7] = ampr
+
+    Ltemp = np.zeros((64, 64), dtype=np.complex_)
+
+    for r in range(8):
+        for q in range(8):
+            if r!=q:
+                Ltemp[r*8+q][r*8+q] = (-1j)*(Hint[r,r] - Hint[q,q])
+
+    Omega = np.zeros((64, 64), dtype=np.complex_)
+
+    for i in range(64):
+        Omega[i, i] = (1j)*drivefreq
+
+    return Ltemp, Omega
+
+# LtempCalculus(0,1)
+# raise ValueError('aaa')
+
+
+
+
+
+
+@njit
+def GetL1(Ltemp, L0, Omega, nmax):
+    """
+    Devuelve Splus0 y Sminus0
+    """
+
+    # Sp = (-1)*(np.matrix(np.linalg.inv(L0 - (nmax+1)*Omega))*0.5*np.matrix(Ltemp))
+    # Sm = (-1)*(np.matrix(np.linalg.inv(L0 + (nmax+1)*Omega))*0.5*np.matrix(Ltemp))
+
+    Sp = (-1)*np.linalg.inv(L0 - (nmax+1)*Omega).dot(0.5*Ltemp)
+    Sm = (-1)*np.linalg.inv(L0 + (nmax+1)*Omega).dot(0.5*Ltemp)
+
+    for n in list(range(nmax+1))[(nmax+1)::-1][0:len(list(range(nmax+1))[(nmax+1)::-1])-1]: #jaja esto solo es para que vaya de nmax a 1 bajando. debe haber algo mas facil pero kcio
+        # Sp = (-1)*(np.matrix(np.linalg.inv(L0 - n*Omega + (0.5*Ltemp*np.matrix(Sp))))*0.5*np.matrix(Ltemp))
+        # Sm = (-1)*(np.matrix(np.linalg.inv(L0 + n*Omega + (0.5*Ltemp*np.matrix(Sm))))*0.5*np.matrix(Ltemp))
+
+        Sp = (-1)*np.linalg.inv(L0 - n*Omega + (0.5*Ltemp.dot(Sp))).dot(0.5*Ltemp)
+        Sm = (-1)*np.linalg.inv(L0 + n*Omega + (0.5*Ltemp.dot(Sm))).dot(0.5*Ltemp)
+
+    # L1 = 0.5*np.matrix(Ltemp)*(np.matrix(Sp) + np.matrix(Sm))
+    L1 = 0.5*Ltemp.dot(Sp + Sm)
+
+    return L1
+
+@njit
+def EffectiveL(gPS, gPD, lwg, lwp):
+
+    """
+    Siendo Heff = H + EffectiveL, calcula dicho EffectiveL que es (-0.5j)*sumatoria(CmDaga*Cm) que luego sirve para calcular el Liouvilliano
+    """
+    Leff = np.zeros((8, 8), dtype=np.complex_)
+
+    Leff[0, 0] = 2*lwg
+    Leff[1, 1] = 2*lwg
+    Leff[2, 2] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[3, 3] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[4, 4] = 2*lwp
+    Leff[5, 5] = 2*lwp
+    Leff[6, 6] = 2*lwp
+    Leff[7, 7] = 2*lwp
+
+    return (-0.5j)*Leff
+
+@njit
+def CalculateSingleMmatrix(gPS, gPD, lwg, lwp):
+    """
+    Si tomamos el Liuvilliano como L = (-j)*(Heff*deltak - Heffdaga*deltak) + sum(Mm),
+    esta funcion calcula dichos Mm, que tienen dimensión 64x64 ya que esa es la dimensión del L. Estas componentes
+    salen de hacer la cuenta a mano conociendo los Cm y considerando que Mm[8*(r-1)+s, 8*(k-1)+j] = Cm[r,l] + Cmdaga[j,s] = Cm[r,l] + Cm[s,j]
+    ya que los componentes de Cm son reales.
+    Esta M es la suma de las 8 matrices M.
+
+    """
+
+    M = np.zeros((64, 64), dtype=np.complex_)
+
+    M[0,27] = (2/3)*gPS
+
+    M[9,18] = (2/3)*gPS
+
+    M[0,18] = (1/3)*gPS
+    M[1,19] = -(1/3)*gPS
+    M[8,26] = -(1/3)*gPS
+    M[9,27] = (1/3)*gPS
+
+    M[36,18] = (1/2)*gPD
+    M[37,19] = (1/np.sqrt(12))*gPD
+    M[44,26] = (1/np.sqrt(12))*gPD
+    M[45,27] = (1/6)*gPD
+
+    M[54,18] = (1/6)*gPD
+    M[55,19] = (1/np.sqrt(12))*gPD
+    M[62,26] = (1/np.sqrt(12))*gPD
+    M[63,27] = (1/2)*gPD
+
+    M[45,18] = (1/3)*gPD
+    M[46,19] = (1/3)*gPD
+    M[53,26] = (1/3)*gPD
+    M[54,27] = (1/3)*gPD
+
+    M[0,0] = 2*lwg
+    M[1,1] = 2*lwg
+    M[8,8] = 2*lwg
+    M[9,9] = 2*lwg
+
+    for k in [36, 37, 38, 39, 44, 45, 46, 47, 52, 53, 54, 55, 60, 61, 62, 63]:
+        M[k,k]=2*lwp
+
+    return M
+
+@njit
+def dopplerBroadening(wlg, wlp, alpha, T, mcalcio = 6.655e-23*1e-3):
+    """
+    Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
+    wlg es la longitud de onda doppler, wlp la longitud de onda repump, T la temperatura del ion en kelvin, y alpha (en rads) el ángulo
+    que forman ambos láseres.
+    """
+
+    kboltzmann = 1.38e-23 #J/K
+
+    gammaD = (2*np.pi)*np.sqrt((1/(wlg*wlg)) + (1/(wlp*wlp)) - 2*(1/(wlg*wlp))*np.cos(alpha))*np.sqrt(kboltzmann*T/(1*mcalcio))
+
+    return gammaD
+
+@njit
+def FullL_MM(rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp = 0,
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, beta=0,
+          drivefreq=2*np.pi*22.135*1e6, T = 0, alpha = 0, circularityprobe=1):
+
+    """
+    Calcula el Liouvilliano total de manera explícita índice a índice. Suma aparte las componentes de las matrices M.
+    Es la más eficiente hasta ahora.
+    """
+
+
+    db = dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)
+
+    lwg = np.sqrt(lwg**2 + (0.83*db)**2)
+    lwp = np.sqrt(lwp**2 + (0.17*db)**2)
+
+    CC = EffectiveL(gPS, gPD, lwg, lwp)
+
+    Heff = H0matrix(Detg, Detp, u) + HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe) + CC
+
+    # Heffdaga = np.matrix(Heff).getH()
+    Heffdaga = np.conj(np.transpose(Heff))
+
+    Lfullpartial = np.zeros((64, 64), dtype=np.complex_)
+
+    for r  in range(8):
+        for q in range(8):
+            for k in range(8):
+                for j in range(8):
+                    if j!=q and r!=k:
+                        pass
+                    elif j==q and r!=k:
+                        if (r < 2 and k > 3) or (k < 2 and r > 3) or (r > 3 and k > 3) or (r==0 and k==1) or (r==1 and k==0) or (r==2 and k==3) or (r==3 and k==2):
+                            #todo esto sale de analizar explicitamente la matriz y tratar de no calcular cosas de más que dan cero
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k])
+                    elif j!=q and r==k:
+                        if (j < 2 and q > 3) or (q < 2 and j > 3) or (j > 3 and q > 3) or (j==0 and q==1) or (j==1 and q==0) or (j==2 and q==3) or (j==3 and q==2):
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(-Heffdaga[j,q])
+                    else:
+                        if Heff[r,k] == Heffdaga[j,q]:
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k]-Heffdaga[j,q])
+
+    M = CalculateSingleMmatrix(gPS, gPD, lwg, lwp)
+    # L0 = np.array(np.matrix(Lfullpartial) + M)
+    L0 = Lfullpartial + M
+
+    #ESTA PARTE ES CUANDO AGREGAS MICROMOCION
+    nmax = 3
+
+    #print(nmax)
+    Ltemp, Omega = LtempCalculus(beta, drivefreq)
+    #print(factor)
+    L1 = GetL1(Ltemp, L0, Omega, nmax)
+    Lfull = L0 + L1 #ESA CORRECCION ESTA EN L1
+    #HASTA ACA
+
+    #NORMALIZACION DE RHO
+    i = 0
+    while i < 64:
+        if i%9 == 0:
+            Lfull[0, i] = 1
+        else:
+            Lfull[0, i] = 0
+        i = i + 1
+
+    return Lfull
+
+
+"""
+Scripts para correr un experimento y hacer el análisis de los datos
+"""
+
+@njit
+def CPTspectrum8levels_MM(sg, sp, gPS, gPD, Detg, u, lwg, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, Circularityprobe, beta, drivefreq, freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+
+    """
+    ESTA ES LA FUNCION QUE ESTAMOS USANDO
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6+0*freqStep*1e6, freqStep*1e6)
+    Detg = 2*np.pi*Detg*1e6
+    #lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+    lwg, lwp = lwg*1e6, lwp*1e6
+
+    rabG = sg*gPS
+    rabP = sp*gPD
+    #u = 2*np.pi*u*1e6
+
+    Fluovector = []
+
+    # tinicial = time.time()
+
+    for Detp in DetProbeVector:
+        L = FullL_MM(rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp, phidoppler, titadoppler, phiprobe, titaprobe, beta, drivefreq, Temp, alpha, Circularityprobe)
+        # if solvemode == 1:
+        rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)],dtype=np.complex_))
+        Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27]))
+        Fluovector.append(Fluo)
+
+        # if solvemode == 2:
+        #     Linv = np.linalg.inv(L)
+        #     rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+        #     Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+        #     Fluovector.append(Fluo)
+
+    # tfinal = time.time()
+
+    # print('Done, Total time: ', round((tfinal-tinicial), 2), "s")
+
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+
+    # if plot:
+    #     plt.xlabel('Probe detuning (MHz)')
+    #     plt.ylabel('Fluorescence (A.U.)')
+    #     plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    #     plt.legend()
+
+    return DetProbeVectorMHz, Fluovector
+
+
+# @njit
+# def lolo():
+#     L = FullL_MM(100,200,12,123,14)
+#     return np.linalg.solve(L, np.array([int(i==0) for i in range(64)],dtype=np.complex_))
+
+# lolo()
+
+# raise ValueError('áaa')
+
+#%%
+if __name__ == "__main__":
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+
+    B = 25 #campo magnetico en gauss
+
+    u = c*B
+
+
+    sg, sr, sp = 0.5, 1.5, 4 #parámetros de saturación del doppler y repump
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+    rabG, rabR, rabP = sg*gPS, sr*gPD, sp*gPD #frecuencias de rabi
+
+    lwg, lwr, lwp = 0.3, 0.3, 0.3 #ancho de linea de los laseres
+
+    Detg = -25
+    Detr = 20 #detuning del doppler y repump
+
+    Temp = 0.0e-3 #temperatura en K
+    alpha = 0*(np.pi/180) #angulo entre los láseres
+
+    phidoppler, titadoppler = 0, 90
+    phirepump, titarepump = 0, 90
+
+    phiprobe, titaprobe = 0, 90
+
+
+    plotCPT = False
+
+    freqMin = -50
+    freqMax = 50
+    freqStep = 5e-2
+
+    # Frequencyvector, Fluovector = CPTspectrum8levels_MM(rabG, rabR, rabP, gPS, gPD, Detg, Detr,   u,  lwg,   lwr,        lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe,    phirepump, titarepump,             freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=plotCPT, solvemode=1)
+    Frequencyvector, Fluovector = PerformExperiment_8levels_MM(0.9,6.2,135591138.92893547,8482300.164692441,-24.5,32500000.0,0.1,0.1,0.001,0,0,90,0,90,2.0,139078306.77442014,-54.39999999999998,26.26666666666671,0.6666666666666856,circularityprobe=1,plot=False,solvemode=1,detpvec=None)
+
+
+    plt.plot(Frequencyvector, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    plt.xlabel('Probe detuning (MHz)')
+    plt.ylabel('Fluorescence (A.U.)')

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/threeLevel_2repumps_AnalysisFunctions.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jul  2 16:30:09 2020
+
+@author: oem
+"""
+
+
+import os
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+from scipy.signal import argrelextrema
+from EITfit.threeLevel_2repumps_linealpol_python_scripts import CPTspectrum8levels, CPTspectrum8levels_fixedRabi
+import random
+from scipy.signal import savgol_filter as sf
+
+
+def CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump):
+    
+    if titadoppler==0:
+        NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
+    elif titadoppler==90:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+    else:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+    
+    PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
+    
+    return [detuningdoppler + dr for dr in NegativeDR], [detuningrepump + dr for dr in PositiveDR]
+
+
+def GetClosestIndex(Vector, value, tolerance=1e-3):
+
+    i = 0
+    while i<len(Vector):
+        if abs(Vector[i] - value) < tolerance:
+            return i
+        else:
+            i = i + 1
+    return GetClosestIndex(Vector, value, tolerance=2*tolerance)
+    
+
+def FindDRFrequencies(Freq, Fluo, TeoDR, entorno=3):
+    """
+    Busca los indices y la frecuencia de los minimos en un entorno cercano al de la DR.
+    Si no encuentra, devuelve el valor teórico.
+    """
+    
+    IndiceDRteo1, IndiceEntornoinicialDRteo1, IndiceEntornofinalDRteo1 = GetClosestIndex(Freq, TeoDR[0]), GetClosestIndex(Freq, TeoDR[0]-entorno), GetClosestIndex(Freq, TeoDR[0]+entorno)
+    IndiceDRteo2, IndiceEntornoinicialDRteo2, IndiceEntornofinalDRteo2 = GetClosestIndex(Freq, TeoDR[1]), GetClosestIndex(Freq, TeoDR[1]-entorno), GetClosestIndex(Freq, TeoDR[1]+entorno)
+    IndiceDRteo3, IndiceEntornoinicialDRteo3, IndiceEntornofinalDRteo3 = GetClosestIndex(Freq, TeoDR[2]), GetClosestIndex(Freq, TeoDR[2]-entorno), GetClosestIndex(Freq, TeoDR[2]+entorno)
+    IndiceDRteo4, IndiceEntornoinicialDRteo4, IndiceEntornofinalDRteo4 = GetClosestIndex(Freq, TeoDR[3]), GetClosestIndex(Freq, TeoDR[3]-entorno), GetClosestIndex(Freq, TeoDR[3]+entorno)
+    IndiceDRteo5, IndiceEntornoinicialDRteo5, IndiceEntornofinalDRteo5 = GetClosestIndex(Freq, TeoDR[4]), GetClosestIndex(Freq, TeoDR[4]-entorno), GetClosestIndex(Freq, TeoDR[4]+entorno)
+    IndiceDRteo6, IndiceEntornoinicialDRteo6, IndiceEntornofinalDRteo6 = GetClosestIndex(Freq, TeoDR[5]), GetClosestIndex(Freq, TeoDR[5]-entorno), GetClosestIndex(Freq, TeoDR[5]+entorno)
+
+    EntornoFreqDR1, EntornoFreqDR2 = Freq[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Freq[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
+    EntornoFreqDR3, EntornoFreqDR4 = Freq[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Freq[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
+    EntornoFreqDR5, EntornoFreqDR6 = Freq[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Freq[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
+
+    EntornoFluoDR1, EntornoFluoDR2 = Fluo[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Fluo[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
+    EntornoFluoDR3, EntornoFluoDR4 = Fluo[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Fluo[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
+    EntornoFluoDR5, EntornoFluoDR6 = Fluo[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Fluo[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
+    
+    IndiceFluoMinimaEntorno1, IndiceFluoMinimaEntorno2 = argrelextrema(np.array(EntornoFluoDR1), np.less)[0], argrelextrema(np.array(EntornoFluoDR2), np.less)[0]
+    IndiceFluoMinimaEntorno3, IndiceFluoMinimaEntorno4 = argrelextrema(np.array(EntornoFluoDR3), np.less)[0], argrelextrema(np.array(EntornoFluoDR4), np.less)[0]
+    IndiceFluoMinimaEntorno5, IndiceFluoMinimaEntorno6 = argrelextrema(np.array(EntornoFluoDR5), np.less)[0], argrelextrema(np.array(EntornoFluoDR6), np.less)[0]
+        
+    try:
+        FreqDR1 = EntornoFreqDR1[int(IndiceFluoMinimaEntorno1)]
+        IndiceDR1 = GetClosestIndex(Freq, FreqDR1)
+    except:
+        FreqDR1 = TeoDR[0]
+        IndiceDR1 = IndiceDRteo1
+    try:
+        FreqDR2 = EntornoFreqDR2[int(IndiceFluoMinimaEntorno2)]
+        IndiceDR2 = GetClosestIndex(Freq, FreqDR2)
+    except:
+        FreqDR2 = TeoDR[1]
+        IndiceDR2 = IndiceDRteo2
+    try:
+        FreqDR3 = EntornoFreqDR3[int(IndiceFluoMinimaEntorno3)]
+        IndiceDR3 = GetClosestIndex(Freq, FreqDR3)
+    except:
+        FreqDR3 = TeoDR[2]
+        IndiceDR3 = IndiceDRteo3
+    try:
+        FreqDR4 = EntornoFreqDR4[int(IndiceFluoMinimaEntorno4)]
+        IndiceDR4 = GetClosestIndex(Freq, FreqDR4)
+    except:
+        FreqDR4 = TeoDR[3]
+        IndiceDR4 = IndiceDRteo4
+    try:
+        FreqDR5 = EntornoFreqDR5[int(IndiceFluoMinimaEntorno5)]
+        IndiceDR5 = GetClosestIndex(Freq, FreqDR5)
+    except:
+        FreqDR5 = TeoDR[4]
+        IndiceDR5 = IndiceDRteo5
+    try:
+        FreqDR6 = EntornoFreqDR6[int(IndiceFluoMinimaEntorno6)]
+        IndiceDR6 = GetClosestIndex(Freq, FreqDR6)
+    except:
+        FreqDR6 = TeoDR[5]
+        IndiceDR6 = IndiceDRteo6
+                
+    return [IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6], [FreqDR1, FreqDR2, FreqDR3, FreqDR4, FreqDR5, FreqDR6]
+
+
+def FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=-100, getindices=False):
+    """
+    Toma los indices donde estan las DR y evalua su fluorescencia. Esos indices son minimos locales en un entorno
+    cercano a las DR teoricas y, si no hay ningun minimo, toma la teorica.
+    Luego, hace el cociente de esa fluorescencia y un factor de normalización segun NormalizationCriterium:
+        1: Devuelve la fluorescencia absoluta de los minimos
+        2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no 
+        hubiera una resonancia oscura y hubiera una recta. Ese valor esta a  DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
+        3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
+        o el valor limite por izquierda de la curva
+        4: Deuelve el cociente entre la fluorescencia del minimo y el valor de fluorescencia a detuning 0 MHz        
+    
+    """
+    IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6 = IndicesDR[0], IndicesDR[1], IndicesDR[2], IndicesDR[3], IndicesDR[4], IndicesDR[5]
+    FluorescenceOfMinimums = [Fluo[IndiceDR1], Fluo[IndiceDR2], Fluo[IndiceDR3], Fluo[IndiceDR4], Fluo[IndiceDR5], Fluo[IndiceDR6]]
+    FrequencyOfMinimums = [Freq[IndiceDR1], Freq[IndiceDR2], Freq[IndiceDR3], Freq[IndiceDR4], Freq[IndiceDR5], Freq[IndiceDR6]]
+    DistanciaFrecuenciaCociente = 25
+
+    if NormalizationCriterium==0:
+        print('che')
+        return FrequencyOfMinimums, FluorescenceOfMinimums
+    
+    if NormalizationCriterium==1:      
+        Fluorescenciacerodetuning = Fluo[GetClosestIndex(Freq, 0)]
+        Fluorescenciaasintotica = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
+        return FrequencyOfMinimums, np.array([Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica])
+        
+    if NormalizationCriterium==2:
+        k = 0
+        while k < len(Freq):
+            if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
+                FluoIzquierda = Fluo[k]
+                indiceizquierda = k
+                print('Izq:', Freq[k])
+                break
+            else:
+                k = k + 1
+        
+        l = 0
+        while l < len(Freq):
+            if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
+                FluoDerecha = Fluo[l]
+                indicederecha = l
+                print('Der: ', Freq[l])
+                break
+            else:
+                l = l + 1
+        FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
+        print(FluoNormDivisor)
+    
+    if NormalizationCriterium==3:
+        #asintotico
+        FluoNormDivisor = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
+
+    if NormalizationCriterium==4:
+        #este te tira la fluorescencia de detuning 0
+        FluoNormDivisor = Fluo[GetClosestIndex(Freq, 0)]
+
+   
+    RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
+    print('Esto: ', RelativeFluorescenceOfMinimums)
+
+    if NormalizationCriterium==2 and getindices==True:
+        return FrequencyOfMinimums, RelativeFluorescenceOfMinimums, indiceizquierda, indicederecha
+    
+    return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
+
+
+def GetFinalMaps(MapasDR1, MapasDR2, MapasDR3, MapasDR4, MapasDR5, MapasDR6):
+    """
+    Nota: esto vale para polarizacion del 397 sigma+ + sigma-. Sino hay que cambiar los coeficientes.
+    La estructura es:
+        MapasDRi = [MapaMedido_criterio1_DRi, MapaMedido_criterio2_DRi, MapaMedido_criterio3_DRi, MapaMedido_criterio4_DRi]
+    """
+    Mapa1 = MapasDR1[0]
+    
+    Mapa2pi = np.sqrt(3)*(MapasDR2[1] + MapasDR5[1])
+    Mapa2smas = np.sqrt(12/2)*MapasDR3[1] + (2/np.sqrt(2))*MapasDR6[1]
+    Mapa2smenos = (2/np.sqrt(2))*MapasDR1[1] + np.sqrt(12/2)*MapasDR4[1]
+    
+    Mapa3pi = np.sqrt(3)*(MapasDR2[2] + MapasDR5[2])
+    Mapa3smas = np.sqrt(12/2)*MapasDR3[2] + (2/np.sqrt(2))*MapasDR6[2]
+    Mapa3smenos = (2/np.sqrt(2))*MapasDR1[2] + np.sqrt(12/2)*MapasDR4[2]
+        
+    return Mapa1, [Mapa2pi, Mapa2smas, Mapa2smenos], [Mapa3pi, Mapa3smas, Mapa3smenos]
+
+
+def CombinateDRwithCG(RelMinMedido1, RelMinMedido2, RelMinMedido3, RelMinMedido4):
+    Fluo1 = RelMinMedido1[0]
+    
+    Fluo2pi = np.sqrt(3)*(RelMinMedido2[1] + RelMinMedido2[4])
+    Fluo2smas = np.sqrt(12/2)*RelMinMedido2[2] + (2/np.sqrt(2))*RelMinMedido2[5]
+    Fluo2smenos = (2/np.sqrt(2))*RelMinMedido2[0] + np.sqrt(12/2)*RelMinMedido2[3]
+    
+    Fluo3pi = np.sqrt(3)*(RelMinMedido3[1] + RelMinMedido3[4])
+    Fluo3smas = np.sqrt(12/2)*RelMinMedido3[2] + (2/np.sqrt(2))*RelMinMedido3[5]
+    Fluo3smenos = (2/np.sqrt(2))*RelMinMedido3[0] + np.sqrt(12/2)*RelMinMedido3[3]
+    
+    return Fluo1, [Fluo2pi, Fluo2smas, Fluo2smenos], [Fluo3pi, Fluo3smas, Fluo3smenos]
+
+
+
+def IdentifyPolarizationCoincidences(theoricalmap, target, tolerance=1e-1):
+    """
+    Busca en un mapa 2D la presencia de un valor target (medido) con tolerancia tolerance.
+    Si lo encuentra, pone un 1. Sino, un 0. Al plotear con pcolor se verá
+    en blanco la zona donde el valor medido se puede hallar.
+    """
+    CoincidenceMatrix = np.zeros((len(theoricalmap), len(theoricalmap[0])))
+    
+    i = 0
+    while i<len(theoricalmap):
+        j = 0
+        while j<len(theoricalmap[0]):
+            if abs(theoricalmap[i][j]-target) < tolerance:
+                CoincidenceMatrix[i][j] = 1
+            j=j+1
+        i=i+1
+    return CoincidenceMatrix
+
+def RetrieveAbsoluteCoincidencesBetweenMaps(MapsVectors):
+    
+    MatrixSum = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
+    AbsoluteCoincidencesMatrix = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
+
+    
+    MatrixMapsVectors = []
+    
+    for i in range(len(MapsVectors)):
+        MatrixMapsVectors.append(np.matrix(MapsVectors[i]))
+
+    for i in range(len(MatrixMapsVectors)):
+        MatrixSum = MatrixSum + MatrixMapsVectors[i]
+        
+    MaxNumberOfCoincidences = np.max(MatrixSum)
+    
+    ListMatrixSum = [list(i) for i in list(np.array(MatrixSum))]
+    
+    for i in range(len(ListMatrixSum)):
+        for j in range(len(ListMatrixSum[0])):
+            if ListMatrixSum[i][j] == MaxNumberOfCoincidences:
+                AbsoluteCoincidencesMatrix[i][j] = 1
+    
+    return AbsoluteCoincidencesMatrix, MaxNumberOfCoincidences
+
+
+def MeasureMeanValueOfEstimatedArea(AbsoluteCoincidencesMap, X, Y):
+    NonZeroIndices = np.nonzero(AbsoluteCoincidencesMap)
+    
+    Xsum = 0
+    Xvec = []
+    Ysum = 0
+    Yvec = []
+    N = len(NonZeroIndices[0])
+    for i in range(N):
+        Xsum = Xsum + X[NonZeroIndices[1][i]]
+        Xvec.append(X[NonZeroIndices[1][i]])
+        Ysum = Ysum + Y[NonZeroIndices[0][i]]
+        Yvec.append(Y[NonZeroIndices[0][i]])
+    Xaverage = Xsum/N
+    Yaverage = Ysum/N
+    Xspread = np.std(Xvec)
+    Yspread = np.std(Yvec)
+    return Xaverage, Yaverage, N, Xspread, Yspread
+    
+
+def MeasureRelativeFluorescenceFromCPT(Freq, Fluo, u, titadoppler, detuningrepump, detuningdoppler, frefasint=-100, entorno=3):
+    ResonanciasTeoricas, ResonanciasPositivas = CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)
+    IndicesDR, FreqsDR = FindDRFrequencies(Freq, Fluo, ResonanciasTeoricas, entorno=entorno)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums0 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=0, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums1 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums2, indiceizquierda, indicederecha = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=2, frecuenciareferenciacriterioasintotico=frefasint, getindices=True)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums3 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=3, frecuenciareferenciacriterioasintotico=frefasint)
+    FrequencyOfMinimums, RelativeFluorescenceOfMinimums4 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=4, frecuenciareferenciacriterioasintotico=frefasint)
+    print('hola')
+    print(RelativeFluorescenceOfMinimums0)
+    return RelativeFluorescenceOfMinimums0, RelativeFluorescenceOfMinimums1, RelativeFluorescenceOfMinimums2, RelativeFluorescenceOfMinimums3, RelativeFluorescenceOfMinimums4, IndicesDR, [indiceizquierda, indicederecha]
+
+
+def GenerateNoisyCPT(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+
+def GenerateNoisyCPT_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+def GenerateNoisyCPT_fit(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqs, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, min(freqs), max(freqs) + freqs[1]-freqs[0], freqs[1]-freqs[0], plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+
+def AddNoiseToCPT(Fluo, noisefactor):
+    return [f+noisefactor*(2*random.random()-1) for f in Fluo]
+
+
+def SmoothNoisyCPT(Fluo, window=11, poly=3):
+    SmoothenFluo = sf(Fluo, window, poly)
+    return SmoothenFluo
+
+
+def GetMinimaInfo(Freq, Fluo, u, titadoppler, detuningdoppler, detuningrepump, MinimumCriterium=2, NormalizationCriterium=1):
+    """
+    FUNCION VIEJA
+    
+    Esta funcion devuelve valores de frecuencias y fluorescencia relativa de los minimos.
+    Minimumcriterion:
+        1: Saca los minimos con funcion argelextrema
+        2: Directamente con las frecuencias teoricas busca las fluorescencias
+    
+    Normalizationcriterium:
+        1: Devuelve la fluorescencia absoluta de los minimos
+        2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no 
+        hubiera una resonancia oscura y hubiera una recta. Ese valor esta a  DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
+        3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
+        o el valor limite por izquierda de la curva
+    
+    """
+    FluorescenceOfMaximum = max(Fluo)
+    FrequencyOfMaximum = Freq[Fluo.index(FluorescenceOfMaximum)]
+
+
+    #criterio para encontrar los minimos
+    #criterio usando minimos de la fluorescencia calculados con la curva
+    if MinimumCriterium == 1:
+        LocationOfMinimums = argrelextrema(np.array(Fluo), np.less)[0]
+        FluorescenceOfMinimums = np.array([Fluo[i] for i in LocationOfMinimums])
+        FrequencyOfMinimums = np.array([Freq[j] for j in LocationOfMinimums])
+
+
+    #criterio con las DR teoricas
+    if MinimumCriterium == 2:
+        FrecuenciasDRTeoricas, FrecuenciasDRTeoricasPositivas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
+        FrequencyOfMinimums = []
+        FluorescenceOfMinimums =[]
+        print(FrecuenciasDRTeoricas)
+        k=0
+        
+        ventanita = 0.001
+        while k < len(Freq):
+            if Freq[k] < FrecuenciasDRTeoricas[0] + ventanita and Freq[k] > FrecuenciasDRTeoricas[0] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[1] + ventanita and Freq[k] > FrecuenciasDRTeoricas[1] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[2] + ventanita and Freq[k] > FrecuenciasDRTeoricas[2] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[3] + ventanita and Freq[k] > FrecuenciasDRTeoricas[3] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[4] + ventanita and Freq[k] > FrecuenciasDRTeoricas[4] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            elif Freq[k] < FrecuenciasDRTeoricas[5] + ventanita and Freq[k] > FrecuenciasDRTeoricas[5] - ventanita:
+                FrequencyOfMinimums.append(Freq[k])
+                FluorescenceOfMinimums.append(Fluo[k])
+            k = k + 1
+        print(FrequencyOfMinimums)
+        if len(FrequencyOfMinimums) !=  len(FrecuenciasDRTeoricas):
+            print('NO ANDA BIEN ESTO PAPI, revisalo')
+        
+            
+    #esto es para establecer un criterio para la fluorescencia relativa
+    DistanciaFrecuenciaCociente = 15
+    
+    if NormalizationCriterium==1:      
+        FluoNormDivisor = 1 
+        
+    if NormalizationCriterium==2:
+        k = 0
+        while k < len(Freq):
+            if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
+                FluoIzquierda = Fluo[k]
+                print('Izq:', Freq[k])
+                break
+            else:
+                k = k + 1
+        
+        l = 0
+        while l < len(Freq):
+            if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
+                FluoDerecha = Fluo[l]
+                print('Der: ', Freq[l])
+                break
+            else:
+                l = l + 1
+        FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
+        print(FluoNormDivisor)
+    
+    if NormalizationCriterium==3:
+        FluoNormDivisor = Fluo[0]
+    
+    RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
+    
+    return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
+
+
+def GetPlotsofFluovsAngle_8levels(FrequencyOfMinimumsVector, RelativeFluorescenceOfMinimumsVector, u, titadoppler, detuningdoppler, detuningrepump, ventana=0.25, taketheoricalDR=False):
+    
+    #primero buscamos las frecuencias referencia que se parezcan a las 6:
+    i = 0
+    FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[0]
+    
+    FrecuenciasDRTeoricas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
+    
+    while i < len(FrequencyOfMinimumsVector):
+        if len(FrequencyOfMinimumsVector[i])==len(FrecuenciasDRTeoricas):
+            FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[i]
+            print('Cool! Taking the DR identified with any curve')
+            break
+        else:
+            i = i + 1
+            if i==len(FrequencyOfMinimumsVector):
+                print('No hay ningun plot con 5 resonancias oscuras. Tomo las teóricas')
+                FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
+    
+    if taketheoricalDR:
+        FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
+        
+        
+            
+    Ventana = abs(ventana*(FrecuenciasReferenciaBase[1] - FrecuenciasReferenciaBase[0])) #ventana separadora de resonancias
+    
+    print('Ventana = ', Ventana)
+    DarkResonance1Frequency = []
+    DarkResonance1Fluorescence = []
+    
+    DarkResonance2Frequency = []
+    DarkResonance2Fluorescence = []
+    
+    DarkResonance3Frequency = []
+    DarkResonance3Fluorescence = []
+   
+    DarkResonance4Frequency = []
+    DarkResonance4Fluorescence = []
+    
+    DarkResonance5Frequency = []
+    DarkResonance5Fluorescence = []
+    
+    DarkResonance6Frequency = []
+    DarkResonance6Fluorescence = []    
+    
+    i = 0
+    
+   
+    while i < len(FrequencyOfMinimumsVector):
+        j = 0
+        FrecuenciasReferencia = [i for i in FrecuenciasReferenciaBase]
+
+
+        while j < len(FrequencyOfMinimumsVector[i]):
+            
+            
+            if abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[0])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[0])-Ventana):
+                DarkResonance1Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance1Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[0] = 0
+            
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[1])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[1])-Ventana):
+                DarkResonance2Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance2Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[1] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[2])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[2])-Ventana):
+                DarkResonance3Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance3Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[2] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[3])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[3])-Ventana):
+                DarkResonance4Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance4Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[3] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[4])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[4])-Ventana):
+                DarkResonance5Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance5Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[4] = 0
+
+            elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[5])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[5])-Ventana):
+                DarkResonance6Frequency.append(FrequencyOfMinimumsVector[i][j])
+                DarkResonance6Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
+                FrecuenciasReferencia[5] = 0
+                
+            else:
+                #print('Algo anduvo mal, por ahi tenes que cambiar la ventana che')
+                pass
+            j = j + 1
+        
+        if np.count_nonzero(FrecuenciasReferencia) > 0:
+            if FrecuenciasReferencia[0] != 0:
+                DarkResonance1Frequency.append(FrecuenciasReferencia[0])
+                DarkResonance1Fluorescence.append()
+                
+            if FrecuenciasReferencia[1] != 0:
+                DarkResonance2Frequency.append(FrecuenciasReferencia[1])
+                DarkResonance2Fluorescence.append(0)
+
+            if FrecuenciasReferencia[2] != 0:
+                DarkResonance3Frequency.append(FrecuenciasReferencia[2])
+                DarkResonance3Fluorescence.append(0)
+
+            if FrecuenciasReferencia[3] != 0:
+                DarkResonance4Frequency.append(FrecuenciasReferencia[3])
+                DarkResonance4Fluorescence.append(0)
+
+            if FrecuenciasReferencia[4] != 0:
+                DarkResonance5Frequency.append(FrecuenciasReferencia[4])
+                DarkResonance5Fluorescence.append(0)
+
+            if FrecuenciasReferencia[5] != 0:
+                DarkResonance6Frequency.append(FrecuenciasReferencia[5])
+                DarkResonance6Fluorescence.append(0)            
+            
+        i = i + 1
+
+    return DarkResonance1Frequency, DarkResonance1Fluorescence, DarkResonance2Frequency, DarkResonance2Fluorescence, DarkResonance3Frequency, DarkResonance3Fluorescence, DarkResonance4Frequency, DarkResonance4Fluorescence, DarkResonance5Frequency, DarkResonance5Fluorescence, DarkResonance6Frequency, DarkResonance6Fluorescence, FrecuenciasReferenciaBase
+
+
+
+def PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
+    
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+    
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    
+    Fluovectors = []
+    
+    for titaprobe in titaprobeVec:
+        
+        tinicial = time.time()
+
+        ProbeDetuningVectorL, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
+
+        tfinal = time.time()
+        
+        print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")   
+        
+        if plot:
+            plt.figure()
+            plt.xlabel('Repump detuning (MHz')
+            plt.ylabel('Fluorescence (A.U.)')
+            plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
+            plt.legend()
+            
+        Fluovectors.append(Fluovector)
+        
+        if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
+            Fluovectors = Fluovector
+   
+    return ProbeDetuningVectorL, Fluovectors
+
+
+
+def PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
+    
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+    
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    
+    Fluovectors = []
+    
+    for titaprobe in titaprobeVec:
+        
+        tinicial = time.time()
+
+        ProbeDetuningVectorL, Fluovector = CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
+
+        tfinal = time.time()
+        
+        print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")   
+        
+        if plot:
+            plt.figure()
+            plt.xlabel('Repump detuning (MHz')
+            plt.ylabel('Fluorescence (A.U.)')
+            plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
+            plt.legend()
+            
+        Fluovectors.append(Fluovector)
+        
+        if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
+            Fluovectors = Fluovector
+   
+    return ProbeDetuningVectorL, Fluovectors
+
+
+

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/threeLevel_2repumps_CPTPlotter.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Sep  1 17:58:39 2020
+
+@author: oem
+"""
+
+
+
+import os
+import numpy as np
+#os.chdir('/home/oem/Nextcloud/G_liaf/liaf-TrampaAnular/Código General/EIT-CPT/Buenos Aires/Experiment Simulations/CPT scripts/Eight Level 2 repumps')
+from threeLevel_2repumps_AnalysisFunctions import CalculoTeoricoDarkResonances_8levels, GetMinimaInfo, GetPlotsofFluovsAngle_8levels, PerformExperiment_8levels, FindDRFrequencies, FindRelativeFluorescencesOfDR, GenerateNoisyCPT, SmoothNoisyCPT, GetFinalMaps, GenerateNoisyCPT_fixedRabi, GenerateNoisyCPT_fit
+import matplotlib.pyplot as plt
+import time
+from threeLevel_2repumps_AnalysisFunctions import MeasureRelativeFluorescenceFromCPT, IdentifyPolarizationCoincidences, RetrieveAbsoluteCoincidencesBetweenMaps, GetClosestIndex
+
+#C:\Users\Usuario\Nextcloud\G_liaf\liaf-TrampaAnular\Código General\EIT-CPT\Buenos Aires\Experiment Simulations\CPT scripts\Eight Level 2 repumps
+
+ub = 9.27e-24
+h = 6.63e-34
+c = (ub/h)*1e-4 #en unidades de MHz/G
+
+#u = 1e6
+u = 33.5e6
+
+B = (u/(2*np.pi))/c
+
+
+#sg, sp = 0.6, 5  #parámetros de control, saturación del doppler y repump
+#rabG, rabP = sg*gPS, sp*gPD #frecuencias de rabi
+
+gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+
+
+
+lw = 0.1
+DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth = lw, lw, lw #ancho de linea de los laseres
+
+DetDoppler = -36 #42
+DetRepumpVec = [DetDoppler+29.6]
+
+Tvec = [0.7] #temperatura en mK
+alpha = 0*(np.pi/180) #angulo entre los láseres
+
+
+phidoppler, titadoppler = 0, 90
+phirepump, titarepump = 0,  0
+phiprobe = 0
+titaprobe = 90
+
+#Calculo las resonancias oscuras teóricas
+#ResonanciasTeoricas, DRPositivas = CalculoTeoricoDarkResonances_8levels(u/(2*np.pi*1e6), titadoppler, DetDoppler, DetRepump)
+
+#Parametros de la simulacion cpt
+center = -45
+span = 80
+freqMin = center-span*0.5
+freqMax = center+span*0.5
+
+""" parametros para tener espectros coherentes
+freqMin = -56
+freqMax = 14
+"""
+freqStep = 1e-1
+
+noiseamplitude = 0
+
+RelMinMedido0Vector = []
+RelMinMedido1Vector = []
+RelMinMedido2Vector = []
+RelMinMedido3Vector = []
+RelMinMedido4Vector = []
+
+#Sr = np.arange(0, 10, 0.2)
+
+#Sg = np.arange(0.01, 1, 0.05)
+#Sp = np.arange(0.1, 6.1, 1)
+
+#Sg = [0.6**2]
+#Sp = [2.3**2]
+
+
+Sg = [1.4]
+Sp = [6]
+
+Sr = [11]
+
+
+
+
+i = 0
+
+save = False
+
+
+showFigures = True
+
+if not showFigures:
+    plt.ioff()
+else:
+    plt.ion()
+    
+fig1, ax1 = plt.subplots()
+
+offsetx = 464
+
+ax1.plot([f-offsetx for f in FreqsDR], CountsDR, 'o')
+
+run = True
+
+Scale = 730
+Offset = 600 #600 para 20k cuentas aprox
+
+MaxCoherenceValue = []
+
+for sg in Sg:
+    for sp in Sp:
+        rabG, rabP = sg*gPS, sp*gPD
+        for Ti in Tvec:
+            T = Ti*1e-3
+            for DetRepump in DetRepumpVec:
+                print(T)
+                for sr in Sr:
+                
+                    rabR = sr*gPD
+                    
+                    #MeasuredFreq, MeasuredFluo = GenerateNoisyCPT(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, [titaprobe], phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=noiseamplitude)
+                    if run:
+                        MeasuredFreq4, MeasuredFluo4 = GenerateNoisyCPT_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, [titaprobe], phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=noiseamplitude)
+                    
+                    #SmoothFluo = SmoothNoisyCPT(MeasuredFluo, window=9, poly=2)
+                    SmoothFluo4 = MeasuredFluo4
+                    
+                    #Scale = max(BestC)/max([100*s for s in SmoothFluo4])
+                    ax1.plot(MeasuredFreq4, [Scale*100*f + Offset for f in SmoothFluo4], label=f'Sr = {sr}')
+                    ax1.axvline(DetDoppler, linestyle='--', linewidth=1)
+                    #if sr != 0:
+                        #ax1.axvline(DetRepump, linestyle='--', linewidth=1)
+
+                    MaxCoherenceValue.append(np.max(SmoothFluo4))
+
+                    #print(titaprobe)
+                    ax1.set_xlabel('Detuning Rebombeo (MHz)')
+                    ax1.set_ylabel('Fluorescencia (AU)')
+                    ax1.set_title(f'B: {round(B, 2)} G, Sdop: {round(sg, 2)}, Sp: {round(sp, 2)}, Sr: {round(sr, 2)}, lw: {lw} MHz, T: {Ti} mK')
+                    #ax1.set_ylim(0, 8)
+                    #ax1.axvline(DetDoppler, linestyle='dashed', color='red', linewidth=1)
+                    #ax1.axvline(DetRepump, linestyle='dashed', color='black', linewidth=1)
+                    #ax1.set_title('Pol Doppler y Repump: Sigma+ Sigma-, Pol Probe: PI')
+                    #ax1.legend()
+                    ax1.grid()
+                    print (f'{i+1}/{len(Sg)*len(Sp)}')
+                    i = i + 1
+                    if save:
+                        plt.savefig(f'Mapa_plots_100k_1mk/CPT_SMSM_sdop{round(sg, 2)}_sp{round(sp, 2)}_sr{round(sr, 2)}.jpg')
+                    ax1.legend()
+                    
+"""
+plt.figure()
+plt.plot(Sr, MaxCoherenceValue, 'o')
+plt.xlabel('Sr')
+plt.ylabel('Coherence')
+"""
+
+"""
+                        
+plt.figure()
+plt.plot(MeasuredFreq, [100*f for f in SmoothFluo], color='darkred')
+plt.xlabel('Desintonía 866 (MHz)')
+plt.ylabel('Fluorescencia (A.U.)')
+plt.axvline(-30, color='darkblue', linewidth=1.2, linestyle='--')
+plt.yticks(np.arange(0.4, 1.8, 0.2))
+plt.ylim(0.5, 1.6)
+plt.grid()
+              
+
+
+plt.figure()
+plt.plot(MeasuredFreq4, [100*f for f in SmoothFluo4], color='darkred')
+plt.xlabel('Desintonía 866 (MHz)')
+plt.ylabel('Fluorescencia (A.U.)')
+plt.axvline(-30, color='darkblue', linewidth=1.2, linestyle='--')
+plt.yticks(np.arange(0.8, 2.4, 0.4))
+plt.grid()
+              
+"""
+
+
+
+
+#%%
+from scipy.optimize import curve_fit
+
+T = 0.5e-3
+
+sg = 0.7
+sp = 6
+sr = 0
+DetDoppler = -14
+DetRepump = 0
+
+FitsSp = []
+FitsOffset = []
+
+Sg = [0.87]
+
+def FitEIT(freqs, SP, offset):
+    MeasuredFreq, MeasuredFluo = GenerateNoisyCPT_fit(0.87, sr, SP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, [titaprobe], phirepump, titarepump, freqs, plot=False, solvemode=1, detpvec=None, noiseamplitude=noiseamplitude)
+    FinalFluo = [f*43000 + 2685 for f in MeasuredFluo]
+    return FinalFluo
+
+freqs = [f-offsetx+32 for f in FreqsDR]
+freqslong = np.arange(min(freqs), max(freqs)+freqs[1]-freqs[0], 0.1*(freqs[1]-freqs[0]))
+
+popt, pcov = curve_fit(FitEIT, freqs, CountsDR, p0=[5, 700], bounds=(0, [10, 1e6]))
+
+FitsSp.append(popt[0])
+FitsOffset.append(popt[1])
+
+print(popt)
+
+FittedEIT = FitEIT(freqslong, *popt)
+
+
+plt.figure()
+plt.errorbar(freqs, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', capsize=2, markersize=2)
+plt.plot(freqslong, FitEIT(freqslong, *popt))
+plt.title(f'Sdop: {round(popt[0], 2)}, Spr: {round(popt[1], 2)}, T: {T*1e3} mK, detDop: {DetDoppler} MHz')
+
+np.savetxt('CPT_measured.txt', np.transpose([freqs, CountsDR]))
+np.savetxt('CPT_fitted.txt', np.transpose([freqslong, FittedEIT]))

analisis/plots/20250516_CPT_bladetrap/Data2/EITfit/threeLevel_2repumps_linealpol_python_scripts.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr  7 22:30:01 2020
+
+@author: nico
+"""
+
+
+import numpy as np
+import time
+import matplotlib.pyplot as plt
+from scipy.signal import argrelextrema
+
+
+"""
+Scripts para el calculo de la curva CPT
+"""
+
+
+def H0matrix(Detg, Detp, u): 
+    """
+    Calcula la matriz H0 en donde dr es el detuning del doppler, dp es el retuning del repump y u es el campo magnético en Hz/Gauss.
+    Para esto se toma la energía del nivel P como 0
+    """
+    
+    eigenEnergies = (Detg-u, Detg+u, -u/3, u/3, Detp-6*u/5, Detp-2*u/5, Detp+2*u/5, Detp+6*u/5) #pagina 26 de Oberst. los lande del calcio son iguales a Bario.
+    H0 = np.diag(eigenEnergies)
+    
+    return H0
+
+
+def HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe):
+    """
+    Calcula la matriz de interacción Hsp + Hpd, en donde rabR es la frecuencia de rabi de la transición Doppler SP, 
+    rabP es la frecuencia de rabi de la transición repump DP, y las componentes ei_r y ei_p son las componentes de la polarización
+    del campo eléctrico incidente de doppler y repump respectivamente. Deben estar normalizadas a 1
+    """
+    HI = np.zeros((8, 8), dtype=np.complex_)
+    
+    i, j = 1, 3
+    HI[i-1, j-1] = (rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1]) 
+    
+    i, j = 1, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 2, 3
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(-1j*phidoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+     
+    i, j = 2, 4
+    HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.cos(titadoppler)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])   
+    
+    
+    i, j = 3, 5
+    HI[i-1, j-1] = -(rabP/2) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 3, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 3, 7
+    HI[i-1, j-1] = rabP/np.sqrt(12) * np.sin(titaprobe)*np.exp(1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    
+    i, j = 4, 6
+    HI[i-1, j-1] = -(rabP/np.sqrt(12)) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 4, 7
+    HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe) 
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    i, j = 4, 8
+    HI[i-1, j-1] = (rabP/2) * np.sin(titaprobe)*np.exp(1j*phiprobe)
+    HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
+    
+    return HI
+
+
+def Lplusminus(detr, detp, phirepump, titarepump, forma=1):
+    Hintplus = np.zeros((8, 8), dtype=np.complex_)
+    Hintminus = np.zeros((8, 8), dtype=np.complex_)
+
+
+    Hintplus[4, 2] = (-1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
+    Hintplus[5, 2] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintplus[6, 2] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
+    
+    Hintplus[5, 3] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
+    Hintplus[6, 3] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintplus[7, 3] = (1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
+    
+
+    Hintminus[2, 4] = (-1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
+    Hintminus[2, 5] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintminus[2, 6] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
+    
+    Hintminus[3, 5] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
+    Hintminus[3, 6] = (-1/np.sqrt(3))*np.cos(titarepump)
+    Hintminus[3, 7] = (1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
+
+    if forma==1:
+        Lplus = np.zeros((64, 64), dtype=np.complex_)
+        Lminus = np.zeros((64, 64), dtype=np.complex_)
+        DeltaBar = np.zeros((64, 64), dtype=np.complex_)
+    
+        for r  in range(8):
+            for q in range(8):
+                for k in range(8):
+                    for j in range(8):
+                        if j==q:
+                            if (k==2 or k==3) and r > 3:
+                                Lplus[r*8+q][k*8+j] = (-1j)*(Hintplus[r,k]) 
+                            if (r==2 or r==3) and k > 3:
+                                Lminus[r*8+q][k*8+j] = (-1j)*(Hintminus[r,k])
+                        elif r==k:
+                            if (q==2 or q==3) and j > 3:
+                                Lplus[r*8+q][k*8+j] = (-1j)*(- Hintplus[j,q]) 
+                            if (j==2 or j==3) and q > 3:
+                                Lminus[r*8+q][k*8+j] = (-1j)*(- Hintminus[j,q])
+    if forma==2:
+        deltaKro = np.diag([1, 1, 1, 1, 1, 1, 1, 1])
+        Lplus = (-1j)*(np.kron(Hintplus, deltaKro) - np.kron(deltaKro, Hintplus))
+        Lminus = (-1j)*(np.kron(Hintminus, deltaKro) - np.kron(deltaKro, Hintminus))
+    
+    DeltaBar = np.zeros((64, 64), dtype=np.complex_)   
+    for i in range(64):
+        
+        DeltaBar[i, i] = (1j)*(detr - detp)
+
+    return np.matrix(Lminus), np.matrix(Lplus), np.matrix(DeltaBar)
+
+
+def GetL1(Lplus, Lminus, DeltaBar, L0, rabR, nmax):
+    """
+    Devuelve Splus0 y Sminus0
+    """
+    
+    Sp = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 - (nmax+1)*DeltaBar))*np.matrix(Lplus))
+    Sm = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 + (nmax+1)*DeltaBar))*np.matrix(Lminus))
+    
+    for n in list(range(nmax+1))[(nmax+1)::-1][0:len(list(range(nmax+1))[(nmax+1)::-1])-1]: #jaja esto solo es para que vaya de nmax a 1 bajando. debe haber algo mas facil pero kcio
+        Sp = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 - n*DeltaBar + rabR*(Lminus*np.matrix(Sp))))*np.matrix(Lplus))
+        Sm = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 + n*DeltaBar + rabR*(Lplus*np.matrix(Sm))))*np.matrix(Lminus))
+    L1 = 0.5*rabR*(np.matrix(Lminus)*np.matrix(Sp) + np.matrix(Lplus)*np.matrix(Sm))
+    
+    return L1
+
+
+def EffectiveL(gPS, gPD, lwg, lwr, lwp):
+
+    """
+    Siendo Heff = H + EffectiveL, calcula dicho EffectiveL que es (-0.5j)*sumatoria(CmDaga*Cm) que luego sirve para calcular el Liouvilliano
+    """
+    Leff = np.zeros((8, 8), dtype=np.complex_)
+    
+    Leff[0, 0] = 2*lwg
+    Leff[1, 1] = 2*lwg
+    Leff[2, 2] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[3, 3] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
+    Leff[4, 4] = 2*(lwr + lwp)
+    Leff[5, 5] = 2*(lwr + lwp)
+    Leff[6, 6] = 2*(lwr + lwp)
+    Leff[7, 7] = 2*(lwr + lwp)
+    
+    return (-0.5j)*Leff
+
+
+def CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp):
+    """
+    Si tomamos el Liuvilliano como L = (-j)*(Heff*deltak - Heffdaga*deltak) + sum(Mm),
+    esta funcion calcula dichos Mm, que tienen dimensión 64x64 ya que esa es la dimensión del L. Estas componentes
+    salen de hacer la cuenta a mano conociendo los Cm y considerando que Mm[8*(r-1)+s, 8*(k-1)+j] = Cm[r,l] + Cmdaga[j,s] = Cm[r,l] + Cm[s,j]
+    ya que los componentes de Cm son reales.
+    Esta M es la suma de las 8 matrices M.
+    
+    """
+    
+    M = np.matrix(np.zeros((64, 64), dtype=np.complex_))
+    
+    M[0,27] = (2/3)*gPS
+    
+    M[9,18] = (2/3)*gPS
+    
+    M[0,18] = (1/3)*gPS
+    M[1,19] = -(1/3)*gPS
+    M[8,26] = -(1/3)*gPS
+    M[9,27] = (1/3)*gPS
+    
+    M[36,18] = (1/2)*gPD
+    M[37,19] = (1/np.sqrt(12))*gPD
+    M[44,26] = (1/np.sqrt(12))*gPD
+    M[45,27] = (1/6)*gPD
+    
+    M[54,18] = (1/6)*gPD
+    M[55,19] = (1/np.sqrt(12))*gPD
+    M[62,26] = (1/np.sqrt(12))*gPD
+    M[63,27] = (1/2)*gPD
+    
+    M[45,18] = (1/3)*gPD
+    M[46,19] = (1/3)*gPD
+    M[53,26] = (1/3)*gPD
+    M[54,27] = (1/3)*gPD
+    
+    M[0,0] = 2*lwg
+    M[1,1] = 2*lwg
+    M[8,8] = 2*lwg
+    M[9,9] = 2*lwg
+    
+    factor1 = 1
+    factor2 = 1
+    factor3 = 1
+    factor4 = 1
+    
+    
+    #M[36, 45] = lwp
+    
+    
+    M[36,36] = 2*(lwr + factor1*lwp)
+    M[37,37] = 2*(lwr + factor1*lwp)
+    M[38,38] = 2*(lwr + factor1*lwp)
+    M[39,39] = 2*(lwr + factor1*lwp)
+    M[44,44] = 2*(lwr + factor2*lwp)
+    M[45,45] = 2*(lwr + factor2*lwp)
+    M[46,46] = 2*(lwr + factor2*lwp)
+    M[47,47] = 2*(lwr + factor2*lwp)
+    M[52,52] = 2*(lwr + factor3*lwp)
+    M[53,53] = 2*(lwr + factor3*lwp)
+    M[54,54] = 2*(lwr + factor3*lwp)
+    M[55,55] = 2*(lwr + factor3*lwp)
+    M[60,60] = 2*(lwr + factor4*lwp)
+    M[61,61] = 2*(lwr + factor4*lwp)
+    M[62,62] = 2*(lwr + factor4*lwp)
+    M[63,63] = 2*(lwr + factor4*lwp)
+
+    return M
+
+def dopplerBroadening(wlg, wlp, alpha, T, mcalcio = 6.655e-23*1e-3):
+    """
+    Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
+    wlg es la longitud de onda doppler, wlp la longitud de onda repump, T la temperatura del ion en kelvin, y alpha (en rads) el ángulo
+    que forman ambos láseres.
+    """
+    
+    kboltzmann = 1.38e-23 #J/K
+    
+    gammaD = (2*np.pi)*np.sqrt((1/(wlg*wlg)) + (1/(wlp*wlp)) - 2*(1/(wlg*wlp))*np.cos(alpha))*np.sqrt(kboltzmann*T/(2*mcalcio))
+
+    return gammaD
+    
+
+def FullL_efficient(rabG, rabR, rabP, gPS = 0, gPD = 0, Detg = 0, Detr = 0, Detp = 0, u = 0, lwg = 0, lwr=0, lwp = 0, 
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, phirepump=0, titarepump=0, T = 0, alpha = 0):
+
+    """
+    Calcula el Liouvilliano total de manera explícita índice a índice. Suma aparte las componentes de las matrices M.
+    Es la más eficiente hasta ahora.
+    """
+    
+    
+    db = dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)
+    
+    #lwr = np.sqrt(lwr**2 + dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)**2)
+    lwg = np.sqrt(lwg**2 + db**2)
+    lwr = np.sqrt(lwr**2 + db**2)    
+    
+    CC = EffectiveL(gPS, gPD, lwg, lwr, lwp)
+    
+    Heff = H0matrix(Detg, Detp, u) + HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe) + CC
+    
+    Heffdaga = np.matrix(Heff).getH()
+    
+    Lfullpartial = np.zeros((64, 64), dtype=np.complex_)
+    
+    for r  in range(8):
+        for q in range(8):
+            for k in range(8):
+                for j in range(8):
+                    if j!=q and r!=k:
+                        pass
+                    elif j==q and r!=k:
+                        if (r < 2 and k > 3) or (k < 2 and r > 3) or (r > 3 and k > 3) or (r==0 and k==1) or (r==1 and k==0) or (r==2 and k==3) or (r==3 and k==2): #todo esto sale de analizar explicitamente la matriz y tratar de no calcular cosas de más que dan cero
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k])
+                    elif j!=q and r==k:
+                        if (j < 2 and q > 3) or (q < 2 and j > 3) or (j > 3 and q > 3) or (j==0 and q==1) or (j==1 and q==0) or (j==2 and q==3) or (j==3 and q==2):
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(-Heffdaga[j,q])
+                    else:
+                        if Heff[r,k] == Heffdaga[j,q]:
+                            pass
+                        else:
+                            Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k]-Heffdaga[j,q])
+
+    M = CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp) 
+    L0 = np.array(np.matrix(Lfullpartial) + M)
+    
+    nmax = 1
+    Lminus, Lplus, DeltaBar = Lplusminus(Detr, Detp, phirepump, titarepump)
+    factor1 = np.exp(1j*0.2*np.pi)
+    factor2 = np.exp(-1j*0.2*np.pi)
+    #print(factor)
+    L1 = GetL1(factor1*Lplus, factor2*Lminus, DeltaBar, L0, rabR, nmax)
+    Lfull = L0 + L1    
+    
+    
+    #NORMALIZACION DE RHO
+    i = 0 
+    while i < 64:
+        if i%9 == 0:
+            Lfull[0, i] = 1
+        else:
+            Lfull[0, i] = 0
+        i = i + 1    
+    
+    return Lfull 
+
+
+"""
+Scripts para correr un experimento y hacer el análisis de los datos
+"""
+
+def CalculoTeoricoDarkResonances(u, titadoppler):
+    if titadoppler==0:
+        NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
+    elif titadoppler==90:
+        NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
+    
+    PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
+    
+    return NegativeDR, PositiveDR
+
+
+def CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, 
+                      freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+    
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180) 
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
+    
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
+    Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
+    lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+    
+    #u = 2*np.pi*u*1e6
+    
+    Fluovector = []
+    
+    tinicial = time.time()
+    
+    for Detp in DetProbeVector:
+        L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
+        
+        if solvemode == 1:
+            rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)]))
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)                
+
+        if solvemode == 2:
+            Linv = np.linalg.inv(L) 
+            rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)
+        
+    tfinal = time.time()
+        
+    print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
+    
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+    
+    if plot:
+        plt.xlabel('Probe detuning (MHz)')
+        plt.ylabel('Fluorescence (A.U.)')
+        plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+        plt.legend()
+        
+    return DetProbeVectorMHz, Fluovector
+
+
+def CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, 
+                      freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
+    """
+    Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
+    
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+
+
+
+    
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180) 
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
+    
+    DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
+    Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
+    #lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
+    lwg, lwr, lwp = lwg*1e6, lwr*1e6, lwp*1e6
+
+    rabG = sg*gPS
+    rabR = sr*gPD
+    rabP = sp*gPD
+    #u = 2*np.pi*u*1e6
+    
+    Fluovector = []
+    
+    tinicial = time.time()
+    
+    for Detp in DetProbeVector:
+
+        L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
+        
+        if solvemode == 1:
+            coh = 5
+            rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)]))
+            #Fluo = np.abs(rhovectorized[coh])
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)                
+
+        if solvemode == 2:
+            Linv = np.linalg.inv(L) 
+            rhovectorized = [Linv[j][0] for j in range(len(Linv))]
+            Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
+            Fluovector.append(Fluo)
+        
+    tfinal = time.time()
+        
+    print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
+    
+    DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
+    
+    if plot:
+        plt.xlabel('Probe detuning (MHz)')
+        plt.ylabel('Fluorescence (A.U.)')
+        plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+        plt.legend()
+        
+    return DetProbeVectorMHz, Fluovector
+#%%
+
+if __name__ == "__main__":
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+    
+    B = 25 #campo magnetico en gauss
+    
+    u = c*B
+    
+    
+    sg, sr, sp = 0.5, 1.5, 4 #parámetros de saturación del doppler y repump
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+    rabG, rabR, rabP = sg*gPS, sr*gPD, sp*gPD #frecuencias de rabi
+    
+    lwg, lwr, lwp = 0.3, 0.3, 0.3 #ancho de linea de los laseres
+    
+    Detg = -25
+    Detr = 20 #detuning del doppler y repump
+    
+    Temp = 0.0e-3 #temperatura en K
+    alpha = 0*(np.pi/180) #angulo entre los láseres
+    
+    phidoppler, titadoppler = 0, 90
+    phirepump, titarepump = 0, 90
+    
+    phiprobe, titaprobe = 0, 90
+    
+    
+    plotCPT = False
+    
+    freqMin = -50
+    freqMax = 50
+    freqStep = 5e-2
+    
+    Frequencyvector, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=plotCPT, solvemode=1)
+    
+    NegativeDR, PositiveDR = CalculoTeoricoDarkResonances(u/(2*np.pi*1e6), titadoppler)
+
+    plt.plot(Frequencyvector, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
+    plt.xlabel('Probe detuning (MHz)')
+    plt.ylabel('Fluorescence (A.U.)')
+    
+    for PDR in PositiveDR:
+        plt.axvline(Detr+PDR, linestyle='--', linewidth=0.5, color='red')
+    for NDR in NegativeDR:
+        plt.axvline(Detg+NDR, linestyle='--', linewidth=0.5, color='blue')
+
+
+
+    
+    #parametros que andan piola:
+        
+    
+    """
+    ub = 9.27e-24
+    h = 6.63e-34
+    c = (ub/h)*1e-4 #en unidades de MHz/G
+    
+    B = 17 #campo magnetico en gauss
+    
+    u = c*B
+    #u = 80e6
+    
+    sr, sp = 0.53, 4.2
+    gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
+    rabR, rabP = sr*gPS, sp*gPD
+    lw = 2*np.pi * 0.33e6
+    lwr, lwp = lw, lw #ancho de linea de los laseres
+    dr_spec = - 2*np.pi* 26e6
+    
+    
+    freqSteps = 500
+    freqMin = -100e6
+    freqMax = 100e6
+    
+    dps = 2*np.pi*np.linspace(freqMin, freqMax, freqSteps)
+    
+    #dps = [-30e6]
+    
+    alfar = 90*(np.pi/180)
+    ex_r, ey_r, ez_r = np.sin(alfar)*np.cos(0), np.sin(alfar)*np.sin(0), np.cos(alfar)
+    
+    
+    alfap = 90*(np.pi/180)
+    ex_p, ey_p, ez_p = np.sin(alfap)*np.cos(0), np.sin(alfap)*np.sin(0), np.cos(alfap)
+    
+    """