no me deja pullear en mi compu si no hago esto

analisis/plots/20211115_CPT_DosLaseres_v04/CPT_plotter_20211115.py

@@ -531,6 +531,15 @@ plt.savefig('CPT_3lasers.svg')
 
 
 
+#%%
 
+i = 0
+for j in range(28):
+    plt.figure()
+    plt.plot([2*f*1e-6 for f in Freqs[j]], Counts[j], 'o-')
+    plt.xlabel('Frecuencia MHz)')
+    plt.ylabel('counts')
+    plt.grid()
+    plt.savefig(f'espectros_4_{j}.png',dpi = 200)
 
 

analisis/plots/20211201_CPT_DosLaseres_v05/CPT_plotter_20211201.py

@@ -358,4 +358,17 @@ plt.xlabel('Frecuencia (MHz)')
 plt.ylabel('counts')
 plt.grid()
 #plt.ylim(2500, 6100)
-plt.legend()
+plt.legend()
\ No newline at end of file
+
+
+#%%
+
+i = 0
+for j in range(28):
+    plt.figure()
+    plt.plot([2*f*1e-6 for f in Freqs[j]], Counts[j], 'o-')
+    plt.xlabel('Frecuencia MHz)')
+    plt.ylabel('counts')
+    plt.grid()
+    plt.savefig(f'espectros_5_{j}.png',dpi = 200)
+

analisis/plots/20220503_EspectrosUVnuevos/UV_spectrums2.py

@@ -9,17 +9,17 @@ import os
 from scipy import interpolate
 
 # Solo levanto algunos experimentos
-Calib_Files = """000007324-UV_Scan_withcalib_Haeffner
+Calib_Files = """000007209-UV_Scan_withcalib_Haeffner
 """ 
 
 
 
 #carpeta pc nico labo escritorio:
 #C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20220503_EspectrosUVnuevos\Data
-
+carpeta = "./Data/"
 def SeeKeys(files):
     for i, fname in enumerate(files.split()):
-        data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+        data = h5py.File(carpeta + fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
         print(fname)
         print(list(data['datasets'].keys()))
 
@@ -33,7 +33,7 @@ for i, fname in enumerate(Calib_Files.split()):
     print(SeeKeys(Calib_Files))
     print(i)
     print(fname)
-    data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+    data = h5py.File(carpeta + fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
     print(list(data['datasets'].keys()))
     Amps.append(np.array(data['datasets']['UV_Amplitudes']))
     Freqs.append(np.array(data['datasets']['UV_Frequencies']))

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/EITfit/MM_eightLevel_2repumps_AnalysisFunctions.py

@@ -12,24 +12,64 @@ import numpy as np
 import time
 import matplotlib.pyplot as plt
 from scipy.signal import argrelextrema
-import MM_eightLevel_2repumps_python_scripts
-from MM_eightLevel_2repumps_python_scripts import CPTspectrum8levels,CPTspectrum8levels_vel
+#import MM_eightLevel_2repumps_python_scripts
+#from MM_eightLevel_2repumps_python_scripts import CPTspectrum8levels,CPTspectrum8levels_vel
+from lolo_modelo_full_8niveles import CPTspectrum8levels_MM as CPTspectrum8levels
 from scipy.optimize import curve_fit
 import random
 from scipy.signal import savgol_filter as sf
 from scipy.stats import norm
-                             
 
+from numba import jit,njit
+
+
+@njit
+def trapz_numba(y, x):
+    result = 0.0
+    n = len(x)
+
+    for i in range(1, n):
+        result += (x[i] - x[i-1]) * (y[i] + y[i-1]) / 2.0
+
+    return result      
+
+@njit
+def trapz_numba_axis(y, x, axis):
+    """
+    Numerical integration using the trapezoidal rule along a specified axis.
+    
+    Parameters:
+    - x (array): Array of x values.
+    - y (array): Array of y values.
+    - axis (int): Axis along which to integrate.
+    
+    Returns:
+    - result (array): Result of the numerical integration along the specified axis.
+    """
+    result = np.zeros(np.shape(y)[1])
+    
+    # Iterate over the specified axis
+    for i in range(len(result)):
+        
+        result[i] = trapz_numba(y[:,i], x)
+    
+    return result
+
+                       
+@njit
 def prob_energia(E,T):
     kboltz = 1.380649e-23
     mcalcio = 6.655e-23*1e-3
     prob =  np.exp(-E/(kboltz*T)) #*E**2 
     return prob
 
+@njit
+def normal_pdf(v,sigma):
+    return np.exp(-v**2/(2*sigma**2)) /np.sqrt(2 * np.pi * sigma**2)
 
 
-
-def PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = False):
+@njit
+def PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = False):
     """
     solvemode=1: resuelve con np.linalg.solve
     solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
