ion statistics

analisis/plots/20221017_IonStatistics/IonStatistics.py

+import h5py
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+import re
+import ast
+from scipy.optimize import curve_fit
+import os
+
+# Solo levanto algunos experimentos
+Stat_files = [8731]
+
+def expo(T, tau, N0, C):
+    global T0
+    return N0*np.exp(-(T-T0)/tau) + C
+
+def pow_from_amp(amp):
+    """Paso de amplitud urukul a potencia medida por Nico"""
+    # Forma altamente ineficiente de hacer esto, pero me salio asi
+    amplitudes_UV = np.flip(np.array([0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.22, 0.24, 0.26, 0.28, 0.30]))
+    assert amp in amplitudes_UV
+    potencias_UV = np.flip(np.array([4, 10, 19, 32, 49, 71, 96, 125, 155, 183, 208, 229]))
+    return potencias_UV[np.where(amplitudes_UV == amp)][0]
+
+"""
+plt.plot(amplitudes_UV, potencias_UV, 'ko-', lw=0.2)
+plt.xlabel("Amplitud Urukul")
+plt.ylabel("Potencia /uW")
+plt.grid()
+"""
+#%%
+
+BINW = 10e-9
+T0 = 0.0e-6
+
+Stat_Heigths = []
+Stat_Bins = []
+
+for i, fname in enumerate(Stat_files):
+    #print(i)
+    #print(fname)
+    data = h5py.File('Data/00000'+str(fname)+'-IonStatistics.h5', 'r')
+    counts = np.array(data['datasets']['counts'])
+    bines = np.arange(counts.min(), counts.max()+BINW, BINW)
+
+    heigs, binsf  = np.histogram(counts, bines[bines>T0])
+    
+    Stat_Heigths.append(heigs)
+    Stat_Bins.append(binsf)
+    
+#%%
+
+#plt.figure()
+#plt.plot(Stat_Bins[0][:-1], Stat_Heigths[0])
+
+plt.figure()
+plt.hist(Stat_Heigths[0], bins=np.arange(100,350, 1), histtype='step')
+
+
+
+

