Merge branch 'master' of https://code.df.uba.ar/nnunez/artiq_experiments

analisis/plots/20230424_CPT_TemperatureSens/CPT_plotter_20230424.py

@@ -14,12 +14,15 @@ from scipy import interpolate
 
 os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20230424_CPT_TemperatureSens/Data/')
 
+#os.chdir('/home/muri/nubeDF/Documents/codigos/artiq_experiments/analisis/plots/20230424_CPT_TemperatureSens/Data/')
+
 CPT_FILES = """000011345-IR_Scan_withcal_optimized
 000011331-IR_Scan_withcal_optimized
 """ 
 
-HEATING_FILES = """000011346-IR_Scan_withcal_optimized
-000011347-IR_Scan_withcal_optimized
+
+HEATING_FILES = """000015318-IR_Scan_withcal_optimized
+000015319-IR_Scan_withcal_optimized
 000011348-IR_Scan_withcal_optimized
 000011349-IR_Scan_withcal_optimized
 000011350-IR_Scan_withcal_optimized
@@ -32,6 +35,13 @@ HEATING_FILES = """000011346-IR_Scan_withcal_optimized
 000011344-IR_Scan_withcal_optimized
 """
 
+CPT_FILES = ''
+
+for i in range(0,20):
+    CPT_FILES = CPT_FILES + f'0000153{39 + i}-IR_Scan_withcal_optimized\n'
+CPT_FILES = CPT_FILES
+
+
 
 def SeeKeys(files):
     for i, fname in enumerate(files.split()):
@@ -74,6 +84,8 @@ UVCPTAmp_heating = []
 No_measures_heating = []
 Freqs_heating = []
 
+#%%
+
 for i, fname in enumerate(HEATING_FILES.split()):
     print(str(i) + ' - ' + fname)
     #print(fname)
@@ -95,12 +107,14 @@ for i, fname in enumerate(HEATING_FILES.split()):
 Ploteo las cpt de referencia / plotting the reference CPT
 """
 
-jvec = [0,1]
+jvec = [0,1,2,3,4,5,6,7]
+jvec = np.arange(0,20)
 
 plt.figure()
 i = 0
 for j in jvec:
-    plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2)
+    #plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j], yerr=np.sqrt(Counts[j]), fmt='o', capsize=2, markersize=2)
+    plt.errorbar([2*f*1e-6 for f in Freqs[j]], Counts[j])
     i = i + 1
 plt.xlabel('Frecuencia (MHz)')
 plt.ylabel('counts')

analisis/plots/20230804_RotationalDopplerShift_v2/RDS_CPT2.py

@@ -262,7 +262,7 @@ secular=False
 
 
 plt.plot([2*f*1e-6 for f in ExtraIR1_Freqs[1][1:]], [1e-3*c for c in ExtraCounts[1][1:]], '-o', markersize=2, label='Colineales')
-plt.xlabel('Frequency (MHz)', fontname='STIXgeneral',fontsize=15)
+plt.xlabel('IR1 Frequency (MHz)', fontname='STIXgeneral',fontsize=15)
 plt.ylabel('kCounts',fontname='STIXgeneral',fontsize=15)
 plt.xlim(425, 455)
 plt.yticks([1,1.5,2,2.5,3],fontname='STIXgeneral',fontsize=15)

analisis/plots/20230912_RotationalDopplerShift_v6/RDS_piezobeamsizes.py

@@ -46,8 +46,8 @@ def SeeKeys(files):
         print(fname)
         print(list(data['datasets'].keys()))
 
-def Lorentzian( x, A, B, x0, gam,C):
-    #C=0
+def Lorentzian( x, A, B, x0, gam):
+    C=0
     return A * gam**2 / ( gam**2 + ( x - x0 )**2) + B - C*(x - x0)
 
 
@@ -598,9 +598,12 @@ Resonancias DD configuracion +2/-2 colineal variando la ubicacion del ion en los
 Alta estadistica
 """
 
+def Lorentzian( x, A, B, x0, gam,C):
 
+    return A * gam**2 / ( gam**2 + ( x - x0 )**2) + B - C*(x - x0)
 
-palette = sns.color_palette("tab10")
+
+palette = sns.color_palette("rocket")
 
 pmlocmedvec = np.arange(0,len(HS_FILES),1)
 
@@ -644,17 +647,17 @@ for med in pmlocmedvec:
         
     fi=775
     if med not in [800]:
-        plt.plot([1e3*(2*f*1e-6-435)-fi for f in PiezoHSFrequencies[med][1:]], [1e-3*c for c in PiezoHSCounts[med][1:]], '-o', markersize=2, alpha=0.7)
-        plt.plot([1e3*(f-435)-fi for f in Freqs],[l*1e-3 for l in Lorentzian(Freqs,*popt)],linewidth=5)
+        plt.plot([1e3*(2*f*1e-6-435)-fi for f in PiezoHSFrequencies[med][1:]], [1e-3*c for c in PiezoHSCounts[med][1:]], '-o', markersize=2, alpha=0.7,color=palette[0])
+        plt.plot([1e3*(f-435)-fi for f in Freqs],[l*1e-3 for l in Lorentzian(Freqs,*popt)],linewidth=5,color=palette[0])
     
     jj=jj+1
-plt.xlabel('Frequency (kHz)',fontname='STIXgeneral',fontsize=50)
-plt.ylabel('kCounts',fontname='STIXgeneral',fontsize=50)
+plt.xlabel('IR1 Frequency (kHz)',fontname='STIXgeneral',fontsize=40)
+plt.ylabel('kCounts',fontname='STIXgeneral',fontsize=40)
 plt.xticks([-50,-25,0,25,50],fontname='STIXgeneral',fontsize=35)
 plt.yticks([3.5,4,4.5],fontname='STIXgeneral',fontsize=35)
 plt.tight_layout()
 plt.grid()
 plt.savefig('/home/nico/Nextcloud/Nico/Doctorado/Charlas/2023 Europe/DDresonancesexperimental_fine.pdf')
 # plt.legend()