@@ -40,36 +80,39 @@ def PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinew
     kp = np.pi * 2/866 *1e9
     if detpvec is None:
         detpvec = np.arange(freqMin,freqMax,freqStep)    
-
+    
     if boltzmann == True:
-        FluorescencesVel = []
+        
         sigmaT = np.sqrt(kboltz*T/m)
         
         velvec = np.linspace(-4*sigmaT,4*sigmaT,100)
-        MBprobVec = norm.pdf(velvec,loc = 0, scale = sigmaT)
+        MBprobVec = normal_pdf(velvec,sigmaT)
         MBprobVec = MBprobVec/np.trapz(MBprobVec,velvec)
+        FluorescencesVel = np.zeros((len(velvec),len(detpvec)))
         for i in range(len(velvec)):
             detVel = detpvec - kp*velvec[i]/(2*np.pi*1e6)
-            _, Fluorescence = CPTspectrum8levels( sg, sp, gPS, gPD, DetDoppler-kg*velvec[i]/(2*np.pi*1e6),u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0,alpha, phidoppler, titadoppler, phiprobe, titaprobe,circularityprobe, betag, betap, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False,detpvec=detVel)
+            _, Fluorescence = CPTspectrum8levels( sg, sp, gPS, gPD, DetDoppler-kg*velvec[i]/(2*np.pi*1e6),u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0,alpha, phidoppler, titadoppler, phiprobe, titaprobe,circularityprobe, betag, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False,detpvec=detVel)
             
-            FluorescencesVel.append(Fluorescence)    
-        FluorescencesVel = np.array(FluorescencesVel)
-        MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
-        Fluorescence = np.trapz(FluorescencesVel*MBprobMat,velvec,axis = 0)
+            FluorescencesVel[i,:] = Fluorescence    
+        
+        #MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
+        MBprobMat = MBprobVec.repeat(FluorescencesVel.shape[1]).reshape((-1, FluorescencesVel.shape[1]))
+        Fluorescence = trapz_numba_axis(FluorescencesVel*MBprobMat,velvec,0)
         ProbeDetuningVectorL = np.arange(freqMin,freqMax,freqStep)
 
     elif dephasing == True:  
-        tinicial = time.time()
-        ProbeDetuningVectorL, Fluorescence = CPTspectrum8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe, betag, betap, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, detpvec = detpvec, plot=False, solvemode=1)
-        tfinal = time.time()
+        print('holaa')
+        #tinicial = time.time()
+        ProbeDetuningVectorL, Fluorescence = CPTspectrum8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, circularityprobe, betag, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, detpvec = detpvec, plot=False, solvemode=1)
+        #tfinal = time.time()
         
 
     else:
-        drivefreq =2*np.pi* 0.6e6
-        FluorescencesVel = []
-        sigmaT = np.sqrt(T*kboltz/m)
-        velvec = np.linspace(-4*sigmaT,4*sigmaT,100)
-        MBprobVec = norm.pdf(velvec,loc = 0, scale = sigmaT)
+        
+        
+        # sigmaT = np.sqrt(T*kboltz/m)
+        # velvec = np.linspace(-4*sigmaT,4*sigmaT,100)
+        # MBprobVec = norm.pdf(velvec,loc = 0, scale = sigmaT)
         
         Evec = np.linspace(0,8*kboltz*T,100)
         velvec = np.zeros(2*len(Evec)) 
@@ -78,22 +121,461 @@ def PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinew
         MBprobVec = prob_energia(Evec, T)
         MBprobVec = MBprobVec/np.trapz(MBprobVec,Evec)
         
-        betagVec = velvec*kg/drivefreq
-        betapVec = velvec*kp/drivefreq
-        for i in range(len(betagVec)): 
-            betag,betap = betagVec[i],betapVec[i]
-                                                         
-            Frequencies, Fluorescence = CPTspectrum8levels( sg, sp, gPS, gPD, DetDoppler,u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0,alpha, phidoppler, titadoppler, phiprobe, titaprobe,circularityprobe, betag, betap, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False,detpvec=detpvec)
-            FluorescencesVel.append(Fluorescence)
-        FluorescencesVel = np.array(FluorescencesVel)
-        MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
-        Fluorescence = np.trapz(FluorescencesVel*MBprobMat,Evec,axis = 0)
+        kgvelVec = velvec*kg
+        
+       
+        FluorescencesVel = np.zeros((len(velvec),len(detpvec)))
+       
+        for i in range(len(kgvelVec)): 
+            #tinicial = time.time()
+            kgvel = kgvelVec[i]
+            #print(i)                                            
+            Frequencies, Fluorescence = CPTspectrum8levels( sg, sp, gPS, gPD, DetDoppler,u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0,alpha, phidoppler, titadoppler, phiprobe, titaprobe,circularityprobe, kgvel, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False,detpvec=detpvec)
+            FluorescencesVel[i,:] = Fluorescence
+            #tfinal = time.time()
+            #print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
+        #MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
+        #print('aca')
+        MBprobMat = MBprobVec.repeat(FluorescencesVel.shape[1]).reshape((-1, FluorescencesVel.shape[1]))
+
+        #print(np.shape(MBprobMat))
+        #print(np.shape(FluorescencesVel))
+        #print('aca 2')
+        Fluorescence = trapz_numba_axis(FluorescencesVel*MBprobMat,Evec,0)
+        #print('aca 3')
         ProbeDetuningVectorL = Frequencies