analisis/plots/20221017_IonStatistics/Simulaciones/OBE3levels_chtime_fitting.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import numpy as np
+import matplotlib.pyplot as plt
+from OBE3levels_functions import *
+from scipy.optimize import curve_fit
+from scipy.stats import poisson, norm
+import os
+import sys
+plt.ion()
+"""
+Este código simula la dinámica temporal de un ion de calcio sometido a un láser UV y dos láseres IR.
+
+"""
+def expo(T, tau, N0, C):
+    #global T0
+    return N0*np.exp(-(T-0e-6)/tau) + C
+
+ADD_NOISE = True
+NOISE_VAR = 30 # Varianza del ruido
+NOISE_ADJ = 78 # Escaleo del ruido versus el maximo
+NOISE_MEAN = 0 # Media del ruido
+NRO_REPS = 100 # Veces a promediar
+
+species = 'Calcium'
+#Parámetros de saturacion de los láseres rebombeo (IR), que son proporcionales a la intensidad de los mismos
+#Para reproducir SP, hay que poner 0 a los dos
+SatParRep1, SatParRep2 = 0, 0
+
+sat_intensity = 43.3e-5 # uW um-2
+g21 = 135591138.92893547
+
+#El número de los láseres IR aplicados en el experimento. Para el experimento de la transicion SP, es 0.
+Number_IR_lasers = 0
+
+#Detunings de los láseres en MHz. Si es 0, están en resonancia. En principio
+#podría ser un parámetro libre a ajustar porque influye mucho
+DetDoppler = -int(sys.argv[1])
+DetRep2 = 0
+DetRep1 = 0
+
+#Anchos de línea de los láseres, en MHz. No son tan cruciales, están entre 0.1 y 0.5 MHz
+DopplerLaserLinewidth, RepumpLaserLinewidth = 1e3, 1
+
+# Para trabajar con las curvas teóricas:
+# range_inicio = np.arange(sat_intensity/10, sat_intensity, sat_intensity/50)
+# range_fin = np.arange(sat_intensity, 30*sat_intensity, sat_intensity/5)
+# laser_intensities = np.append(range_inicio, range_fin)
+
+# Para trabajar con nuestros datos, saco las intensidades, a un radio de 100um:
+laser_intensities = np.array([0.00662085, 0.00582507, 0.0049338, 0.00397887, 0.00305577, 0.00226, 0.00155972, 0.00101859, 0.00060479, 0.00031831])
+
+
+OmegaRabi2 = (g21**2) * (laser_intensities / sat_intensity)
+SatParDopplerVec = OmegaRabi2 / ((DetDoppler*2*np.pi*1e6)**2 + g21**2)
+
+#Parámetros de la simulación, en us
+tfinOBE = 1
+pasoOBE = 1e-5
+
+#Estado electrónico inicial de la simulacion. Posibilidades: |S> = 1, |P> = 2, |D> = 3
+initial_state = 1
+
+#Corre y simula la dinámica. Devuelve elementos de matriz densidad
+# fig1, ax1 = plt.subplots()
+
+alltaus = list()
+allN0 = list()
+allmeans = list()
+allstds = list()
+
+temporal_taus = list()
+temporal_N0 = list()
+
+for i, SatParDoppler in enumerate(SatParDopplerVec):
+    t, RealRho11, RealRho22, RealRho33, AbsRho12, AbsRho13, AbsRho23 = Solve3LevelBloch(EquationSystemModel, species, SatParDoppler, SatParRep1, SatParRep2, DetDoppler, DetRep1, DetRep2, DopplerLaserLinewidth, RepumpLaserLinewidth, no_IR_lasers=Number_IR_lasers, initcond=initial_state, tfin=tfinOBE, paso=pasoOBE, printprogress=False)
+
+    k = int(0.1*len(t)) #parametro para que fitee solo la parte final
+    k = np.argmin(np.abs(t - 0.2e-6))
+
+    temporal_taus = list()
+    incert_estad = 0
+
+    for nro_rep in range(NRO_REPS):
+        poblacion_fit = RealRho22[k:].copy()
+        if ADD_NOISE:
+            # Genero ruido similar al medido
+            ruido_poisson = poisson.rvs(NOISE_VAR, size=len(poblacion_fit))
+            # Reescaleo el ruido para que sea compatible con las simulaciones
+            ruido_poisson *= np.max(poblacion_fit)/NOISE_ADJ
+
+            poblacion_fit += ruido_poisson
+            # poblacion_fit += norm.rvs(0, NOISE_VAR, size=len(poblacion_fit))*np.max(poblacion_fit)/NOISE_ADJ
+            poblacion_fit[np.where(poblacion_fit < 0 )] = 0 # Limpio posibles valores negativos
+
+        t_fit = t[k+1:]
+
+        popt, pcov = curve_fit(expo, t_fit, poblacion_fit, p0=(1e-6, 1e-2, 1e-2))
+
+        temporal_taus.append(popt[0])
+        incert_estad += np.sqrt(pcov[0,0])
+        # temporal_N0.append(popt[1])
+        if not ADD_NOISE: # si no quiere que agregue ruido, entonces solo hace el for una vez
+            break
+
+        print(f" {nro_rep}  ", end='\r')
+
+    if ADD_NOISE:
+        allmeans.append(np.mean(temporal_taus))
+        allstds.append(np.std(temporal_taus) + incert_estad/NRO_REPS)
+    else:
+        allmeans.append(np.mean(temporal_taus))
+
+
+    #ploteo en funcion de 100*Rho22 para que me devuelva en porcentaje la población del excitado que es proporcional a la fluorescencia detectada
+    # lbl='Detuning Rep1:' + str(DetRep1) + ' MHz'
+    # ax1.plot(t[:-1]*1e6, RealRho22, '.', label=lbl, markersize=1)
+    # ax1.plot(t[:-1]*1e6, 100*RealRho22, '.', label=lbl, markersize=1)
+    # ax1.plot(t_fit*1e6, 100*expo(t_fit, *popt))
+    print(f'({i:02d}/{len(SatParDopplerVec)}) Done {SatParDoppler}'.ljust(50, ' '))
+
+# ax1.legend(SatParDopplerVec)
+# ax1.set_xlabel('Time (us)')
+# ax1.set_ylabel('100*Rho22')
+
+#%% #####################################################################################
+### Guardo valores al CSV para analizar luego
+# nDet=f"Det{np.abs(DetDoppler)}"
+# CSV_FNAME = f"error_teorico_tau_{nDet}.csv"
+# try:
+#     import pandas as pd
+#     data = pd.read_csv(CSV_FNAME)
+#     data["meanval"] = allmeans
+#     data["stdval"] = allstds
+#     data.to_csv(CSV_FNAME, index=False)
+#     print(f"SAVED {nDet}".ljust(50, ' '))
+# except FileNotFoundError:
+#     data = pd.DataFrame({'I': laser_intensities, "meanval": allmeans, "stdval": allstds})
+#     data.to_csv(CSV_FNAME, index=False)
+#     print(f"CREATED CSV FILE".ljust(50, ' '))
+#     print(f"SAVED DETUNING {-DetDoppler}".ljust(50, ' '))
+
+#%% ####################################################################################
+### Plots de los valores de intensidades medidas para distintos parametros
+### Cada axes arma un eje horizontal distinto, es para comparar las posibilidades
+
+# fig4, ax4 = plt.subplots()
+# allmeans = [t*1e6 for t in allmeans]
+
+# ax4b = ax4.twinx()
+# ax4.plot(SatParDopplerVec, allmeans, '-o', lw=0.4)
+# ax4b.plot(SatParDopplerVec, allN0, 'k-^', lw=0.4)
+# ax4.set_xlabel("Parametro de saturación")
+# ax4.set_ylabel("Tau (circulo)")
+# ax4b.set_ylabel("Alturas (triang)")
+# ax4.set_title(f"SP; IR_Lasers=0; DetuningDoppler={DetDoppler}; LineWidth={DopplerLaserLinewidth}")
+
+# ax4by = ax4.twiny()
+# ax4by.plot(laser_intensities, allmeans, '-', lw=0)
+# ax4by.xaxis.set_label_position('bottom')
+# ax4by.xaxis.tick_bottom()
+# ax4by.spines['bottom'].set_position(("axes", -0.27))
+# ax4by.set_xlabel(r"Intensidad [ $\mu$W/$\mu$m$^2$ ]")
+
+# ax4cy = ax4.twiny()
+# ax4cy.plot(laser_intensities*np.pi*(35**2), allmeans, '-', lw=0)
+# ax4cy.xaxis.set_label_position('bottom')
+# ax4cy.xaxis.tick_bottom()
+# ax4cy.spines['bottom'].set_position(("axes", -0.58))
+# ax4cy.set_xlabel(r"Potencia (r=35$\mu$m) [ $\mu$W ]")
+
+# ax4dy = ax4.twiny()
+# ax4dy.plot(laser_intensities*np.pi*(50**2), allmeans, '-', lw=0)
+# ax4dy.xaxis.set_label_position('bottom')
+# ax4dy.xaxis.tick_bottom()
+# ax4dy.spines['bottom'].set_position(("axes", -0.9))
+# ax4dy.set_xlabel(r"Potencia (r=50$\mu$m) [ $\mu$W ]")