-#plt.title(f'Ancho: {round(1e3*popt[3],2)} kHz')
+print(f'Ancho: {round(1e3*popt[3],2)} kHz')
 

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/EITfit/MM_eightLevel_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
+import MM_eightLevel_2repumps_python_scripts
+from MM_eightLevel_2repumps_python_scripts import CPTspectrum8levels,CPTspectrum8levels_vel
+from scipy.optimize import curve_fit
+import random
+from scipy.signal import savgol_filter as sf
+from scipy.stats import norm
+
+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):
+    """
+    solvemode=1: resuelve con np.linalg.solve
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    kboltz = 1.380649e-23
+    m = 6.6551079e-26
+    kg = np.pi * 2/397 *1e9
+    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)
+        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)
+            
+            FluorescencesVel.append(Fluorescence)    
+        FluorescencesVel = np.array(FluorescencesVel)
+        MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
+        Fluorescence = np.trapz(FluorescencesVel*MBprobMat,velvec,axis = 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()
+        
+
+    else:
+        drivefreq =2*np.pi* 1e3
+        FluorescencesVel = []
+        sigmaT = np.sqrt(T*kboltz/m)
+        velvec = np.linspace(-4*sigmaT,4*sigmaT,100)
+        MBprobVec = norm.pdf(velvec,loc = 0, scale = sigmaT)
+        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=None)
+            FluorescencesVel.append(Fluorescence)
+        FluorescencesVel = np.array(FluorescencesVel)
+        MBprobMat = np.tile(MBprobVec,(FluorescencesVel.shape[1],1)).T
+        Fluorescence = np.trapz(FluorescencesVel*MBprobMat,velvec,axis = 0)
+        ProbeDetuningVectorL = Frequencies
+    
+    #print('Done, Total time: ', round((tfinal-tinicial), 2), "s")   
+          
+    return ProbeDetuningVectorL, Fluorescence
+
+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
+    solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
+    """
+    tinicial = time.time()
+    ProbeDetuningVectorL, Fluovector = CPTspectrum8levels_vel(velvect,titavect,phivect,velprob,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(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(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_vel(velvect,titavect,phivect,velprob,sg, sp, gPS, gPD, DetDoppler, u, DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, drivefreq,beta, freqMin, freqMax, freqStep, circularityprobe=1, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
+    Frequencyvector, Fluovector = 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 = circularityprobe, plot=False, solvemode=1, detpvec=None)
+    NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
+    return Frequencyvector, NoisyFluovector
+
+
+
+def SmoothNoisyCPT(Fluo, window=11, poly=3):
+    SmoothenFluo = sf(Fluo, window, poly)
+    return SmoothenFluo
+
+
+
+def fitCPT_8levels(Freq,Fluo,sg,sp,gPS,gPD,DetDoppler,u,DopplerLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe,beta, drivefreq, freqMin, freqMax, freqStep, circularityprobe=1, plot=False,dephasing = False, boltzmann = False ):
+
+    def SpectrumForFit(Freq,sg,sp,T,DetDoppler):
+        freq,spectra = PerformExperiment_8levels(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=Freq)
+        return spectra
+
+    
+    popt,pcov = curve_fit(SpectrumForFit,Freq,Fluo,p0 = [sg,sp,T,DetDoppler],bounds = ([0.01,0.01,0.00001,-50e6],[1,20,1,30e6]))
+    fitted_spectra = SpectrumForFit(Freq,*popt)
+    return Freq, fitted_spectra,popt,pcov
+
+

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/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
+from MM_eightLevel_2repumps_AnalysisFunctions import PerformExperiment_8levels, GenerateNoisyCPT_vel, SmoothNoisyCPT, PerformExperiment_8levels_vel
+import matplotlib.pyplot as plt
+import time
+import h5py
+
+def MB_prob(vel,T):
+    return (m/(2*np.pi*kboltz*T))**(3/2)*2* vel**2 *np.exp(-vel**2 * m/(2*kboltz*T))
+
+
+
+def MB_prob_1d(vel,T):
+    return np.sqrt(m/(2*np.pi*kboltz*T))*np.exp(-vel**2 * m/(2*kboltz*T))
+
+#plt.rcParams.update({
+#    "text.usetex": True,
+#    "font.family": "CM Roman"
+#})
+
+
+
+plt.rcParams["axes.prop_cycle"] = plt.cycler('color', ['#DBAD1F', '#C213DB', '#DB4F2A', '#0500DB', '#09DB9B', '#B4001B', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf'])
+
+
+ub = 9.27e-24 #magneton de bohr
+h = 6.63e-34 #cte de planck
+c = (ub/h)*1e-4 #en unidades de MHz/G
+kboltz = 1.380649e-23
+
+m = 6.6551079e-26 
+
+u = 32e6 #esto equivale aprox al campo que tenemos
+
+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.1 #ancho de linea de los laseres en MHz
+DopplerLaserLinewidth, ProbeLaserLinewidth = lw, lw #ancho de linea de los laseres
+
+DetDoppler = -15 #detuning doppler en MHz
+
+T = 0.1e-3 #temperatura en K
+alpha = 0 #angulo entre los láseres
+
+#estos son los angulos de la polarizacion de los laseres respecto al campo magnetico
+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 
+
+#Parametros de la simulacion cpt todo en MHz
+center = -15
+span = 30
+freqMin = center-span*0.5
+freqMax = center+span*0.5
+freqStep = 0.5
+
+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.3
+sp = 3
+
+drivefreq=2*np.pi*22.135*1e6
+
+kg = np.pi * 2/397 *1e9
+kp = np.pi * 2/866 *1e9
+
+#betavec = np.arange(0,1.1,0.1)
+betavec=[0]
+
+
+
+#fig1, ax1 = plt.subplots()
+
+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.28]
+s23vec = np.linspace(0.5,15,20)
+#s23vec = [s23vec[2],s23vec[6],s23vec[10],s23vec[18]]
+s23vec = [10]
+
+Tmat = np.zeros((len(s12vec),len(s23vec)))
+
+Tvec_ajuste = np.zeros(len(Tvec))
+phivect = 0
+
+T = 10e-3
+
+titavect = np.linspace(0,2*np.pi,100)    
+  
+betag = 0
+betap = 0 
+alpha = 0
+
+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)
+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)
+
+#%%
+plt.plot(Freq,Fluo_boltz,label = 'Velocities')
+plt.plot(Freq,Fluo_sb,label = 'Sideband')
+plt.plot(Freq,Fluo_deph, label = 'Dephasing')
+
+plt.legend()
+
+#%%
+
+os.chdir('/home/muri/nubeDF/Documents/codigos/artiq_experiments/analisis/plots/20230920_CPT_TemperatureSens_v2/Data/')
+
+CPT_FILES = """000011345-IR_Scan_withcal_optimized
+000011331-IR_Scan_withcal_optimized
+""" 
+
+for i in range(0,9):
+    CPT_FILES = CPT_FILES + f'0000153{19 + i}-IR_Scan_withcal_optimized\n'
+CPT_FILES = CPT_FILES
+
+
+
+def SeeKeys(files):
+    for i, fname in enumerate(files.split()):
+        data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+        print(fname)
+        print(list(data['datasets'].keys()))
+
+print(SeeKeys(CPT_FILES))
+
+Counts = []
+Freqs = []
+
+AmpTisa = []
+UVCPTAmp = []
+No_measures = []
+
+for i, fname in enumerate(CPT_FILES.split()):
+    print(str(i) + ' - ' + fname)
+    #print(fname)
+    data = h5py.File(fname+'.h5', 'r') # Leo el h5: Recordar que nuestros datos estan en 'datasets'
+
+    # Aca hago algo repugnante para poder levantar los strings que dejamos
+    # que además tenian un error de tipeo al final. Esto no deberá ser necesario
+    # cuando se solucione el error este del guardado.
+    Freqs.append(np.array(data['datasets']['IR1_Frequencies']))
+    Counts.append(np.array(data['datasets']['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']))
+
+

analisis/plots/20230920_CPT_TemperatureSens_v2/Data/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]))