-    
-    #print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
+        
+        #print('aca 4')
           
     return ProbeDetuningVectorL, Fluorescence
 
+
+@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/866*397
+    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 = int(beta)
+
+    #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(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,detpvec = None):
+
+    """
+    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)
+    DetProbeVector = 2*np.pi*detpvec*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, np.array(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')
+
+'''
 def PerformExperiment_8levels_vel(velvect,titavect,phivect,velprob,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
@@ -145,5 +627,5 @@ def try_fitCPT_8levels(A,bgnd,f0,Freq,Fluo,sg,sp,gPS,gPD,DetDoppler,u,DopplerLas
         return spectra*A + bgnd
     return SpectrumForFit(Freq, sg, sp, T, DetDoppler, A, bgnd, f0)
     
-    
+'''  
     
\ No newline at end of file

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/EITfit/MM_eightLevel_2repumps_CPTPlotter.py

@@ -10,13 +10,16 @@ Created on Tue Sep  1 17:58:39 2020
 
 import os
 import numpy as np
-from MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels, GenerateNoisyCPT_vel, SmoothNoisyCPT, PerformExperiment_8levels_vel,fitCPT_8levels,try_fitCPT_8levels
+from lolo_modelo_full_8niveles import Fluorescencia_RK, Fluorescencia_prueba
+from MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels#, GenerateNoisyCPT_vel, SmoothNoisyCPT, PerformExperiment_8levels_vel,fitCPT_8levels,try_fitCPT_8levels
 import matplotlib.pyplot as plt
 import time
 import h5py
-
+from numba import njit,jit
 os.chdir('/home/muri/nubeDF/Documents/codigos/artiq_experiments/analisis/plots/20231218_CPT_muri/Data')
 
+os.chdir('/home/muri/nubeDF/Documents/codigos/artiq_experiments/analisis/plots/20230920_CPT_TemperatureSens_v2/Data')
+
 CPT_FILES = """000016262-IR_Scan_withcal_optimized
 000016239-IR_Scan_withcal_optimized
 000016240-IR_Scan_withcal_optimized
@@ -25,9 +28,37 @@ CPT_FILES = """000016262-IR_Scan_withcal_optimized
 000016255-IR_Scan_withcal_optimized
 000016256-IR_Scan_withcal_optimized
 000016257-IR_Scan_withcal_optimized
-""" 
+"""
 
 
+CPT_FILES = """000015300-IR_Scan_withcal_optimized
+000015301-IR_Scan_withcal_optimized
+000015302-IR_Scan_withcal_optimized
+"""
+'''
+CPT_FILES = """000016645-IR_Scan_withcal_optimized
+000016646-IR_Scan_withcal_optimized
+000016647-IR_Scan_withcal_optimized
+000016648-IR_Scan_withcal_optimized
+000016649-IR_Scan_withcal_optimized
+000016650-IR_Scan_withcal_optimized
+000016655-IR_Scan_withcal_optimized
+000016656-IR_Scan_withcal_optimized
+000016657-IR_Scan_withcal_optimized
+000016658-IR_Scan_withcal_optimized
+000016659-IR_Scan_withcal_optimized
+000016660-IR_Scan_withcal_optimized
+000016661-IR_Scan_withcal_optimized
+000016662-IR_Scan_withcal_optimized
+000016664-IR_Scan_withcal_optimized
+000016666-IR_Scan_withcal_optimized
+000016669-IR_Scan_withcal_optimized
+000016670-IR_Scan_withcal_optimized
+000016672-IR_Scan_withcal_optimized
+000016673-IR_Scan_withcal_optimized
+"""
+'''
+
 def SeeKeys(files):
     for i, fname in enumerate(files.split()):
         data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
@@ -46,6 +77,7 @@ AmpTisa = []
 UVCPTAmp = []
 No_measures = []
 Voltages = []
+Heating = []
 
 for i, fname in enumerate(CPT_FILES.split()):
     print(str(i) + ' - ' + fname)
@@ -56,11 +88,13 @@ for i, fname in enumerate(CPT_FILES.split()):
     # que además tenian un error de tipeo al final. Esto no deberá ser necesario
     # cuando se solucione el error este del guardado.
     Freqs.append(np.array(data['datasets']['IR1_Frequencies']))
-    Counts.append(np.array(data['datasets']['data_array']))
+    #Counts.append(np.array(data['datasets']['data_array']))
+    Counts.append(np.array(data['datasets']['counts_spectrum']))
     #AmpTisa.append(np.array(data['datasets']['TISA_CPT_amp']))
     UVCPTAmp.append(np.array(data['datasets']['UV_CPT_amp']))
     No_measures.append(np.array(data['datasets']['no_measures']))
-    Voltages.append(np.array(data['datasets']['scanning_voltages']))
+    Heating.append(np.array(data['datasets']['heating']))
+    #Voltages.append(np.array(data['datasets']['scanning_voltages']))
 
 def Split(array,n):
     length=len(array)/n
@@ -80,9 +114,8 @@ def Split(array,n):
 CountsSplit = []
 CountsSplit.append(Split(Counts[0],len(Freqs[0])))
 
-
-CountsSplit_2ions = []
-CountsSplit_2ions.append(Split(Counts[4],len(Freqs[4])))
+#CountsSplit_2ions = []
+#CountsSplit_2ions.append(Split(Counts[4],len(Freqs[4])))
 
 #%%
 
@@ -92,7 +125,8 @@ Las que valen la pena son de la 1 a la 9.
 En particular, la 4 tiene poca micromoción: 
 """
 