analisis/plots/20221017_IonStatistics/Simulaciones/OBE3levels_functions.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May 15 18:11:00 2020
+
+@author: oem
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from odeintw import odeintw
+import time
+#matplotlib.use('Qt5Agg')
+
+
+"""
+odeintw es un wrapper de la funcion odeint de scipy.integrate que sirve para integrar
+sistemas de ecuaciones con parámetros complejos.
+En este script se resuelven las ecuaciones
+de Bloch de un sistema de 3 niveles en configuración lambda (1 y 3 son ground, 2 es excitado)
+para ion de Bario o de Calcio.
+Este script considera que hay dos láseres repump prendidos con sus respectivas frecuencias de
+rabi y detuning.
+Además, se usan dos distintos sistemas de referencia rotantes, condensados en las funciones
+RepumpReferenceSystemModel 1 y 2. La 2 es más útil para prender de a poco el láser repump2, y
+cuando su frec. de rabi es cero, no aparece en las ecuaciones.
+https://github.com/WarrenWeckesser/odeintw
+
+
+"""
+
+
+#CHEQUEADOS LOS DAMPINGS
+
+
+def EquationSystemModel(z, t, rabi12, rabi23, g21, g23, glg, glr, detr, detg, no_IR_lasers):
+    """
+    Se plantean acá las ecuaciones dinámicas de la matriz densidad (OBE: Optical Bloch Equations)
+    para la interacción entre un átomo y un campo electromagnético.
+        -rabi12, rabi23: frecuencias de rabi, proporcionales a la intensidad de los láseres
+        -g21, g23: anchos de línea naturales de las transiciones
+        -glg, glr: anchos de línea de los láseres doppler (UV) y repump (IR) respectivamente
+        -detg, detr: detunings de transición doppler (UV) y repump (IR) respectivamente
+        -no_IR_lasers: número de láseres IR aplicados
+    """
+
+    rho33, rho22, rho11, rho12, rho13, rho23 = z
+
+    rho21 = np.conjugate(rho12)
+    rho31 = np.conjugate(rho13)
+    rho32 = np.conjugate(rho23)
+
+    drho11dt = (1j)*0.5*rabi12*(rho12 - rho21) + rho22*g21
+    drho22dt = (-1j)*0.5*rabi12*(rho12 - rho21) - 1j*0.5*(np.conjugate(rabi23)*rho32 - rabi23*rho23) - rho22*(g21 + g23)
+    drho33dt = (-1j)*0.5*(rabi23*rho23 - np.conjugate(rabi23)*rho32) + g23*rho22
+
+    drho12dt = (-1j)*0.5*rabi12*(rho22 - rho11) - (1j)*(detg*rho12 - 0.5*rabi23*rho13) - rho12*(0.5*g21 + 0.5*g23 + glg)
+    drho13dt = (-1j)*(detg*rho13 + 0.5*rabi12*rho23 - 0.5*np.conjugate(rabi23)*rho12 - detr*rho13) - rho13*(glg + no_IR_lasers*glr)
+    drho23dt = (-1j)*0.5*rabi12*rho13 - (1j)*0.5*np.conjugate(rabi23)*(rho33-rho22) + (1j)*detr*rho23 - rho23*(0.5*g21 + 0.5*g23 + no_IR_lasers*glr)
+
+    dzdt = [drho33dt, drho22dt, drho11dt, drho12dt, drho13dt, drho23dt]
+
+    return dzdt
+
+
+def RepumpReferenceSystemModel(rabi23repump1, rabi23repump2, Detuningrepump1, Detuningrepump2, t):
+    """
+    Cuando hay más de un láser en una transición se producen batidos entre los mismos, haciendo que
+    la frecuencia de rabi sea un parámetro que varía en el tiempo. Para resolver eso numéricamente se
+    implementa de esta manera, calculando la frecuencia de rabi efectiva antes de resolver las OBE.
+    Al resolver el sistema, se pasa a un marco de referencia rotante respecto a los láseres. Como los dos láseres
+    IR excitan la misma transición, hay que elegir una de ambas frecuencias para rotar. Esto elige al repump1
+    como tal referencia, pero es posible elegir a repump2, o incluso a un promedio entre ambas frecuencias (por eso se llama detuningmedio).
+    Es decir, esto considera a la matriz de rotación como:
+        U = exp(-i.wg.t)|1><1| + |2><2| + exp(-i.wr1.t)|3><3|
+    """
+
+    rabi23 = [rabi23repump1 + rabi23repump2*np.exp(-1*(Detuningrepump1-Detuningrepump2)*1j*tiemp) for tiemp in t]
+    DetuningMedioRepump = Detuningrepump1
+    return rabi23, DetuningMedioRepump
+
+
+def GetTransitionLinewidths(species):
+    """
+    Devuelve los anchos de línea naturales según la especie iónica
+    """
+
+
+    if species == 'Calcium':
+        g21, g23, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
+    elif species == 'Barium':
+        g21, g23, = 2*np.pi*15.1e6, 2*np.pi*5.3e6 #anchos de linea de las transiciones
+    else:
+        print('Species not implemented')
+        raise
+
+    return g21, g23
+
+
+def Solve3LevelBloch(EquationSystemModel, species, SatParDoppler, SatParRep1, SatParRep2, DetDoppler, DetRep1, DetRep2,
+                     glg, glr, no_IR_lasers, initcond=1, tfin=10, paso=1e-3, printprogress=False):
+    """
+    Resuelve las ecuaciones de bloch del modelo dado. Los inputs son:
+        glg, glr: anchos de línea de doppler y repump en MHz
+        SatPar: parámetro de saturación de cada láser
+        Det: detunings en MHz
+        rabi12, rabi23 = sr*g21, sp*g23 #frecuencias de rabi
+    """
+
+    t = np.arange(0, tfin*1e-6, paso*1e-6)
+
+    g21, g23 = GetTransitionLinewidths(species)
+
+
+    #convierto todo a MHz
+    DopplerLaserLinewidth = glg*2*np.pi*1e6
+    RepumpLaserLinewidth = glr*2*np.pi*1e6
+
+    DetuningDoppler = DetDoppler*2*np.pi*1e6
+    Detuningrepump1 = DetRep1*2*np.pi*1e6
+    Detuningrepump2 = DetRep2*2*np.pi*1e6
+
+    rabi12 = np.sqrt(SatParDoppler*(4*(DetuningDoppler**2) + g21**2)/2)
+    rabi23repump1 = np.sqrt(SatParRep1*(4*(Detuningrepump1**2) + g23**2)/2)
+    rabi23repump2 = np.sqrt(SatParRep2*(4*(Detuningrepump2**2) + g23**2)/2)
+
+    rabi23vst, _ = RepumpReferenceSystemModel(rabi23repump1, rabi23repump2, Detuningrepump1, Detuningrepump2, t)
+
+    rho33sol = []
+    rho22sol = []
+    rho11sol = []
+
+    rho12sol = []
+    rho13sol = []
+    rho23sol = []
+
+
+    #condicion inicial: poblacion en el nivel 1
+    if initcond==1:
+        z0 = [0+0j, 0+0j, 1+0j, 0+0j, 0+0j, 0+0j] #formato: [rho33, rho22, rho11, rho12, rho13, rho23]
+    elif initcond==2:
+        z0 = [0+0j, 1+0j, 0+0j, 0+0j, 0+0j, 0+0j]
+    elif initcond==3:
+        z0 = [1+0j, 0+0j, 0+0j, 0+0j, 0+0j, 0+0j]
+    else:
+        print('Mala condicion inicial')
+        raise
+
+    #z0 = [0.5+0j, 0+0j, 0.5+0j, -0.5+0j, 0+0j, 0+0j] DARK STATE
+
+    for i in range(1,len(t)):
+        # span for next time step
+        tspan = [t[i-1],t[i]]
+        # solve for next step
+        #z = odeintw(EquationSystemModel, z0, tspan, args=(rabi12, rabi23vst[i], g21, g23, DopplerLaserLinewidth, RepumpLaserLinewidth, DetuningMedioRepump, DetuningDoppler), full_output=True)
+        z = odeintw(EquationSystemModel, z0, tspan, args=(rabi12, rabi23repump1, g21, g23, DopplerLaserLinewidth, RepumpLaserLinewidth, Detuningrepump1, DetuningDoppler, no_IR_lasers), full_output=True)
+        # store solution for plotting
+        rho22sol.append(z[0][1][1])
+        rho33sol.append(z[0][1][0])
+        rho11sol.append(z[0][1][2])
+
+        rho12sol.append(z[0][1][3])
+        rho13sol.append(z[0][1][4])
+        rho23sol.append(z[0][1][5])
+
+        # next initial condition
+        z0 = z[0][1]
+        if printprogress:
+            print(i, '/', len(t), end="\r", flush=True)
+
+    #Tomo la parte real de los diagonales y el módulo de las coherencias
+    RealRho22 = [r.real for r in rho22sol]
+    RealRho11 = [r.real for r in rho11sol]
+    RealRho33 = [r.real for r in rho33sol]
+    AbsRho12 = [np.abs(r) for r in rho12sol]
+    AbsRho13 = [np.abs(r) for r in rho13sol]
+    AbsRho23 = [np.abs(r) for r in rho23sol]
+
+    return t, RealRho11, RealRho22, RealRho33, AbsRho12, AbsRho13, AbsRho23
+
+
+
+def RepumpSpectrum(EquationSystemModel, species, SatParDoppler, SatParRep1, SatParRep2, DetDoppler, DetRep2,
+                     glg, glr, tfinOBE=10, pasoOBE=1e-3, frecini=-100, frecfin=100, pasofrec=1, plotSpectrum=True, sweepingrepump=True):
+    """
+    Hace un barrido en la frecuencia del repump.
+    Considerando que puede haber oscilaciones de la población por haber dos repumps,
+    """
+
+    g21, g23 = GetTransitionLinewidths(species)
+
+    FluorescenceVector = []
+
+    t = np.arange(0, tfinOBE*1e-6, pasoOBE*1e-6)
+
+    DopplerLaserLinewidth = glg*2*np.pi*1e6
+    RepumpLaserLinewidth = glr*2*np.pi*1e6
+
+    DetuningDoppler = DetDoppler*2*np.pi*1e6
+    Detuningrepump2 = DetRep2*2*np.pi*1e6
+
+    rabi12 = np.sqrt(SatParDoppler*(4*(DetuningDoppler**2) + g21**2)/2)
+    rabi23repump1 = np.sqrt(SatParRep1*(4*(Detuningrepump1**2) + g23**2)/2)
+    rabi23repump2 = np.sqrt(SatParRep2*(4*(Detuningrepump2**2) + g23**2)/2)
+
+    Repump1DetuningVector = np.arange(frecini*2*np.pi*1e6, frecfin*2*np.pi*1e6, pasofrec*2*np.pi*1e6)
+
+    counter = 1
+
+    FluorescenceVector = []
+
+    for DetRep1 in Repump1DetuningVector:
+        tini = time.time()
+
+        rho33sol = []
+        rho22sol = []
+        rho11sol = []
+
+        rho12sol = []
+        rho13sol = []
+        rho23sol = []
+
+        rabi23vst, DetuningMedioRepump = RepumpReferenceSystemModel(rabi23repump1, rabi23repump2, DetRep1, Detuningrepump2, t)
+
+        z0 = [0+0j, 0+0j, 1+0j, 0+0j, 0+0j, 0+0j] #formato: [rho33, rho22, rho11, rho12, rho13, rho23]
+
+
+        for i in range(1,len(t)):
+            # span for next time step
+            tspan = [t[i-1],t[i]]
+            # solve for next step
+            z = odeintw(EquationSystemModel,z0,tspan, args=(rabi12, rabi23vst[i], g21, g23, DopplerLaserLinewidth, RepumpLaserLinewidth, DetuningMedioRepump, DetuningDoppler), full_output=True)
+            # store solution for plotting
+
+            rho22sol.append(z[0][1][1])
+            rho33sol.append(z[0][1][0])
+            rho11sol.append(z[0][1][2])
+
+            rho12sol.append(z[0][1][3])
+            rho13sol.append(z[0][1][4])
+            rho23sol.append(z[0][1][5])
+
+            # next initial condition
+            z0 = z[0][1]
+    #        print(i, '/', len(t))
+        print
+
+        FluorescenceVector.append(np.real(rho12sol[-1]))
+
+        tfin = time.time()
+        print(counter, '/', len(Repump1DetuningVector), ', elapsed time: ', round(tfin-tini, 2), ' s')
+        counter = counter + 1
+
+    if plotSpectrum:
+        plt.plot([1e-6*d/(2*np.pi) for d in Repump1DetuningVector], [100*f for f in FluorescenceVector], 'o')
+        plt.xlabel('Repump detuning (MHz')
+        plt.ylabel('Fluorescence (A.U.)')
+        plt.title('Doppler detuning: ' + str(DetDoppler) + ' MHz, Repump 2 detuning: ' + str(DetRep2) + ' MHz')
+
+    return Repump1DetuningVector, FluorescenceVector
+
+
+#Funciones auxiliares para el análisis
+def FindFluoAtTime(fluo, time, timetarget, pasoti):
+    """
+    Encuentra la fluorescencia para determinado tiempo (en us).
+    El pasoti es la tolerancia temporal para el cual busca esa fluorescencia
+    """
+    ti = timetarget*1e-6
+    pasoti = pasoti*1e-6
+    i = 0
+    fluoattivec = []
+    while i < len(fluo):
+        if abs(time[i]-ti) < 2*pasoti:
+            fluoattivec.append(fluo[i])
+        i = i + 1
+    return np.mean(fluoattivec)
+
+def FindTimeForFluo(fluo, time, fluotarget, tolerancefluo=1e-3):
+    i = 0
+    timesforfluo = []
+    while i < len(time)-1:
+        if np.abs(fluo[i]-fluotarget) < 2*tolerancefluo:
+            timesforfluo.append(time[i])
+        i = i + 1
+    if len(timesforfluo) == 0:
+        print('La fluorescencia no toma ese valor')
+        return 0
+    elif len(timesforfluo) == 1:
+        return timesforfluo[0]
+    else:
+#        return np.mean(timesforfluo[1:])
+        return timesforfluo[-1]
+
+def FindCharacteristicDecayTime(fluo, time):
+    pass
+
+
+def IntegrateWindow(fluo, time, ti, window):
+    ti = ti*1e-6
+    window = window*1e-6
+    i = 0
+    fluowindow = []
+    while i<len(fluo):
+        if time[i]>ti and time[i]<=ti+window:
+            fluowindow.append(fluo[i])
+        i = i + 1
+        if time[i]>ti+window:
+            return np.mean(fluowindow)
+
+
+def Lorentz(x, x0, gamma, ymax):
+    scale = ymax*np.pi*gamma*0.5
+    return (scale/np.pi)*(0.5*gamma)/((x-x0)**2 + (0.5*gamma)**2)
+
+def rho23model(rabip, gamma, gammas, detuningp, rho0, t):
+    rho23 = ((1j)*rho0*rabip/(2*gamma))*(np.exp(-gamma*t) - np.exp(-(gammas+(rabip**2)/(4*gamma))*t -(1j)*detuningp*t))
+    return [np.real(rho23), np.abs(np.imag(rho23))]
+
+
+def dopplerBroadening(T, wldopp=0.397e-6, wlrep=0.866e-6, alpha=0, mcalcio = 6.655e-23*1e-3):
+    """
+    Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
+    wlr 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/(wldopp*wldopp)) + (1/(wlrep*wlrep)) - 2*(1/(wldopp*wlrep))*np.cos(alpha))*np.sqrt(kboltzmann*T/(2*mcalcio))
+
+    return gammaD
+

analisis/plots/20221017_IonStatistics/Simulaciones/error_teorico_tau_det20.csv

+I,meanval,stdval
+0.00662085,2.63563456542463e-07,3.62795752e-08
+0.00582507,2.674966392578803e-07,3.99478821e-08
+0.0049338,2.7426723419290064e-07,3.66249293e-08
+0.00397887,2.8263714050526124e-07,3.83554920e-08
+0.00305577,2.9770905441813384e-07,3.90291265e-08
+0.00226,3.1786383950851434e-07,4.27893803e-08
+0.00155972,3.562791436115684e-07,4.67339215e-08
+0.00101859,4.198979464288562e-07,5.28900825e-08
+0.00060479,5.457916989936301e-07,6.95843326e-08
+0.00031831,8.275611779916097e-07,1.30369576e-07

analisis/plots/20221017_IonStatistics/Simulaciones/medidos_powers_taus.csv

+Pow,Tau,N0
+208.0,4.107136373675053e-07,242.3949275128121
+183.0,3.994370315889412e-07,255.700355179914
+155.0,4.4882562779198994e-07,226.84461233128172
+125.0,4.727871415748306e-07,217.84774878119222
+96.0,5.254285404397582e-07,222.57732196795308
+71.0,6.079817946589069e-07,189.5426628192835
+49.0,7.274218641328625e-07,168.69185328968453
+32.0,8.823248220314772e-07,140.40570782757854
+19.0,1.280491601444429e-06,90.5188484395718
+10.0,1.96934320559994e-06,60.092743064494314

analisis/plots/20221017_IonStatistics/Simulaciones/plot_decaytimes_vs_intensity_detunings.py

+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+#plt.ion()
+
+dataREAL = pd.read_csv("medidos_powers_taus.csv")
+dataSIM  = pd.read_csv("sim_intensidad_taus_full.csv")
+errsSIM  = pd.read_csv("error_teorico_tau_det20.csv")
+I_sat = 43.3*1e-5 #uW/um2
+
+#%% #####################################################################################
+### Plot de lineas teoricas tau-intensidad variando detuning + datos medidos
+
+fig2, ax2 = plt.subplots()
+detuning_lines = ['Det10', 'Det20', 'Det30', 'Det40', 'Det50', 'Det60', 'Det70', 'Det80', 'Det90']
+
+# Curvas de las simulaciones con tau en unidades de microseg
+for det in detuning_lines:
+    ax2.plot(dataSIM.I, dataSIM[det]*1e6, 'k-', ms=1, lw=0.5)
+
+# Curvas de las mediciones cambiando el radio posible con tau en microseg
+
+for cambio in [0]: # esto es para shiftear la potencia medida
+    potencias = dataREAL.Pow - cambio
+    for rad in [85.7]: # esto evalua con distintos radios
+        ax2.errorbar([p/(np.pi*(rad**2)) for p in UVpotVec], Taus, yerr=np.mean(1e6*errsSIM.stdval.values),
+                    fmt='.--', ms=6, lw=.5, \
+                    color="#3949ab",
+                    capsize=3)
+    
+    
+        #ax2.errorbar(potencias/(np.pi*(rad**2)), dataREAL.Tau*1e6, \
+        #             yerr = 1e6*errsSIM.stdval.values, \
+        #             fmt='.--', ms=6, lw=.5, \
+        #             color="#3949ab",
+        #             capsize=3,
+        #             label=f'r={rad}; cambio={cambio}')
+
+ax2.set_xlim([1.2e-4, 0.2e-1])
+ax2.set_ylim([2e-1, 0.8e1])
+ax2.set_ylabel(r"$\tau$ [$\mu$s]")
+#ax2.set_xlabel(r"$I = P / (\pi\,r^2)$ [$\mu$W/$\mu$m$^2$]")
+ax2.set_xlabel(r'UV intensity ($\mu$W/$\mu$m$^2$)')
+# ax2.set_title(r"datos(punteadas) $r \in [30; 190]\ \mu m$   |   simulacion(llenas) $\Delta \in -[10; 90]$ MHz")
+#ax2.legend(markerscale=2)
+ax2.grid(which='minor', alpha=0.2)
+ax2.set_xscale("log")
+ax2.set_yscale("log")
+
+