-jvec = [4] # de la 1 a la 9 vale la pena, despues no
+CountsSplit = Counts
+jvec = [0,1,2] # de la 1 a la 9 vale la pena, despues no
 
 drive=22.1
 
@@ -101,7 +135,7 @@ Frequencies = Freqs[0]
 plt.figure()
 i = 0
 for j in jvec:
-    #plt.errorbar([2*f*1e-6-419-13 for f in Frequencies], CountsSplit[0][j], yerr=np.sqrt(CountsSplit[0][j]), fmt='o', capsize=2, markersize=2)
+    plt.errorbar([2*f*1e-6-440 for f in Frequencies], CountsSplit[j], yerr=np.sqrt(CountsSplit[j]), fmt='o', capsize=2, markersize=2,label = '{:}'.format(j))
     i = i + 1
 plt.xlabel('Frecuencia (MHz)')
 plt.ylabel('counts')
@@ -111,7 +145,7 @@ plt.grid()
     #plt.axvline(dr+drive)
 plt.legend()
 
-
+#%%
 plt.rcParams["axes.prop_cycle"] = plt.cycler('color', ['#DBAD1F', '#C213DB', '#DB4F2A', '#0500DB', '#09DB9B', '#B4001B', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf'])
 
 
@@ -128,12 +162,12 @@ 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.5 #ancho de linea de los laseres en MHz
+lw = 0.1 #ancho de linea de los laseres en MHz
 DopplerLaserLinewidth, ProbeLaserLinewidth = lw, lw #ancho de linea de los laseres
 
-DetDoppler = -24.5 #detuning doppler en MHz
+DetDoppler = -10 #detuning doppler en MHz
 
-T = 0.1e-3 #temperatura en K
+T = 0 #temperatura en K
 alpha = 0 #angulo entre los láseres
 
 #estos son los angulos de la polarizacion de los laseres respecto al campo magnetico
@@ -156,11 +190,14 @@ print(freqMin,freqMax,freqStep)
 noiseamplitude = 0
 
 #parametros de saturacion de los laseres. g: doppler. p: probe (un rebombeo que scanea), r: repump (otro rebombeo fijo)
-sg = 0.25
-sp = 5.0
-
+sg = 0.4
+sp = 8
 drivefreq=2*np.pi*22.135*1e6
 
+rabG = gPS*sg
+rabP =  gPD*sp
+
+
 kg = np.pi * 2/397 *1e9
 kp = np.pi * 2/866 *1e9
 
@@ -174,24 +211,8 @@ betavec=[0]
 FrequenciesVec = []
 FluorescencesVec = []
 
-alphavec = [0]
-Tvec = np.linspace(0.5e-3,100e-3,15)
-#Tvec = np.linspace(0.001,0.1,20)
-s12vec = np.linspace(0.2,0.6,20)
-#s12vec = [s12vec[2],s12vec[4],s12vec[10],s12vec[18]]
-s12vec = [0.51]
-s23vec = np.linspace(0.5,15,20)
-#s23vec = [s23vec[2],s23vec[6],s23vec[10],s23vec[18]]
-s23vec = [3.5]
-
-Tmat = np.zeros((len(s12vec),len(s23vec)))
+T = 20e-3
 
-Tvec_ajuste = np.zeros(len(Tvec))
-phivect = 0
-
-T = 8e-3
-
-titavect = np.linspace(0,2*np.pi,100)    
   
 betag = 0
 betap = 0 
@@ -202,21 +223,42 @@ SCALE = 1
 
 
 OFFSET = 0
-OFFSET = 1.55401398e+02
+#OFFSET = 1.55401398e+02
+
+drivefreq =2*np.pi*0.6*1e6
+betag = 0
 
+dt = 1e-9
+tfinal = 2e-5
+N = int(tfinal/dt)
 
 
-Freq,Fluo_sb = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = False)
-Freq,Fluo_boltz = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = True)
-T = 0.8e-3
-Freq,Fluo_deph = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, betap, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = True, boltzmann = False)
+Freq,Fluo_sb = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = False)
+#Freq,Fluo_boltz = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = False, boltzmann = True)
 
+#Freq,Fluo_deph = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None,dephasing = True)
+Freq,Fluo0 = PerformExperiment_8levels(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, 0, alpha, phidoppler, titadoppler, phiprobe, titaprobe, betag, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, dephasing  = True, boltzmann = False)
+Freq, Fluo_full = Fluorescencia_RK(rabG, rabP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, phidoppler, titadoppler, phiprobe, titaprobe, betag,drivefreq ,freqMin=freqMin, freqMax=freqMax, freqStep=freqStep )
+#Freq, Fluo_full = Fluorescencia_prueba(dt,N,rabG, rabP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, phidoppler, titadoppler, phiprobe, titaprobe, betag,drivefreq ,freqMin=freqMin, freqMax=freqMax, freqStep=freqStep )
 
-plt.plot(Freq,SCALE*Fluo_boltz+OFFSET,label = 'Semiclassic')
-plt.plot(Freq,SCALE*Fluo_sb+OFFSET,label = 'Sideband')
-plt.plot(Freq,SCALE*Fluo_deph+OFFSET, label = 'Dephasing')
 
+plt.figure()
+plt.plot(Freq,SCALE*Fluo0+OFFSET, label = 'T = 0')
+#plt.plot(Freq,SCALE*Fluo_deph+OFFSET,':', label = 'Dephasing')
+#plt.plot(Freq,SCALE*Fluo_boltz+OFFSET ,label = 'Instantaneous')
+plt.plot(Freq,SCALE*Fluo_sb+OFFSET,'--',label = 'Sideband')
+
+
+
+plt.xlabel('Frecuencia (MHz)')
+plt.ylabel('Fluorescencia (U.A.)')
+plt.grid()
+#for dr in drs:
+#    plt.axvline(dr)
+    #plt.axvline(dr+drive)
 plt.legend()
+plt.savefig('20mK.png',dpi = 300)
+
 
 
 
@@ -246,11 +288,12 @@ u = 32.5e6
 
 #B = (u/(2*np.pi))/c
 
-correccion = 13
+correccion = 30
 
 offsetxpi = 419+correccion
+offsetxpi = 443
 DetDoppler = -11.5-correccion
-
+DetDoppler = -25
 gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 
 alpha = 0
 
@@ -258,44 +301,56 @@ alpha = 0
 drivefreq = 2*np.pi*22.135*1e6
 
 
-selectedcurve = 4
+selectedcurve = 0
 
-FreqsDR = [2*f*1e-6-offsetxpi for f in Freqs[0]]
-CountsDR = CountsSplit[0][selectedcurve]
-CountsDR[100]=0.5*(CountsDR[99]+CountsDR[101])
-CountsDR[105]=0.5*(CountsDR[104]+CountsDR[106])
+FreqsDR = np.array([2*f*1e-6-offsetxpi for f in Freqs[0]])
+CountsDR = CountsSplit[selectedcurve]
+#CountsDR[100]=0.5*(CountsDR[99]+CountsDR[101])
+#CountsDR[105]=0.5*(CountsDR[104]+CountsDR[106])
 
 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, SG, SP, SCALE1, OFFSET, TEMP,f0):
+@njit
+def FitEIT_MM_single(freqs, SG, SP, SCALE1, OFFSET, TEMP,LW,f0,DetDoppler):
 #def FitEIT_MM(freqs, SG, SP, SCALE1, OFFSET, BETA1):
     #BETA = 1.8
     # SG = 0.6
     # SP = 8.1
     # TEMP = 0.2e-3
     BETA1 = 0
-    BETA2 = 0
-    Detunings, Fluorescence1 = PerformExperiment_8levels(SG, SP, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, TEMP, alpha, phidoppler, titadoppler, phiprobe, titaprobe,  BETA1, BETA2, drivefreq, min(freqs), max(freqs)+(freqs[1]-freqs[0]), freqs[1]-freqs[0], circularityprobe=CircPr, plot=False, solvemode=1, detpvec=np.array(freqs)-f0)
+    freqs = freqs - f0
+    ProbeLaserLinewidth = LW
+    DopplerLaserLinewidth = LW
+    Detunings, Fluorescence1 = PerformExperiment_8levels(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,detpvec= freqs, solvemode=1,dephasing = True)
 
-    ScaledFluo1 = np.array([f*SCALE1 + OFFSET for f in Fluorescence1])
-    return ScaledFluo1
+    ScaledFluo1 = Fluorescence1*SCALE1 + OFFSET
     #return ScaledFluo1
+    return ScaledFluo1
+    
 
 popt_1 = [4.82483692e-01, 8.18685518e+00, 6.48238914e+04, 1.54835760e+02, 5.56566519e-03, 5.24515673e-18]
+popt_1 = [4.82479600e-01, 8.18689939e+00, 6.48254299e+04, 1.54802000e+02,5.56602117e-03, 2.60257233e-18]
+popt_1 =[5.28722824e-01, 8.45443080e+00, 5.45497905e+04, 2.17427415e+02,5.81901418e-14, 1.73736544e+00, 1.19825202e-17]
+popt_1 =[5.28722824e-01, 8.45443080e+00, 5.45497905e+04, 2.17427415e+02,0.5e-3, 0.5, 1.19825202e-17]
+popt_1 = [4.46477838e-01, 8.65075889e+00, 5.09989857e+04, 2.72249632e+02,4.97930986e-04, 5.00010000e-01, 6.18443753e-01]
+popt_1 = [ 0.2,  5,  50000,  150,1e-3,  0.1,  -15.68376255e+00, -1.40024175e+01]
 
+#medicion 7
+popt_1 = [ 5.05319738e-01,  8.25036863e+00,  1.15504507e+04,  8.78454857e+01,6.66853533e-04,  9.99000000e-02, -1.73285475e+01, -1.09633879e+01]
+
+popt_1 = [0.5, 15, 11000, 80,1e-3, 0.1, -18,-10]
 do_fit = True
 if do_fit:
-    #popt_1, pcov_1 = curve_fit(FitEIT_MM_single, FreqsDR, CountsDR, p0=[0.2, 5, 6.89e4, 1.2723, 8e-3,1], bounds=((0, 0, 0, 0, 0, 0), (2, 20, 7e4, 5e4, 15e-3,5)))
+    #popt_1, pcov_1 = curve_fit(FitEIT_MM_single, FreqsDR, CountsDR, p0=[0.3, 7, 11000, 150,5e-3, 0.1, -18,-7], bounds=((0, 0, 0, 0, 0.1e-3,0.0999, -100,-15), (2, 20, 1e5, 5e4, 5e-3,0.10001,-15,-5)))
     FittedEITpi_1 = FitEIT_MM_single(freqslong, *popt_1)
 
 beta1 = popt_1[4]
-errorbeta1 = np.sqrt(pcov_1[4,4])
+#errorbeta1 = np.sqrt(pcov_1[4,4])
 temp1 = popt_1[5]
-errortemp1 = np.sqrt(pcov_1[5,5])
+#errortemp1 = np.sqrt(pcov_1[5,5])
 
 plt.figure()
 plt.errorbar(FreqsDR, CountsDR, yerr=2*np.sqrt(CountsDR), fmt='o', color='darkgreen', alpha=0.5, capsize=2, markersize=2)
@@ -324,7 +379,6 @@ plt.grid()
 
 
 
-
 
 
 

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/EITfit/MM_eightLevel_2repumps_python_scripts.py

@@ -315,7 +315,7 @@ def FullL(rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp
     
     if betag !=0 and betap !=0:
         #ESTA PARTE ES CUANDO AGREGAS MICROMOCION
-        nmax = 2*int(betag)
+        nmax = int(betag)
         #print(nmax)
         Ltemp, Omega = LtempCalculus(betag,betap, drivefreq)
         #print(factor)

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/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
+
+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,complex128,float64
+
+@njit
+def PerformExperiment_8levels_MM(sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, kgvel, 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, kgvel, drivefreq,
+                                freqMin=freqMin, freqMax=freqMax, freqStep=freqStep,
+                                plot=False, solvemode=1,detpvec = detpvec)
+    ProbeDetuningVectorL, Fluovector = rta
+    return ProbeDetuningVectorL, Fluovector
+
+
+
+#%% Estos son los auxiliares ###################################################
+
+"""
+Esta parte es la del modelo
+"""
+
+@njit
+def trapz_numba(y, x):
+    result = 0.0
+    n = len(x)
+
+    for i in range(1, n):
+        result += (x[i] - x[i-1]) * (y[i] + y[i-1]) / 2.0
+
+    return result      
+
+
+@njit
+def dy_dx(t,rho,rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp = 0,
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, kgvel=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.
+    """
+    
+    w1 = 0.6e6*2*np.pi
+    w2 = 0.8e6*2*np.pi
+    w3 = 1.2e6*2*np.pi
+    
+    norm = np.array([1,1,1],dtype = float64)/np.sqrt(3)
+    temp_part  =np.array([np.cos(w1*t), np.cos(w2*t), np.cos(w3*t)])
+    
+    Detg = Detg - kgvel*np.dot(norm,temp_part)    
+    Detp = Detp - kgvel*397/866*np.dot(norm,temp_part)
+    
+    
+    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=complex128)
+
+    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
+    der = np.dot(L0,rho)
+    return der
+
+
+
+@njit
+def runge_kutta_order_4(dy_dx, N, h,rabG, rabP, gPS = 0, gPD = 0, Detg = 0, Detp = 0, u = 0, lwg = 0, lwp = 0,
+          phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, kgvel=0,
+          drivefreq=2*np.pi*22.135*1e6, T = 0, alpha = 0, circularityprobe=1):
+    """
+    Solve the initial value problem using the fourth-order Runge-Kutta method.
+
+    Parameters:
+    - dy_dx: A function representing the derivative of y with respect to x.
+    - y0: Initial value of y at x0.
+    - x0: Initial value of x.
+    - xn: Final value of x.
+    - h: Step size.
+
+    Returns:
+    - x_values: List of x values.
+    - y_values: List of corresponding y values.
+    """
+    y0 = np.zeros(64,dtype = complex128).T
+    y0[0] = 1
+    x_values = np.zeros(N+1)
+    #y_values = []
+
+    y_values = np.zeros((N + 1, y0.shape[0]), dtype=complex128)
+    x =0
+    y = y0
+    x_values[0] = x  
+    y_values[0] = y
+    
+    
+    
+    for i in range(N):
+        
+        k1 =  dy_dx(x, y,rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq, T, alpha, circularityprobe)
+        k2 =  dy_dx(x + 0.5 * h, y + h * 0.5 * k1,rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq, T, alpha, circularityprobe)
+        k3 =  dy_dx(x + 0.5 * h, y + h * 0.5 * k2,rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq, T, alpha, circularityprobe)
+        k4 =  dy_dx(x + h, y + h * k3,rabG, rabP, gPS, gPD, Detg, Detp, u, lwg, lwp,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq, T, alpha, circularityprobe)
+
+        y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
+        x += h
+
+        x_values[i+1] = x  
+        y_values[i+1] = y
+
+    return x_values, y_values
+
+@njit
+def Fluorescencia_RK(rabG, rabP, gPS, gPD , Detg, u , lwg, lwp ,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq=2*np.pi*0.6*1e6,freqMin=-50, freqMax=50, freqStep=1):
+         
+    
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+
+    
+    Detg = 2*np.pi*1e6*Detg
+    lwg, lwp = 2*np.pi*1e6*lwg, 2*np.pi*1e6*lwp
+    
+    tf = 2e-5
+    dt = 1e-9
+    N = int(tf/dt)
+ 
+    detuningvector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
+    Fluorescence = np.zeros(len(detuningvector),dtype = complex128)
+    i = 0
+    for detr in detuningvector:
+        tiempos,rho = runge_kutta_order_4(dy_dx,N,dt,rabG, rabP, gPS, gPD , Detg , detr , u , lwg, lwp ,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq)
+        rho_est = np.mean(rho[-int(2*2*np.pi/(drivefreq*dt)):,18]+rho[-int(2*2*np.pi/(drivefreq*dt)):,27])        
+        Fluorescence[i] =rho_est
+        i = i+1
+    Frequencies = detuningvector/1e6/2/np.pi
+    
+    return Frequencies,Fluorescence
+
+def Fluorescencia_prueba(dt,N,rabG, rabP, gPS, gPD , Detg, u , lwg, lwp ,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq=2*np.pi*0.6*1e6,detuningvector = [0]):
+         
+    phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
+    phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
+    
+    Detg = 2*np.pi*1e6*Detg
+    lwg, lwp = 2*np.pi*1e6*lwg, 2*np.pi*1e6*lwp
+    
+    
+
+    Fluorescence = np.zeros(len(detuningvector))
+    i = 0
+    for detr in detuningvector:
+        tiempos,rho = runge_kutta_order_4(dy_dx,N,dt,rabG, rabP, gPS, gPD , Detg , detr , u , lwg, lwp ,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq)
+        rho_est = np.mean(rho[-int(2*2*np.pi/(drivefreq*dt)):,18]+rho[-int(2*2*np.pi/(drivefreq*dt)):,27])        
+        plt.plot(tiempos[:],rho[:,18]+rho[:,27])
+        Fluorescence[i] =rho_est
+        i = i+1
+    Frequencies = detuningvector/1e6/2/np.pi
+    
+    return Frequencies,Fluorescence,tiempos,rho
+
+
+
+@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=complex128)
+
+
+    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(kgvel:float, drivefreq:float, forma=1):
+    Hint = np.zeros((8, 8), dtype=np.complex_)
+
+    ampg=kgvel
+    ampr=kgvel/866*397
+    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=complex128)
+
+    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=complex128)
+
+    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, kgvel=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=complex128)
+
+    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
+    if kgvel != 0:
+        #ESTA PARTE ES CUANDO AGREGAS MICROMOCION
+        nmax = int(kgvel/drivefreq)
+    
+        #print(nmax)
+        Ltemp, Omega = LtempCalculus(kgvel, drivefreq)
+        #print(factor)
+        L1 = GetL1(Ltemp, L0, Omega, nmax)
+        Lfull = L0 + L1 #ESA CORRECCION ESTA EN L1
+        #HASTA ACA
+    else:
+        Lfull = L0
+    #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, kgvel, drivefreq, freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1,detpvec = None):
+
+    """
+    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)
+    #DetProbeVector = 2*np.pi*detpvec*1e6
+    
+    Detg = 2*np.pi*Detg*1e6
+    
+   
+    lwg, lwp = 2*np.pi*lwg*1e6, 2*np.pi*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, kgvel, drivefreq, Temp, alpha, Circularityprobe)
+        # if solvemode == 1:
+        rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)],dtype=complex128))
+        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, np.array(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
+    kboltz = 1.38e-23
+    m = 6.65e-26
+    u = 32e6 #esto equivale aprox al campo que tenemos
+
+    B = (u/(2*np.pi))/c
+    
+
+    
+    sg, sr, sp = 0.5, 8, 8 #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.1, 0.1, 0.1 #ancho de linea de los laseres
+    lwg,lwr,lwp = 0.1,0.1,0.1
+    kg = 2*np.pi/397e-9
+    Detg = -10
+    Detr = 0 #detuning del doppler y repump
+    Detp = -18
+    
+    Temp = 20e-3 #temperatura en K
+    vel = np.sqrt(kboltz*Temp/m)
+    
+    kgvel = kg*vel
+    
+    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 = -30
+    freqMax = 10
+    freqStep = 0.5
+    
+    
+    
+    dt = 1e-9
+    tfinal = 5e-5
+    N = int(tfinal/dt)
+    
+    drivefreq = 0.6e6*2*np.pi
+    Circularityprobe = 1
+    
+    detvec = np.arange(freqMin,freqMax,freqStep)*2*np.pi*1e6
+    #detvec = np.array([-21.25])*1e6*2*np.pi
+    
+    
+    
+   
+    
+    Frequencyvector_MM, Fluovector_MM = CPTspectrum8levels_MM(sg, sp, gPS, gPD, Detg, u, lwg, lwp, 0, alpha, phidoppler, titadoppler, phiprobe, titaprobe, Circularityprobe, kgvel, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep)
+    Frequencyvector_MM, Fluovector_0 = CPTspectrum8levels_MM(sg, sp, gPS, gPD, Detg, u, lwg, lwp, 0, alpha, phidoppler, titadoppler, phiprobe, titaprobe, Circularityprobe, 0, drivefreq, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep)
+    plt.figure()
+    tic = time.time()
+    #Frequencyvector, Fluovector,tiempo,rho = Fluorescencia_prueba(dt,N,rabG, rabP, gPS, gPD, Detg, u, lwg, lwp, phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq =2*np.pi*0.6*1e6,detuningvector=detvec)
+    Frequencyvector, Fluovector= Fluorescencia_RK(rabG, rabP, gPS, gPD, Detg, u, lwg, lwp, phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq =drivefreq,freqMin=freqMin, freqMax=freqMax, freqStep=freqStep )
+    #tiempo,rho = runge_kutta_order_4(dy_dx, N, dt,rabG, rabP, gPS, gPD , Detg , Detp , u , lwg, lwp ,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq=2*np.pi*0.6*1e6)
+    tac = time.time()
+    
+    #plt.plot(tiempo,rho[:,18]+rho[:,27])
+    #plt.plot(tiempo,rho[:,18])
+    #plt.plot(tiempo,rho[:,27])
+    plt.figure()
+    plt.plot(Frequencyvector_MM, [100*f for f in Fluovector_MM],label = 'Sidebands')
+    plt.plot(Frequencyvector_MM, [100*f for f in Fluovector_0],label = 'T = 0 mK')
+    plt.plot(Frequencyvector, [100*f for f in Fluovector], label='OBS')
+    plt.grid()
+    plt.legend()    
+    plt.xlabel('Probe detuning (MHz)')
+    plt.ylabel('Fluorescence (A.U.)')
+
+    #Lfull = FullL_MM(rabG, rabP, gPS, gPD, Detg, detvec[0], u, lwg , lwp ,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq, Temp , alpha, Circularityprobe)
+    #Mfull = dy_dx(0,0,rabG, rabP, gPS, gPD, Detg*1e6*2*np.pi, detvec[-1], u, lwg , lwp ,phidoppler, titadoppler, phiprobe, titaprobe, kgvel,drivefreq, Temp , alpha, Circularityprobe)