Commit 4d32f5a2 authored by Nicolas Nunez Barreto's avatar Nicolas Nunez Barreto

peds

parent 371bb08a
......@@ -28,9 +28,9 @@
fit-margin-left="0"
fit-margin-right="0"
fit-margin-bottom="0"
inkscape:zoom="1.4142136"
inkscape:cx="2526.139"
inkscape:cy="-697.56084"
inkscape:zoom="1"
inkscape:cx="2455.5"
inkscape:cy="-763.5"
inkscape:window-width="1846"
inkscape:window-height="1016"
inkscape:window-x="0"
......@@ -2324,13 +2324,6 @@ YAmFwWDYy/8DHTJV9KTZCpkAAAAASUVORK5CYII=
cy="-365.55286"
rx="2.9858103"
ry="2.675787" />
<ellipse
style="fill:#aa0000;stroke:#aa0000;stroke-width:2;stroke-linejoin:round"
id="path880-7-1-76-3"
cx="350.41907"
cy="-365.55286"
rx="2.9858103"
ry="2.675787" />
<ellipse
style="fill:#aa0000;stroke:#aa0000;stroke-width:2;stroke-linejoin:round"
id="path880-7-1-76-8"
......@@ -2345,27 +2338,6 @@ YAmFwWDYy/8DHTJV9KTZCpkAAAAASUVORK5CYII=
cy="-365.88306"
rx="2.9858103"
ry="2.675787" />
<ellipse
style="fill:#aa0000;stroke:#aa0000;stroke-width:2;stroke-linejoin:round"
id="path880-7-1-76-8-85"
cx="391.52667"
cy="-365.7836"
rx="2.9858103"
ry="2.675787" />
<ellipse
style="fill:#aa0000;stroke:#aa0000;stroke-width:2;stroke-linejoin:round"
id="path880-7-1-76-8-0"
cx="379.66742"
cy="-365.55286"
rx="2.9858103"
ry="2.675787" />
<ellipse
style="fill:#aa0000;stroke:#aa0000;stroke-width:2;stroke-linejoin:round"
id="path880-7-1-76-1"
cx="364.76254"
cy="-365.55286"
rx="2.9858103"
ry="2.675787" />
<ellipse
style="fill:#aa0000;stroke:#aa0000;stroke-width:2;stroke-linejoin:round"
id="path880-7-1-76-2"
......@@ -2380,13 +2352,6 @@ YAmFwWDYy/8DHTJV9KTZCpkAAAAASUVORK5CYII=
cy="-365.55286"
rx="2.9858103"
ry="2.675787" />
<ellipse
style="fill:#aa0000;stroke:#aa0000;stroke-width:2;stroke-linejoin:round"
id="path880-7-1-76-9"
cx="338.44809"
cy="-365.55286"
rx="2.9858103"
ry="2.675787" />
<path
style="fill:none;stroke:#000000;stroke-width:1.465;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#Arrow1Lend-2-1)"
d="m 293.72037,-415.09093 c 33.9602,-18.35849 69.11882,-0.53749 69.11882,-0.53749"
......
......@@ -283,11 +283,14 @@ def sensgauss(angle):
radius = np.arange(0.1e-6,10e-6,0.01e-6)
ang = 1
plt.figure()
plt.semilogy([r*1e6 for r in radius],sens(radius), linewidth=3, color='coral', label='D-D CPT +2/-2 OAM ')
plt.xlabel('OAM radius (um)')
plt.ylabel('Thermometry sensitivity')
plt.axhline(8.6e6, color='indigo', linestyle='--', label='S-P CPT gaussian sensitivity')
plt.axhline(8.6e6, color='indigo', linestyle='--', label='S-P CPT gaussian sensitivity')
plt.grid()
plt.legend()
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 2 16:30:09 2020
@author: oem
"""
import os
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.signal import argrelextrema
from EITfit.threeLevel_2repumps_linealpol_python_scripts import CPTspectrum8levels, CPTspectrum8levels_fixedRabi
import random
from scipy.signal import savgol_filter as sf
def CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump):
if titadoppler==0:
NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
elif titadoppler==90:
NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
else:
NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
return [detuningdoppler + dr for dr in NegativeDR], [detuningrepump + dr for dr in PositiveDR]
def GetClosestIndex(Vector, value, tolerance=1e-3):
i = 0
while i<len(Vector):
if abs(Vector[i] - value) < tolerance:
return i
else:
i = i + 1
return GetClosestIndex(Vector, value, tolerance=2*tolerance)
def FindDRFrequencies(Freq, Fluo, TeoDR, entorno=3):
"""
Busca los indices y la frecuencia de los minimos en un entorno cercano al de la DR.
Si no encuentra, devuelve el valor teórico.
"""
IndiceDRteo1, IndiceEntornoinicialDRteo1, IndiceEntornofinalDRteo1 = GetClosestIndex(Freq, TeoDR[0]), GetClosestIndex(Freq, TeoDR[0]-entorno), GetClosestIndex(Freq, TeoDR[0]+entorno)
IndiceDRteo2, IndiceEntornoinicialDRteo2, IndiceEntornofinalDRteo2 = GetClosestIndex(Freq, TeoDR[1]), GetClosestIndex(Freq, TeoDR[1]-entorno), GetClosestIndex(Freq, TeoDR[1]+entorno)
IndiceDRteo3, IndiceEntornoinicialDRteo3, IndiceEntornofinalDRteo3 = GetClosestIndex(Freq, TeoDR[2]), GetClosestIndex(Freq, TeoDR[2]-entorno), GetClosestIndex(Freq, TeoDR[2]+entorno)
IndiceDRteo4, IndiceEntornoinicialDRteo4, IndiceEntornofinalDRteo4 = GetClosestIndex(Freq, TeoDR[3]), GetClosestIndex(Freq, TeoDR[3]-entorno), GetClosestIndex(Freq, TeoDR[3]+entorno)
IndiceDRteo5, IndiceEntornoinicialDRteo5, IndiceEntornofinalDRteo5 = GetClosestIndex(Freq, TeoDR[4]), GetClosestIndex(Freq, TeoDR[4]-entorno), GetClosestIndex(Freq, TeoDR[4]+entorno)
IndiceDRteo6, IndiceEntornoinicialDRteo6, IndiceEntornofinalDRteo6 = GetClosestIndex(Freq, TeoDR[5]), GetClosestIndex(Freq, TeoDR[5]-entorno), GetClosestIndex(Freq, TeoDR[5]+entorno)
EntornoFreqDR1, EntornoFreqDR2 = Freq[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Freq[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
EntornoFreqDR3, EntornoFreqDR4 = Freq[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Freq[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
EntornoFreqDR5, EntornoFreqDR6 = Freq[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Freq[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
EntornoFluoDR1, EntornoFluoDR2 = Fluo[IndiceEntornoinicialDRteo1:IndiceEntornofinalDRteo1], Fluo[IndiceEntornoinicialDRteo2:IndiceEntornofinalDRteo2]
EntornoFluoDR3, EntornoFluoDR4 = Fluo[IndiceEntornoinicialDRteo3:IndiceEntornofinalDRteo3], Fluo[IndiceEntornoinicialDRteo4:IndiceEntornofinalDRteo4]
EntornoFluoDR5, EntornoFluoDR6 = Fluo[IndiceEntornoinicialDRteo5:IndiceEntornofinalDRteo5], Fluo[IndiceEntornoinicialDRteo6:IndiceEntornofinalDRteo6]
IndiceFluoMinimaEntorno1, IndiceFluoMinimaEntorno2 = argrelextrema(np.array(EntornoFluoDR1), np.less)[0], argrelextrema(np.array(EntornoFluoDR2), np.less)[0]
IndiceFluoMinimaEntorno3, IndiceFluoMinimaEntorno4 = argrelextrema(np.array(EntornoFluoDR3), np.less)[0], argrelextrema(np.array(EntornoFluoDR4), np.less)[0]
IndiceFluoMinimaEntorno5, IndiceFluoMinimaEntorno6 = argrelextrema(np.array(EntornoFluoDR5), np.less)[0], argrelextrema(np.array(EntornoFluoDR6), np.less)[0]
try:
FreqDR1 = EntornoFreqDR1[int(IndiceFluoMinimaEntorno1)]
IndiceDR1 = GetClosestIndex(Freq, FreqDR1)
except:
FreqDR1 = TeoDR[0]
IndiceDR1 = IndiceDRteo1
try:
FreqDR2 = EntornoFreqDR2[int(IndiceFluoMinimaEntorno2)]
IndiceDR2 = GetClosestIndex(Freq, FreqDR2)
except:
FreqDR2 = TeoDR[1]
IndiceDR2 = IndiceDRteo2
try:
FreqDR3 = EntornoFreqDR3[int(IndiceFluoMinimaEntorno3)]
IndiceDR3 = GetClosestIndex(Freq, FreqDR3)
except:
FreqDR3 = TeoDR[2]
IndiceDR3 = IndiceDRteo3
try:
FreqDR4 = EntornoFreqDR4[int(IndiceFluoMinimaEntorno4)]
IndiceDR4 = GetClosestIndex(Freq, FreqDR4)
except:
FreqDR4 = TeoDR[3]
IndiceDR4 = IndiceDRteo4
try:
FreqDR5 = EntornoFreqDR5[int(IndiceFluoMinimaEntorno5)]
IndiceDR5 = GetClosestIndex(Freq, FreqDR5)
except:
FreqDR5 = TeoDR[4]
IndiceDR5 = IndiceDRteo5
try:
FreqDR6 = EntornoFreqDR6[int(IndiceFluoMinimaEntorno6)]
IndiceDR6 = GetClosestIndex(Freq, FreqDR6)
except:
FreqDR6 = TeoDR[5]
IndiceDR6 = IndiceDRteo6
return [IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6], [FreqDR1, FreqDR2, FreqDR3, FreqDR4, FreqDR5, FreqDR6]
def FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=-100, getindices=False):
"""
Toma los indices donde estan las DR y evalua su fluorescencia. Esos indices son minimos locales en un entorno
cercano a las DR teoricas y, si no hay ningun minimo, toma la teorica.
Luego, hace el cociente de esa fluorescencia y un factor de normalización segun NormalizationCriterium:
1: Devuelve la fluorescencia absoluta de los minimos
2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no
hubiera una resonancia oscura y hubiera una recta. Ese valor esta a DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
o el valor limite por izquierda de la curva
4: Deuelve el cociente entre la fluorescencia del minimo y el valor de fluorescencia a detuning 0 MHz
"""
IndiceDR1, IndiceDR2, IndiceDR3, IndiceDR4, IndiceDR5, IndiceDR6 = IndicesDR[0], IndicesDR[1], IndicesDR[2], IndicesDR[3], IndicesDR[4], IndicesDR[5]
FluorescenceOfMinimums = [Fluo[IndiceDR1], Fluo[IndiceDR2], Fluo[IndiceDR3], Fluo[IndiceDR4], Fluo[IndiceDR5], Fluo[IndiceDR6]]
FrequencyOfMinimums = [Freq[IndiceDR1], Freq[IndiceDR2], Freq[IndiceDR3], Freq[IndiceDR4], Freq[IndiceDR5], Freq[IndiceDR6]]
DistanciaFrecuenciaCociente = 25
if NormalizationCriterium==0:
print('che')
return FrequencyOfMinimums, FluorescenceOfMinimums
if NormalizationCriterium==1:
Fluorescenciacerodetuning = Fluo[GetClosestIndex(Freq, 0)]
Fluorescenciaasintotica = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
return FrequencyOfMinimums, np.array([Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica, Fluorescenciacerodetuning/Fluorescenciaasintotica])
if NormalizationCriterium==2:
k = 0
while k < len(Freq):
if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
FluoIzquierda = Fluo[k]
indiceizquierda = k
print('Izq:', Freq[k])
break
else:
k = k + 1
l = 0
while l < len(Freq):
if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
FluoDerecha = Fluo[l]
indicederecha = l
print('Der: ', Freq[l])
break
else:
l = l + 1
FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
print(FluoNormDivisor)
if NormalizationCriterium==3:
#asintotico
FluoNormDivisor = Fluo[GetClosestIndex(Freq, frecuenciareferenciacriterioasintotico)]
if NormalizationCriterium==4:
#este te tira la fluorescencia de detuning 0
FluoNormDivisor = Fluo[GetClosestIndex(Freq, 0)]
RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
print('Esto: ', RelativeFluorescenceOfMinimums)
if NormalizationCriterium==2 and getindices==True:
return FrequencyOfMinimums, RelativeFluorescenceOfMinimums, indiceizquierda, indicederecha
return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
def GetFinalMaps(MapasDR1, MapasDR2, MapasDR3, MapasDR4, MapasDR5, MapasDR6):
"""
Nota: esto vale para polarizacion del 397 sigma+ + sigma-. Sino hay que cambiar los coeficientes.
La estructura es:
MapasDRi = [MapaMedido_criterio1_DRi, MapaMedido_criterio2_DRi, MapaMedido_criterio3_DRi, MapaMedido_criterio4_DRi]
"""
Mapa1 = MapasDR1[0]
Mapa2pi = np.sqrt(3)*(MapasDR2[1] + MapasDR5[1])
Mapa2smas = np.sqrt(12/2)*MapasDR3[1] + (2/np.sqrt(2))*MapasDR6[1]
Mapa2smenos = (2/np.sqrt(2))*MapasDR1[1] + np.sqrt(12/2)*MapasDR4[1]
Mapa3pi = np.sqrt(3)*(MapasDR2[2] + MapasDR5[2])
Mapa3smas = np.sqrt(12/2)*MapasDR3[2] + (2/np.sqrt(2))*MapasDR6[2]
Mapa3smenos = (2/np.sqrt(2))*MapasDR1[2] + np.sqrt(12/2)*MapasDR4[2]
return Mapa1, [Mapa2pi, Mapa2smas, Mapa2smenos], [Mapa3pi, Mapa3smas, Mapa3smenos]
def CombinateDRwithCG(RelMinMedido1, RelMinMedido2, RelMinMedido3, RelMinMedido4):
Fluo1 = RelMinMedido1[0]
Fluo2pi = np.sqrt(3)*(RelMinMedido2[1] + RelMinMedido2[4])
Fluo2smas = np.sqrt(12/2)*RelMinMedido2[2] + (2/np.sqrt(2))*RelMinMedido2[5]
Fluo2smenos = (2/np.sqrt(2))*RelMinMedido2[0] + np.sqrt(12/2)*RelMinMedido2[3]
Fluo3pi = np.sqrt(3)*(RelMinMedido3[1] + RelMinMedido3[4])
Fluo3smas = np.sqrt(12/2)*RelMinMedido3[2] + (2/np.sqrt(2))*RelMinMedido3[5]
Fluo3smenos = (2/np.sqrt(2))*RelMinMedido3[0] + np.sqrt(12/2)*RelMinMedido3[3]
return Fluo1, [Fluo2pi, Fluo2smas, Fluo2smenos], [Fluo3pi, Fluo3smas, Fluo3smenos]
def IdentifyPolarizationCoincidences(theoricalmap, target, tolerance=1e-1):
"""
Busca en un mapa 2D la presencia de un valor target (medido) con tolerancia tolerance.
Si lo encuentra, pone un 1. Sino, un 0. Al plotear con pcolor se verá
en blanco la zona donde el valor medido se puede hallar.
"""
CoincidenceMatrix = np.zeros((len(theoricalmap), len(theoricalmap[0])))
i = 0
while i<len(theoricalmap):
j = 0
while j<len(theoricalmap[0]):
if abs(theoricalmap[i][j]-target) < tolerance:
CoincidenceMatrix[i][j] = 1
j=j+1
i=i+1
return CoincidenceMatrix
def RetrieveAbsoluteCoincidencesBetweenMaps(MapsVectors):
MatrixSum = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
AbsoluteCoincidencesMatrix = np.zeros((len(MapsVectors[0]), len(MapsVectors[0][0])))
MatrixMapsVectors = []
for i in range(len(MapsVectors)):
MatrixMapsVectors.append(np.matrix(MapsVectors[i]))
for i in range(len(MatrixMapsVectors)):
MatrixSum = MatrixSum + MatrixMapsVectors[i]
MaxNumberOfCoincidences = np.max(MatrixSum)
ListMatrixSum = [list(i) for i in list(np.array(MatrixSum))]
for i in range(len(ListMatrixSum)):
for j in range(len(ListMatrixSum[0])):
if ListMatrixSum[i][j] == MaxNumberOfCoincidences:
AbsoluteCoincidencesMatrix[i][j] = 1
return AbsoluteCoincidencesMatrix, MaxNumberOfCoincidences
def MeasureMeanValueOfEstimatedArea(AbsoluteCoincidencesMap, X, Y):
NonZeroIndices = np.nonzero(AbsoluteCoincidencesMap)
Xsum = 0
Xvec = []
Ysum = 0
Yvec = []
N = len(NonZeroIndices[0])
for i in range(N):
Xsum = Xsum + X[NonZeroIndices[1][i]]
Xvec.append(X[NonZeroIndices[1][i]])
Ysum = Ysum + Y[NonZeroIndices[0][i]]
Yvec.append(Y[NonZeroIndices[0][i]])
Xaverage = Xsum/N
Yaverage = Ysum/N
Xspread = np.std(Xvec)
Yspread = np.std(Yvec)
return Xaverage, Yaverage, N, Xspread, Yspread
def MeasureRelativeFluorescenceFromCPT(Freq, Fluo, u, titadoppler, detuningrepump, detuningdoppler, frefasint=-100, entorno=3):
ResonanciasTeoricas, ResonanciasPositivas = CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)
IndicesDR, FreqsDR = FindDRFrequencies(Freq, Fluo, ResonanciasTeoricas, entorno=entorno)
FrequencyOfMinimums, RelativeFluorescenceOfMinimums0 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=0, frecuenciareferenciacriterioasintotico=frefasint)
FrequencyOfMinimums, RelativeFluorescenceOfMinimums1 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=1, frecuenciareferenciacriterioasintotico=frefasint)
FrequencyOfMinimums, RelativeFluorescenceOfMinimums2, indiceizquierda, indicederecha = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=2, frecuenciareferenciacriterioasintotico=frefasint, getindices=True)
FrequencyOfMinimums, RelativeFluorescenceOfMinimums3 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=3, frecuenciareferenciacriterioasintotico=frefasint)
FrequencyOfMinimums, RelativeFluorescenceOfMinimums4 = FindRelativeFluorescencesOfDR(Freq, Fluo, IndicesDR, detuningdoppler, NormalizationCriterium=4, frecuenciareferenciacriterioasintotico=frefasint)
print('hola')
print(RelativeFluorescenceOfMinimums0)
return RelativeFluorescenceOfMinimums0, RelativeFluorescenceOfMinimums1, RelativeFluorescenceOfMinimums2, RelativeFluorescenceOfMinimums3, RelativeFluorescenceOfMinimums4, IndicesDR, [indiceizquierda, indicederecha]
def GenerateNoisyCPT(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
Frequencyvector, Fluovector = PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
return Frequencyvector, NoisyFluovector
def GenerateNoisyCPT_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None)
NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
return Frequencyvector, NoisyFluovector
def GenerateNoisyCPT_fit(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqs, plot=False, solvemode=1, detpvec=None, noiseamplitude=0.001):
Frequencyvector, Fluovector = PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, min(freqs), max(freqs) + freqs[1]-freqs[0], freqs[1]-freqs[0], plot=False, solvemode=1, detpvec=None)
NoisyFluovector = [fluo+noiseamplitude*(2*random.random()-1) for fluo in Fluovector]
return Frequencyvector, NoisyFluovector
def AddNoiseToCPT(Fluo, noisefactor):
return [f+noisefactor*(2*random.random()-1) for f in Fluo]
def SmoothNoisyCPT(Fluo, window=11, poly=3):
SmoothenFluo = sf(Fluo, window, poly)
return SmoothenFluo
def GetMinimaInfo(Freq, Fluo, u, titadoppler, detuningdoppler, detuningrepump, MinimumCriterium=2, NormalizationCriterium=1):
"""
FUNCION VIEJA
Esta funcion devuelve valores de frecuencias y fluorescencia relativa de los minimos.
Minimumcriterion:
1: Saca los minimos con funcion argelextrema
2: Directamente con las frecuencias teoricas busca las fluorescencias
Normalizationcriterium:
1: Devuelve la fluorescencia absoluta de los minimos
2: Devuelve el cociente entre la fluorescencia del minimo y un valor medio entre dos puntos lejanos, como si no
hubiera una resonancia oscura y hubiera una recta. Ese valor esta a DistanciaFrecuenciaCociente del detuning del azul (el punto medio entre las dos DR en este caso)
3: Devuelve el cociente entre la fluorescencia del minimo y el valor a -100 MHz (si se hizo de -100 a 100),
o el valor limite por izquierda de la curva
"""
FluorescenceOfMaximum = max(Fluo)
FrequencyOfMaximum = Freq[Fluo.index(FluorescenceOfMaximum)]
#criterio para encontrar los minimos
#criterio usando minimos de la fluorescencia calculados con la curva
if MinimumCriterium == 1:
LocationOfMinimums = argrelextrema(np.array(Fluo), np.less)[0]
FluorescenceOfMinimums = np.array([Fluo[i] for i in LocationOfMinimums])
FrequencyOfMinimums = np.array([Freq[j] for j in LocationOfMinimums])
#criterio con las DR teoricas
if MinimumCriterium == 2:
FrecuenciasDRTeoricas, FrecuenciasDRTeoricasPositivas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
FrequencyOfMinimums = []
FluorescenceOfMinimums =[]
print(FrecuenciasDRTeoricas)
k=0
ventanita = 0.001
while k < len(Freq):
if Freq[k] < FrecuenciasDRTeoricas[0] + ventanita and Freq[k] > FrecuenciasDRTeoricas[0] - ventanita:
FrequencyOfMinimums.append(Freq[k])
FluorescenceOfMinimums.append(Fluo[k])
elif Freq[k] < FrecuenciasDRTeoricas[1] + ventanita and Freq[k] > FrecuenciasDRTeoricas[1] - ventanita:
FrequencyOfMinimums.append(Freq[k])
FluorescenceOfMinimums.append(Fluo[k])
elif Freq[k] < FrecuenciasDRTeoricas[2] + ventanita and Freq[k] > FrecuenciasDRTeoricas[2] - ventanita:
FrequencyOfMinimums.append(Freq[k])
FluorescenceOfMinimums.append(Fluo[k])
elif Freq[k] < FrecuenciasDRTeoricas[3] + ventanita and Freq[k] > FrecuenciasDRTeoricas[3] - ventanita:
FrequencyOfMinimums.append(Freq[k])
FluorescenceOfMinimums.append(Fluo[k])
elif Freq[k] < FrecuenciasDRTeoricas[4] + ventanita and Freq[k] > FrecuenciasDRTeoricas[4] - ventanita:
FrequencyOfMinimums.append(Freq[k])
FluorescenceOfMinimums.append(Fluo[k])
elif Freq[k] < FrecuenciasDRTeoricas[5] + ventanita and Freq[k] > FrecuenciasDRTeoricas[5] - ventanita:
FrequencyOfMinimums.append(Freq[k])
FluorescenceOfMinimums.append(Fluo[k])
k = k + 1
print(FrequencyOfMinimums)
if len(FrequencyOfMinimums) != len(FrecuenciasDRTeoricas):
print('NO ANDA BIEN ESTO PAPI, revisalo')
#esto es para establecer un criterio para la fluorescencia relativa
DistanciaFrecuenciaCociente = 15
if NormalizationCriterium==1:
FluoNormDivisor = 1
if NormalizationCriterium==2:
k = 0
while k < len(Freq):
if Freq[k] < detuningdoppler-DistanciaFrecuenciaCociente + 2 and Freq[k] > detuningdoppler-DistanciaFrecuenciaCociente - 2:
FluoIzquierda = Fluo[k]
print('Izq:', Freq[k])
break
else:
k = k + 1
l = 0
while l < len(Freq):
if Freq[l] < detuningdoppler+DistanciaFrecuenciaCociente + 2 and Freq[l] > detuningdoppler+DistanciaFrecuenciaCociente - 2:
FluoDerecha = Fluo[l]
print('Der: ', Freq[l])
break
else:
l = l + 1
FluoNormDivisor = 0.5*(FluoDerecha+FluoIzquierda)
print(FluoNormDivisor)
if NormalizationCriterium==3:
FluoNormDivisor = Fluo[0]
RelativeFluorescenceOfMinimums = np.array([Fluore/FluoNormDivisor for Fluore in FluorescenceOfMinimums])
return FrequencyOfMinimums, RelativeFluorescenceOfMinimums
def GetPlotsofFluovsAngle_8levels(FrequencyOfMinimumsVector, RelativeFluorescenceOfMinimumsVector, u, titadoppler, detuningdoppler, detuningrepump, ventana=0.25, taketheoricalDR=False):
#primero buscamos las frecuencias referencia que se parezcan a las 6:
i = 0
FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[0]
FrecuenciasDRTeoricas = [darkresonance for darkresonance in CalculoTeoricoDarkResonances_8levels(u, titadoppler, detuningdoppler, detuningrepump)[0]]
while i < len(FrequencyOfMinimumsVector):
if len(FrequencyOfMinimumsVector[i])==len(FrecuenciasDRTeoricas):
FrecuenciasReferenciaBase = FrequencyOfMinimumsVector[i]
print('Cool! Taking the DR identified with any curve')
break
else:
i = i + 1
if i==len(FrequencyOfMinimumsVector):
print('No hay ningun plot con 5 resonancias oscuras. Tomo las teóricas')
FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
if taketheoricalDR:
FrecuenciasReferenciaBase = FrecuenciasDRTeoricas
Ventana = abs(ventana*(FrecuenciasReferenciaBase[1] - FrecuenciasReferenciaBase[0])) #ventana separadora de resonancias
print('Ventana = ', Ventana)
DarkResonance1Frequency = []
DarkResonance1Fluorescence = []
DarkResonance2Frequency = []
DarkResonance2Fluorescence = []
DarkResonance3Frequency = []
DarkResonance3Fluorescence = []
DarkResonance4Frequency = []
DarkResonance4Fluorescence = []
DarkResonance5Frequency = []
DarkResonance5Fluorescence = []
DarkResonance6Frequency = []
DarkResonance6Fluorescence = []
i = 0
while i < len(FrequencyOfMinimumsVector):
j = 0
FrecuenciasReferencia = [i for i in FrecuenciasReferenciaBase]
while j < len(FrequencyOfMinimumsVector[i]):
if abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[0])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[0])-Ventana):
DarkResonance1Frequency.append(FrequencyOfMinimumsVector[i][j])
DarkResonance1Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
FrecuenciasReferencia[0] = 0
elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[1])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[1])-Ventana):
DarkResonance2Frequency.append(FrequencyOfMinimumsVector[i][j])
DarkResonance2Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
FrecuenciasReferencia[1] = 0
elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[2])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[2])-Ventana):
DarkResonance3Frequency.append(FrequencyOfMinimumsVector[i][j])
DarkResonance3Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
FrecuenciasReferencia[2] = 0
elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[3])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[3])-Ventana):
DarkResonance4Frequency.append(FrequencyOfMinimumsVector[i][j])
DarkResonance4Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
FrecuenciasReferencia[3] = 0
elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[4])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[4])-Ventana):
DarkResonance5Frequency.append(FrequencyOfMinimumsVector[i][j])
DarkResonance5Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
FrecuenciasReferencia[4] = 0
elif abs(FrequencyOfMinimumsVector[i][j]) < (abs(FrecuenciasReferencia[5])+Ventana) and abs(FrequencyOfMinimumsVector[i][j]) >= (abs(FrecuenciasReferencia[5])-Ventana):
DarkResonance6Frequency.append(FrequencyOfMinimumsVector[i][j])
DarkResonance6Fluorescence.append(RelativeFluorescenceOfMinimumsVector[i][j])
FrecuenciasReferencia[5] = 0
else:
#print('Algo anduvo mal, por ahi tenes que cambiar la ventana che')
pass
j = j + 1
if np.count_nonzero(FrecuenciasReferencia) > 0:
if FrecuenciasReferencia[0] != 0:
DarkResonance1Frequency.append(FrecuenciasReferencia[0])
DarkResonance1Fluorescence.append()
if FrecuenciasReferencia[1] != 0:
DarkResonance2Frequency.append(FrecuenciasReferencia[1])
DarkResonance2Fluorescence.append(0)
if FrecuenciasReferencia[2] != 0:
DarkResonance3Frequency.append(FrecuenciasReferencia[2])
DarkResonance3Fluorescence.append(0)
if FrecuenciasReferencia[3] != 0:
DarkResonance4Frequency.append(FrecuenciasReferencia[3])
DarkResonance4Fluorescence.append(0)
if FrecuenciasReferencia[4] != 0:
DarkResonance5Frequency.append(FrecuenciasReferencia[4])
DarkResonance5Fluorescence.append(0)
if FrecuenciasReferencia[5] != 0:
DarkResonance6Frequency.append(FrecuenciasReferencia[5])
DarkResonance6Fluorescence.append(0)
i = i + 1
return DarkResonance1Frequency, DarkResonance1Fluorescence, DarkResonance2Frequency, DarkResonance2Fluorescence, DarkResonance3Frequency, DarkResonance3Fluorescence, DarkResonance4Frequency, DarkResonance4Fluorescence, DarkResonance5Frequency, DarkResonance5Fluorescence, DarkResonance6Frequency, DarkResonance6Fluorescence, FrecuenciasReferenciaBase
def PerformExperiment_8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
"""
Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
solvemode=1: resuelve con np.linalg.solve
solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
"""
Fluovectors = []
for titaprobe in titaprobeVec:
tinicial = time.time()
ProbeDetuningVectorL, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
tfinal = time.time()
print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")
if plot:
plt.figure()
plt.xlabel('Repump detuning (MHz')
plt.ylabel('Fluorescence (A.U.)')
plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
plt.legend()
Fluovectors.append(Fluovector)
if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
Fluovectors = Fluovector
return ProbeDetuningVectorL, Fluovectors
def PerformExperiment_8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobeVec, phirepump, titarepump, freqMin, freqMax, freqStep, plot=False, solvemode=1, detpvec=None):
"""
Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
solvemode=1: resuelve con np.linalg.solve
solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
"""
Fluovectors = []
for titaprobe in titaprobeVec:
tinicial = time.time()
ProbeDetuningVectorL, Fluovector = CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, DetDoppler, DetRepump, u, DopplerLaserLinewidth, RepumpLaserLinewidth, ProbeLaserLinewidth, T, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=False, solvemode=1)
tfinal = time.time()
print('Done angle ', titarepump, ' Total time: ', round((tfinal-tinicial), 2), "s")
if plot:
plt.figure()
plt.xlabel('Repump detuning (MHz')
plt.ylabel('Fluorescence (A.U.)')
plt.plot(ProbeDetuningVectorL, Fluovector, label=str(titarepump)+'º tita repump, T: ' + str(T*1e3) + ' mK')
plt.legend()
Fluovectors.append(Fluovector)
if len(titaprobeVec) == 1: #esto es para que no devuelva un vector de vectores si solo fijamos un angulo
Fluovectors = Fluovector
return ProbeDetuningVectorL, Fluovectors
#!/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]))
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 7 22:30:01 2020
@author: nico
"""
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.signal import argrelextrema
"""
Scripts para el calculo de la curva CPT
"""
def H0matrix(Detg, Detp, u):
"""
Calcula la matriz H0 en donde dr es el detuning del doppler, dp es el retuning del repump y u es el campo magnético en Hz/Gauss.
Para esto se toma la energía del nivel P como 0
"""
eigenEnergies = (Detg-u, Detg+u, -u/3, u/3, Detp-6*u/5, Detp-2*u/5, Detp+2*u/5, Detp+6*u/5) #pagina 26 de Oberst. los lande del calcio son iguales a Bario.
H0 = np.diag(eigenEnergies)
return H0
def HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe):
"""
Calcula la matriz de interacción Hsp + Hpd, en donde rabR es la frecuencia de rabi de la transición Doppler SP,
rabP es la frecuencia de rabi de la transición repump DP, y las componentes ei_r y ei_p son las componentes de la polarización
del campo eléctrico incidente de doppler y repump respectivamente. Deben estar normalizadas a 1
"""
HI = np.zeros((8, 8), dtype=np.complex_)
i, j = 1, 3
HI[i-1, j-1] = (rabG/np.sqrt(3)) * np.cos(titadoppler)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 1, 4
HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(1j*phidoppler)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 2, 3
HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.sin(titadoppler)*np.exp(-1j*phidoppler)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 2, 4
HI[i-1, j-1] = -(rabG/np.sqrt(3)) * np.cos(titadoppler)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 3, 5
HI[i-1, j-1] = -(rabP/2) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 3, 6
HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 3, 7
HI[i-1, j-1] = rabP/np.sqrt(12) * np.sin(titaprobe)*np.exp(1j*phiprobe)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 4, 6
HI[i-1, j-1] = -(rabP/np.sqrt(12)) * np.sin(titaprobe)*np.exp(-1j*phiprobe)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 4, 7
HI[i-1, j-1] = -(rabP/np.sqrt(3)) * np.cos(titaprobe)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
i, j = 4, 8
HI[i-1, j-1] = (rabP/2) * np.sin(titaprobe)*np.exp(1j*phiprobe)
HI[j-1, i-1] = np.conjugate(HI[i-1, j-1])
return HI
def Lplusminus(detr, detp, phirepump, titarepump, forma=1):
Hintplus = np.zeros((8, 8), dtype=np.complex_)
Hintminus = np.zeros((8, 8), dtype=np.complex_)
Hintplus[4, 2] = (-1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
Hintplus[5, 2] = (-1/np.sqrt(3))*np.cos(titarepump)
Hintplus[6, 2] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
Hintplus[5, 3] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
Hintplus[6, 3] = (-1/np.sqrt(3))*np.cos(titarepump)
Hintplus[7, 3] = (1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
Hintminus[2, 4] = (-1/2)*np.sin(titarepump)*np.exp(-1j*phirepump)
Hintminus[2, 5] = (-1/np.sqrt(3))*np.cos(titarepump)
Hintminus[2, 6] = (1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(1j*phirepump)
Hintminus[3, 5] = (-1/(2*np.sqrt(3)))*np.sin(titarepump)*np.exp(-1j*phirepump)
Hintminus[3, 6] = (-1/np.sqrt(3))*np.cos(titarepump)
Hintminus[3, 7] = (1/2)*np.sin(titarepump)*np.exp(1j*phirepump)
if forma==1:
Lplus = np.zeros((64, 64), dtype=np.complex_)
Lminus = np.zeros((64, 64), dtype=np.complex_)
DeltaBar = np.zeros((64, 64), dtype=np.complex_)
for r in range(8):
for q in range(8):
for k in range(8):
for j in range(8):
if j==q:
if (k==2 or k==3) and r > 3:
Lplus[r*8+q][k*8+j] = (-1j)*(Hintplus[r,k])
if (r==2 or r==3) and k > 3:
Lminus[r*8+q][k*8+j] = (-1j)*(Hintminus[r,k])
elif r==k:
if (q==2 or q==3) and j > 3:
Lplus[r*8+q][k*8+j] = (-1j)*(- Hintplus[j,q])
if (j==2 or j==3) and q > 3:
Lminus[r*8+q][k*8+j] = (-1j)*(- Hintminus[j,q])
if forma==2:
deltaKro = np.diag([1, 1, 1, 1, 1, 1, 1, 1])
Lplus = (-1j)*(np.kron(Hintplus, deltaKro) - np.kron(deltaKro, Hintplus))
Lminus = (-1j)*(np.kron(Hintminus, deltaKro) - np.kron(deltaKro, Hintminus))
DeltaBar = np.zeros((64, 64), dtype=np.complex_)
for i in range(64):
DeltaBar[i, i] = (1j)*(detr - detp)
return np.matrix(Lminus), np.matrix(Lplus), np.matrix(DeltaBar)
def GetL1(Lplus, Lminus, DeltaBar, L0, rabR, nmax):
"""
Devuelve Splus0 y Sminus0
"""
Sp = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 - (nmax+1)*DeltaBar))*np.matrix(Lplus))
Sm = (-1)*(0.5*rabR)*(np.matrix(np.linalg.inv(L0 + (nmax+1)*DeltaBar))*np.matrix(Lminus))
for n in list(range(nmax+1))[(nmax+1)::-1][0:len(list(range(nmax+1))[(nmax+1)::-1])-1]: #jaja esto solo es para que vaya de nmax a 1 bajando. debe haber algo mas facil pero kcio
Sp = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 - n*DeltaBar + rabR*(Lminus*np.matrix(Sp))))*np.matrix(Lplus))
Sm = (-1)*(rabR)*(np.matrix(np.linalg.inv(L0 + n*DeltaBar + rabR*(Lplus*np.matrix(Sm))))*np.matrix(Lminus))
L1 = 0.5*rabR*(np.matrix(Lminus)*np.matrix(Sp) + np.matrix(Lplus)*np.matrix(Sm))
return L1
def EffectiveL(gPS, gPD, lwg, lwr, lwp):
"""
Siendo Heff = H + EffectiveL, calcula dicho EffectiveL que es (-0.5j)*sumatoria(CmDaga*Cm) que luego sirve para calcular el Liouvilliano
"""
Leff = np.zeros((8, 8), dtype=np.complex_)
Leff[0, 0] = 2*lwg
Leff[1, 1] = 2*lwg
Leff[2, 2] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
Leff[3, 3] = ((2/3)+(1/3))*gPS + ((1/2) + (1/6) + (1/3))*gPD
Leff[4, 4] = 2*(lwr + lwp)
Leff[5, 5] = 2*(lwr + lwp)
Leff[6, 6] = 2*(lwr + lwp)
Leff[7, 7] = 2*(lwr + lwp)
return (-0.5j)*Leff
def CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp):
"""
Si tomamos el Liuvilliano como L = (-j)*(Heff*deltak - Heffdaga*deltak) + sum(Mm),
esta funcion calcula dichos Mm, que tienen dimensión 64x64 ya que esa es la dimensión del L. Estas componentes
salen de hacer la cuenta a mano conociendo los Cm y considerando que Mm[8*(r-1)+s, 8*(k-1)+j] = Cm[r,l] + Cmdaga[j,s] = Cm[r,l] + Cm[s,j]
ya que los componentes de Cm son reales.
Esta M es la suma de las 8 matrices M.
"""
M = np.matrix(np.zeros((64, 64), dtype=np.complex_))
M[0,27] = (2/3)*gPS
M[9,18] = (2/3)*gPS
M[0,18] = (1/3)*gPS
M[1,19] = -(1/3)*gPS
M[8,26] = -(1/3)*gPS
M[9,27] = (1/3)*gPS
M[36,18] = (1/2)*gPD
M[37,19] = (1/np.sqrt(12))*gPD
M[44,26] = (1/np.sqrt(12))*gPD
M[45,27] = (1/6)*gPD
M[54,18] = (1/6)*gPD
M[55,19] = (1/np.sqrt(12))*gPD
M[62,26] = (1/np.sqrt(12))*gPD
M[63,27] = (1/2)*gPD
M[45,18] = (1/3)*gPD
M[46,19] = (1/3)*gPD
M[53,26] = (1/3)*gPD
M[54,27] = (1/3)*gPD
M[0,0] = 2*lwg
M[1,1] = 2*lwg
M[8,8] = 2*lwg
M[9,9] = 2*lwg
factor1 = 1
factor2 = 1
factor3 = 1
factor4 = 1
#M[36, 45] = lwp
M[36,36] = 2*(lwr + factor1*lwp)
M[37,37] = 2*(lwr + factor1*lwp)
M[38,38] = 2*(lwr + factor1*lwp)
M[39,39] = 2*(lwr + factor1*lwp)
M[44,44] = 2*(lwr + factor2*lwp)
M[45,45] = 2*(lwr + factor2*lwp)
M[46,46] = 2*(lwr + factor2*lwp)
M[47,47] = 2*(lwr + factor2*lwp)
M[52,52] = 2*(lwr + factor3*lwp)
M[53,53] = 2*(lwr + factor3*lwp)
M[54,54] = 2*(lwr + factor3*lwp)
M[55,55] = 2*(lwr + factor3*lwp)
M[60,60] = 2*(lwr + factor4*lwp)
M[61,61] = 2*(lwr + factor4*lwp)
M[62,62] = 2*(lwr + factor4*lwp)
M[63,63] = 2*(lwr + factor4*lwp)
return M
def dopplerBroadening(wlg, wlp, alpha, T, mcalcio = 6.655e-23*1e-3):
"""
Calcula el broadening extra semiclásico por temperatura considerando que el ion atrapado se mueve.
wlg es la longitud de onda doppler, wlp la longitud de onda repump, T la temperatura del ion en kelvin, y alpha (en rads) el ángulo
que forman ambos láseres.
"""
kboltzmann = 1.38e-23 #J/K
gammaD = (2*np.pi)*np.sqrt((1/(wlg*wlg)) + (1/(wlp*wlp)) - 2*(1/(wlg*wlp))*np.cos(alpha))*np.sqrt(kboltzmann*T/(2*mcalcio))
return gammaD
def FullL_efficient(rabG, rabR, rabP, gPS = 0, gPD = 0, Detg = 0, Detr = 0, Detp = 0, u = 0, lwg = 0, lwr=0, lwp = 0,
phidoppler=0, titadoppler=0, phiprobe=0, titaprobe=0, phirepump=0, titarepump=0, T = 0, alpha = 0):
"""
Calcula el Liouvilliano total de manera explícita índice a índice. Suma aparte las componentes de las matrices M.
Es la más eficiente hasta ahora.
"""
db = dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)
#lwr = np.sqrt(lwr**2 + dopplerBroadening(0.397e-6, 0.866e-6, alpha, T)**2)
lwg = np.sqrt(lwg**2 + db**2)
lwr = np.sqrt(lwr**2 + db**2)
CC = EffectiveL(gPS, gPD, lwg, lwr, lwp)
Heff = H0matrix(Detg, Detp, u) + HImatrix(rabG, rabP, phidoppler, titadoppler, phiprobe, titaprobe) + CC
Heffdaga = np.matrix(Heff).getH()
Lfullpartial = np.zeros((64, 64), dtype=np.complex_)
for r in range(8):
for q in range(8):
for k in range(8):
for j in range(8):
if j!=q and r!=k:
pass
elif j==q and r!=k:
if (r < 2 and k > 3) or (k < 2 and r > 3) or (r > 3 and k > 3) or (r==0 and k==1) or (r==1 and k==0) or (r==2 and k==3) or (r==3 and k==2): #todo esto sale de analizar explicitamente la matriz y tratar de no calcular cosas de más que dan cero
pass
else:
Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k])
elif j!=q and r==k:
if (j < 2 and q > 3) or (q < 2 and j > 3) or (j > 3 and q > 3) or (j==0 and q==1) or (j==1 and q==0) or (j==2 and q==3) or (j==3 and q==2):
pass
else:
Lfullpartial[r*8+q][k*8+j] = (-1j)*(-Heffdaga[j,q])
else:
if Heff[r,k] == Heffdaga[j,q]:
pass
else:
Lfullpartial[r*8+q][k*8+j] = (-1j)*(Heff[r,k]-Heffdaga[j,q])
M = CalculateSingleMmatrix(gPS, gPD, lwg, lwr, lwp)
L0 = np.array(np.matrix(Lfullpartial) + M)
nmax = 1
Lminus, Lplus, DeltaBar = Lplusminus(Detr, Detp, phirepump, titarepump)
factor1 = np.exp(1j*0.2*np.pi)
factor2 = np.exp(-1j*0.2*np.pi)
#print(factor)
L1 = GetL1(factor1*Lplus, factor2*Lminus, DeltaBar, L0, rabR, nmax)
Lfull = L0 + L1
#NORMALIZACION DE RHO
i = 0
while i < 64:
if i%9 == 0:
Lfull[0, i] = 1
else:
Lfull[0, i] = 0
i = i + 1
return Lfull
"""
Scripts para correr un experimento y hacer el análisis de los datos
"""
def CalculoTeoricoDarkResonances(u, titadoppler):
if titadoppler==0:
NegativeDR = [(-7/5)*u, (-3/5)*u, (-1/5)*u, (1/5)*u, (3/5)*u, (7/5)*u]
elif titadoppler==90:
NegativeDR = [(-11/5)*u, (-7/5)*u, (-3/5)*u, (3/5)*u, (7/5)*u, (11/5)*u]
PositiveDR = [(-8/5)*u, (-4/5)*u, 0, (4/5)*u, (8/5)*u]
return NegativeDR, PositiveDR
def CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump,
freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
"""
Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
solvemode=1: resuelve con np.linalg.solve
solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
"""
phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
#u = 2*np.pi*u*1e6
Fluovector = []
tinicial = time.time()
for Detp in DetProbeVector:
L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
if solvemode == 1:
rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)]))
Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
Fluovector.append(Fluo)
if solvemode == 2:
Linv = np.linalg.inv(L)
rhovectorized = [Linv[j][0] for j in range(len(Linv))]
Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
Fluovector.append(Fluo)
tfinal = time.time()
print('Done, Total time: ', round((tfinal-tinicial), 2), "s")
DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
if plot:
plt.xlabel('Probe detuning (MHz)')
plt.ylabel('Fluorescence (A.U.)')
plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
plt.legend()
return DetProbeVectorMHz, Fluovector
def CPTspectrum8levels_fixedRabi(sg, sr, sp, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump,
freqMin=-100, freqMax=100, freqStep=1e-1, plot=False, solvemode=1):
"""
Hace un experimento barriendo ángulos de repump con el angulo de doppler fijo.
solvemode=1: resuelve con np.linalg.solve
solvemode=2: resuelve invirtiendo L con la funcion np.linalg.inv
"""
phidoppler, titadoppler = phidoppler*(np.pi/180), titadoppler*(np.pi/180)
phiprobe, titaprobe = phiprobe*(np.pi/180), titaprobe*(np.pi/180)
phirepump, titarepump = phirepump*(np.pi/180), titarepump*(np.pi/180)
DetProbeVector = 2*np.pi*np.arange(freqMin*1e6, freqMax*1e6, freqStep*1e6)
Detg, Detr = 2*np.pi*Detg*1e6, 2*np.pi*Detr*1e6
#lwg, lwr, lwp = 2*np.pi*lwg*1e6, 2*np.pi*lwr*1e6, 2*np.pi*lwp*1e6
lwg, lwr, lwp = lwg*1e6, lwr*1e6, lwp*1e6
rabG = sg*gPS
rabR = sr*gPD
rabP = sp*gPD
#u = 2*np.pi*u*1e6
Fluovector = []
tinicial = time.time()
for Detp in DetProbeVector:
L = FullL_efficient(rabG, rabR, rabP, gPS, gPD, Detg, Detr, Detp, u, lwg, lwr, lwp, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, Temp, alpha)
if solvemode == 1:
coh = 5
rhovectorized = np.linalg.solve(L, np.array([int(i==0) for i in range(64)]))
#Fluo = np.abs(rhovectorized[coh])
Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
Fluovector.append(Fluo)
if solvemode == 2:
Linv = np.linalg.inv(L)
rhovectorized = [Linv[j][0] for j in range(len(Linv))]
Fluo = np.real(rhovectorized[18] + np.real(rhovectorized[27])) #estos son los rho33 + rho44
Fluovector.append(Fluo)
tfinal = time.time()
print('Done, Total time: ', round((tfinal-tinicial), 2), "s")
DetProbeVectorMHz = np.arange(freqMin, freqMax, freqStep)
if plot:
plt.xlabel('Probe detuning (MHz)')
plt.ylabel('Fluorescence (A.U.)')
plt.plot(DetProbeVectorMHz, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
plt.legend()
return DetProbeVectorMHz, Fluovector
#%%
if __name__ == "__main__":
ub = 9.27e-24
h = 6.63e-34
c = (ub/h)*1e-4 #en unidades de MHz/G
B = 25 #campo magnetico en gauss
u = c*B
sg, sr, sp = 0.5, 1.5, 4 #parámetros de saturación del doppler y repump
gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6 #anchos de linea de las transiciones
rabG, rabR, rabP = sg*gPS, sr*gPD, sp*gPD #frecuencias de rabi
lwg, lwr, lwp = 0.3, 0.3, 0.3 #ancho de linea de los laseres
Detg = -25
Detr = 20 #detuning del doppler y repump
Temp = 0.0e-3 #temperatura en K
alpha = 0*(np.pi/180) #angulo entre los láseres
phidoppler, titadoppler = 0, 90
phirepump, titarepump = 0, 90
phiprobe, titaprobe = 0, 90
plotCPT = False
freqMin = -50
freqMax = 50
freqStep = 5e-2
Frequencyvector, Fluovector = CPTspectrum8levels(rabG, rabR, rabP, gPS, gPD, Detg, Detr, u, lwg, lwr, lwp, Temp, alpha, phidoppler, titadoppler, phiprobe, titaprobe, phirepump, titarepump, freqMin=freqMin, freqMax=freqMax, freqStep=freqStep, plot=plotCPT, solvemode=1)
NegativeDR, PositiveDR = CalculoTeoricoDarkResonances(u/(2*np.pi*1e6), titadoppler)
plt.plot(Frequencyvector, [100*f for f in Fluovector], label=str(titaprobe) + 'º, T: ' + str(Temp*1e3) + ' mK')
plt.xlabel('Probe detuning (MHz)')
plt.ylabel('Fluorescence (A.U.)')
for PDR in PositiveDR:
plt.axvline(Detr+PDR, linestyle='--', linewidth=0.5, color='red')
for NDR in NegativeDR:
plt.axvline(Detg+NDR, linestyle='--', linewidth=0.5, color='blue')
#parametros que andan piola:
"""
ub = 9.27e-24
h = 6.63e-34
c = (ub/h)*1e-4 #en unidades de MHz/G
B = 17 #campo magnetico en gauss
u = c*B
#u = 80e6
sr, sp = 0.53, 4.2
gPS, gPD, = 2*np.pi*21.58e6, 2*np.pi*1.35e6
rabR, rabP = sr*gPS, sp*gPD
lw = 2*np.pi * 0.33e6
lwr, lwp = lw, lw #ancho de linea de los laseres
dr_spec = - 2*np.pi* 26e6
freqSteps = 500
freqMin = -100e6
freqMax = 100e6
dps = 2*np.pi*np.linspace(freqMin, freqMax, freqSteps)
#dps = [-30e6]
alfar = 90*(np.pi/180)
ex_r, ey_r, ez_r = np.sin(alfar)*np.cos(0), np.sin(alfar)*np.sin(0), np.cos(alfar)
alfap = 90*(np.pi/180)
ex_p, ey_p, ez_p = np.sin(alfap)*np.cos(0), np.sin(alfap)*np.sin(0), np.cos(alfap)
"""
import h5py
import matplotlib.pyplot as plt
import numpy as np
import sys
import re
import ast
from scipy.optimize import curve_fit
import os
from scipy import interpolate
#Mediciones barriendo angulo del TISA y viendo kicking de resonancias oscuras
#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20220106_CPT_DosLaseres_v08_TISA_DR\Data
os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20230817_RotationalDopplerShift_v5/Data')
"""
en este codigo ploteo espectros CPT de resonancias D-D para configuracion +2/+2 y +2/-2 (usando pentaprisma)
"""
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
LEFT_FILES = """VaryingBeamlocation/Left/000014728-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014729-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014730-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014731-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014732-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014733-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014735-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014737-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014739-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014740-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014741-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014742-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014743-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014744-IR_Scan_withcal_optimized
VaryingBeamlocation/Left/000014745-IR_Scan_withcal_optimized
"""
RIGHT_FILES = """VaryingBeamlocation/Right/000014725-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014749-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014750-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014751-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014752-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014755-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014757-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014758-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014760-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014761-IR_Scan_withcal_optimized
VaryingBeamlocation/Right/000014762-IR_Scan_withcal_optimized
"""
EXTRA_FILES = """VaryingBeamlocation/Extra/000014746-IR_Scan_withcal_optimized
VaryingBeamlocation/Extra/000014747-IR_Scan_withcal_optimized
VaryingBeamlocation/Extra/000014763-IR_Scan_withcal_optimized
VaryingBeamlocation/Extra/000014764-IR_Scan_withcal_optimized
"""
def Split(array,n):
length=len(array)/n
splitlist = []
jj = 0
while jj<length:
partial = []
ii = 0
while ii < n:
partial.append(array[jj*n+ii])
ii = ii + 1
splitlist.append(partial)
jj = jj + 1
return splitlist
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(LEFT_FILES))
print(SeeKeys(RIGHT_FILES))
print(SeeKeys(EXTRA_FILES))
LeftCounts = []
LeftFrequencies = []
RightCounts = []
RightFrequencies = []
ExtraCounts = []
ExtraFrequencies = []
for i, fname in enumerate(LEFT_FILES.split()):
print(str(i) + ' - ' + fname)
data = h5py.File(fname+'.h5', 'r')
LeftCounts.append(np.array(data['datasets']['counts_spectrum']))
LeftFrequencies.append(np.array(data['datasets']['IR1_Frequencies']))
for i, fname in enumerate(RIGHT_FILES.split()):
print(str(i) + ' - ' + fname)
data = h5py.File(fname+'.h5', 'r')
RightCounts.append(np.array(data['datasets']['counts_spectrum']))
RightFrequencies.append(np.array(data['datasets']['IR1_Frequencies']))
for i, fname in enumerate(EXTRA_FILES.split()):
print(str(i) + ' - ' + fname)
data = h5py.File(fname+'.h5', 'r')
ExtraCounts.append(np.array(data['datasets']['counts_spectrum']))
ExtraFrequencies.append(np.array(data['datasets']['IR1_Frequencies']))
def ErrorDRdepth(p, f, b):
ep = np.sqrt(p)
ef = np.sqrt(f)
eb = np.sqrt(b)
derivadap = 1/((f-b)**2)
derivadaf = ((p-b)/((f-b)**2))**2
derivadab = ((p-f)/((f-b)**2))**2
return 2*np.sqrt(derivadap*ep*ep + derivadaf*ef*ef + derivadab*eb*eb)
#%%
import seaborn as sns
"""
Resonancias DD configuracion +2/-2 colineal variando la ubicacion del ion en los haces
LOBULO LEFT DEL OAM
"""
palette = sns.color_palette("tab10")
#pmleftlocmedvec = [0,1,2,3,4,5,6,7,8,9,12]
#pmleftlocmedvec = [0,2,3,4,5,7,8,9,12] #esto anda bien
pmleftlocmedvec = [0,1,3,4,5,7,8,9,12]
"""
La 9, 10 y 11 son de la misma medicion. Hay que mergearlas
Lo mismo con la 12, 13 y 14.
Por ahora obvio esto pero despues hay que hacerlo.
"""
#pmlocmedvec = []
# idxvecdr1 = [125,125,125,125,125,125,125,125,125,125,125]
# idxvecdr2 = [406,406,406,406,406,406,406,406,406,406,406]
idxvecdr1 = [125,125,125,125,125,125,125,125,125]
idxvecdr2 = [406,406,406,406,406,406,406,406,406]
#pmlocmedvec = [2]
#LocVecsint = np.arange(1,6,1)
plt.figure()
bkg = 150
pmleftdepthsdr1=[]
errorpmleftdepthsdr1=[]
pmleftdepthsdr2=[]
errorpmleftdepthsdr2=[]
Intensityleft = []
errorIntensityleft = []
idxtest = 125
jj=0
for med in pmleftlocmedvec:
print(med)
if med == 1:
pmleftdepthsdr1.append(1-(LeftCounts[med][1:][idxvecdr1[jj]]-bkg)/(np.mean(LeftCounts[med][1:][70:90])-bkg))
pmleftdepthsdr2.append(1-(LeftCounts[med][1:][idxvecdr2[jj]]-bkg)/(np.mean(LeftCounts[med][1:][70:90])-bkg))
errorpmleftdepthsdr1.append(ErrorDRdepth(LeftCounts[med][1:][idxvecdr1[jj]],np.mean(LeftCounts[med][1:][70:90]), bkg))
errorpmleftdepthsdr2.append(ErrorDRdepth(LeftCounts[med][1:][idxvecdr2[jj]],np.mean(LeftCounts[med][1:][70:90]), bkg))
Intensityleft.append(np.mean(LeftCounts[med][1:][70:90])-bkg)
#errorIntensityleft.append(np.sqrt(np.mean(LeftCounts[med][1:][70:90]))+np.sqrt(bkg))
errorIntensityleft.append(2*np.std(LeftCounts[med][1:][70:90]))
else:
pmleftdepthsdr1.append(1-(LeftCounts[med][1:][idxvecdr1[jj]]-bkg)/(np.mean(LeftCounts[med][1:][0:20])-bkg))
pmleftdepthsdr2.append(1-(LeftCounts[med][1:][idxvecdr2[jj]]-bkg)/(np.mean(LeftCounts[med][1:][0:20])-bkg))
errorpmleftdepthsdr1.append(ErrorDRdepth(LeftCounts[med][1:][idxvecdr1[jj]],np.mean(LeftCounts[med][1:][0:20]), bkg))
errorpmleftdepthsdr2.append(ErrorDRdepth(LeftCounts[med][1:][idxvecdr2[jj]],np.mean(LeftCounts[med][1:][0:20]), bkg))
Intensityleft.append(np.mean(LeftCounts[med][1:][0:20])-bkg)
#errorIntensityleft.append(np.sqrt(np.mean(LeftCounts[med][1:][0:20]))+np.sqrt(bkg))
errorIntensityleft.append(2*np.std(LeftCounts[med][1:][0:20]))
plt.plot([2*f*1e-6 for f in LeftFrequencies[med][1:]], [c for c in LeftCounts[med][1:]], '-o', markersize=2, alpha=0.7)
plt.plot([2*f*1e-6 for f in LeftFrequencies[med][1:]][idxtest], [c for c in LeftCounts[med][1:]][idxtest], 'o',markersize=14)
jj=jj+1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('Counts')
plt.xlim(432, 446.5)
plt.grid()
plt.legend()
#plt.title('Espectros para distintas geometrías')
plt.figure()
plt.errorbar(Intensityleft, pmleftdepthsdr1, yerr=errorpmleftdepthsdr1, xerr=errorIntensityleft, label='DR left', fmt='o',capsize=2, markersize=8)
#plt.errorbar(Intensity, pmleftdepthsdr2, yerr=errorpmleftdepthsdr2, label='DR right', fmt='o',capsize=2, markersize=8)
plt.xlabel('Intensity')
plt.ylabel('DR Relative depth')
#plt.xticks([1,2,3,4,5])
plt.xlim(200,3200)
plt.ylim(0,1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()
#%%
import seaborn as sns
"""
Resonancias DD configuracion +2/-2 colineal variando la ubicacion del ion en los haces
LOBULO RIGHT DEL OAM
"""
palette = sns.color_palette("tab10")
"""
La 8, 9 y 10 son de la misma medicion. Hay que mergearlas
Por ahora obvio esto pero despues hay que hacerlo.
"""
pmrightlocmedvec = [0,4,3,2,1,5,6,7,8]
#pmrightlocmedvec = [8,9,10]
idxvecdr1 = [120,125,123,123,123,124,125,125,125]
idxvecdr2 = [406,406,405,404,404,407,407,406,406]
#pmlocmedvec = [2]
#LocVecsint = np.arange(1,6,1)
plt.figure()
bkg = 150
pmrightdepthsdr1=[]
errorpmrightdepthsdr1=[]
pmrightdepthsdr2=[]
errorpmrightdepthsdr2=[]
Intensityright = []
errorIntensityright = []
idxtest = 125
jj=0
for med in pmrightlocmedvec:
print(med)
pmrightdepthsdr1.append(1-(RightCounts[med][1:][idxvecdr1[jj]]-bkg)/(np.mean(RightCounts[med][1:][0:20])-bkg))
pmrightdepthsdr2.append(1-(RightCounts[med][1:][idxvecdr2[jj]]-bkg)/(np.mean(RightCounts[med][1:][0:20])-bkg))
errorpmrightdepthsdr1.append(ErrorDRdepth(RightCounts[med][1:][idxvecdr1[jj]],np.mean(RightCounts[med][1:][0:20]), bkg))
errorpmrightdepthsdr2.append(ErrorDRdepth(RightCounts[med][1:][idxvecdr2[jj]],np.mean(RightCounts[med][1:][0:20]), bkg))
Intensityright.append(np.mean(RightCounts[med][1:][0:20])-bkg)
#errorIntensityright.append(np.sqrt(np.mean(RightCounts[med][1:][0:20]))+np.sqrt(bkg))
errorIntensityright.append(2*np.std(RightCounts[med][1:][0:20]))
plt.plot([2*f*1e-6 for f in RightFrequencies[med][1:]], [c for c in RightCounts[med][1:]], '-o', markersize=2, alpha=0.7)
#plt.plot([2*f*1e-6 for f in RightFrequencies[med][1:]][idxtest], [c for c in RightCounts[med][1:]][idxtest], 'o',markersize=14)
jj=jj+1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('Counts')
#plt.xlim(432, 446.5)
#plt.xlim(443, 446)
plt.grid()
plt.legend()
#plt.title('Espectros para distintas geometrías')
plt.figure()
plt.errorbar(Intensityright, pmrightdepthsdr1, yerr=errorpmrightdepthsdr1, xerr=errorIntensityright, label='DR right', fmt='o',capsize=2, markersize=8)
#plt.errorbar(Intensityright, pmrightdepthsdr2, yerr=errorpmrightdepthsdr2, label='DR right', fmt='o',capsize=2, markersize=8)
plt.xlabel('Intensity')
plt.ylabel('DR Relative depth')
#plt.xticks([1,2,3,4,5])
plt.xlim(200,3200)
plt.ylim(0,1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()
#%%
"""
Mergeamos los dos lobulos
"""
#orderleft = [5,6,7,8,9,10,4,3,2,1,0]
orderleft = [4,5,6,7,8,3,2,1,0]
orderright = [4,3,2,1,0,5,6,7,8]
IntLeftFinal = [x for _, x in sorted(zip(orderleft, Intensityleft))]
IntRightFinal = [x for _, x in sorted(zip(orderright, Intensityright))]
PmLeftDR1Final = [x for _, x in sorted(zip(orderleft, pmleftdepthsdr1))]
PmRightDR1Final = [x for _, x in sorted(zip(orderright, pmrightdepthsdr1))]
PmLeftDR2Final = [x for _, x in sorted(zip(orderleft, pmleftdepthsdr2))]
PmRightDR2Final = [x for _, x in sorted(zip(orderright, pmrightdepthsdr2))]
errorIntLeftFinal = [x for _, x in sorted(zip(orderleft, errorIntensityleft))]
errorIntRightFinal = [x for _, x in sorted(zip(orderright, errorIntensityright))]
errorPmLeftDR1Final = [x for _, x in sorted(zip(orderleft, errorpmleftdepthsdr1))]
errorPmRightDR1Final = [x for _, x in sorted(zip(orderright, errorpmrightdepthsdr1))]
errorPmLeftDR2Final = [x for _, x in sorted(zip(orderleft, errorpmleftdepthsdr2))]
errorPmRightDR2Final = [x for _, x in sorted(zip(orderright, errorpmrightdepthsdr2))]
IntensityTotal = IntLeftFinal+IntRightFinal
errorIntensityTotal = errorIntLeftFinal+errorIntRightFinal
pmdepthsdr1 = PmLeftDR1Final+PmRightDR1Final
pmdepthsdr2 = PmLeftDR2Final+PmRightDR2Final
errorpmdepthsdr1 = errorPmLeftDR1Final+errorPmRightDR1Final
errorpmdepthsdr2 = errorPmLeftDR2Final+errorPmRightDR2Final
plt.figure()
plt.errorbar(np.arange(+0.5-len(pmdepthsdr1)/2,0.5+len(IntensityTotal)/2,1),[i/np.max(IntensityTotal) for i in IntensityTotal], yerr=[i/np.max(IntensityTotal) for i in errorIntensityTotal], label='Intensity', fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0.5-len(pmdepthsdr1)/2,0.5+len(pmdepthsdr1)/2,1),pmdepthsdr1, yerr=errorpmdepthsdr1, label='DR right', fmt='o',capsize=2, markersize=8)
#plt.errorbar(np.arange(0,len(pmdepthsdr2),1),pmdepthsdr2, yerr=errorpmdepthsdr2, label='DR right', fmt='o',capsize=2, markersize=8)
plt.grid()
#%%
from scipy.special import genlaguerre
#%%
import seaborn as sns
"""
Resonancias DD configuracion +2/-2 colineal variando la ubicacion del ion en los haces
Ploteo mediciones con alta estadistica para mirar (las extra) y ver cuán finitas son
"""
def Gaussian(x,A,B,x0,gam):
return A*np.exp((-1)*(x-x0)**2/(2*gam*gam)) + B
def Lorentzian( x, A, B, x0, gam ):
return A * gam**2 / ( gam**2 + ( x - x0 )**2) + B
palette = sns.color_palette("tab10")
pmextralocmedvec = [2]
plt.figure()
ftrap = 22.1
DR1 = 435.8
DR2 = 444.2
bkg = 150
jj=0
for med in pmextralocmedvec:
Freqs = [2*f*1e-6 for f in ExtraFrequencies[med][1:]]
Counts = [c for c in ExtraCounts[med][1:]]
#popt, pcov = curve_fit(Gaussian, Freqs, Counts, p0=(-200,2100,444.2,0.05))
popt, pcov = curve_fit(Lorentzian, Freqs, Counts, p0=(-2000,2100,444.2,0.05), bounds=((-10000,0,0,0),(0,1e4, 1e4, 1)))
plt.plot(Freqs, Counts, '-o', color=palette[med],markersize=2, alpha=0.5)
#plt.plot(Freqs, Gaussian(Freqs, *popt), '-o', color=palette[med],markersize=2, alpha=1, label=f'{LocVecs[jj]}')
plt.plot(Freqs, Lorentzian(Freqs, *popt), '-o', color=palette[med],markersize=2, alpha=1)
jj=jj+1
plt.xlabel('Frecuencia (MHz)')
plt.ylabel('Counts')
plt.xlim(443.9, 444.5)
plt.grid()
#plt.legend()
plt.title(f'Ancho de la resonancia: {round(1e3*popt[3],1)} kHz')
import h5py
import matplotlib.pyplot as plt
import numpy as np
import sys
import re
import ast
from scipy.optimize import curve_fit
import os
from scipy import interpolate
#Mediciones barriendo angulo del TISA y viendo kicking de resonancias oscuras
#C:\Users\Usuario\Documents\artiq\artiq_experiments\analisis\plots\20220106_CPT_DosLaseres_v08_TISA_DR\Data
os.chdir('/home/nico/Documents/artiq_experiments/analisis/plots/20230817_RotationalDopplerShift_v5/Data')
"""
en este codigo ploteo espectros CPT de resonancias D-D para configuracion +2/+2 y +2/-2 (usando pentaprisma)
"""
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
def Split(array,n):
length=len(array)/n
splitlist = []
jj = 0
while jj<length:
partial = []
ii = 0
while ii < n:
partial.append(array[jj*n+ii])
ii = ii + 1
splitlist.append(partial)
jj = jj + 1
return splitlist
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()))
Piezo1Counts = []
Piezo1Frequencies = []
PIEZO1_FILES = np.arange(791, 826,1)
for i in PIEZO1_FILES:
#print(str(i) + ' - ' + fname)
data = h5py.File(f'VaryingBeamlocation/Piezo/000014{i}-IR_Scan_withcal_optimized'+'.h5', 'r')
Piezo1Counts.append(np.array(data['datasets']['counts_spectrum']))
Piezo1Frequencies.append(np.array(data['datasets']['IR1_Frequencies']))
Piezo2Counts = []
Piezo2Frequencies = []
PIEZO2_FILES = list(np.arange(834, 841,1))+list(np.arange(842, 872,1))
for i in PIEZO2_FILES:
#print(str(i) + ' - ' + fname)
data = h5py.File(f'VaryingBeamlocation/Piezo2/000014{i}-IR_Scan_withcal_optimized'+'.h5', 'r')
Piezo2Counts.append(np.array(data['datasets']['counts_spectrum']))
Piezo2Frequencies.append(np.array(data['datasets']['IR1_Frequencies']))
def ErrorDRdepth(p, f, b):
ep = np.sqrt(p)
ef = np.sqrt(f)
eb = np.sqrt(b)
derivadap = 1/((f-b)**2)
derivadaf = ((p-b)/((f-b)**2))**2
derivadab = ((p-f)/((f-b)**2))**2
return 2*np.sqrt(derivadap*ep*ep + derivadaf*ef*ef + derivadab*eb*eb)
#%%
import seaborn as sns
"""
Resonancias DD configuracion +2/-2 colineal variando la ubicacion del ion en los haces
TODO EL OAM
"""
palette = sns.color_palette("tab10")
pmlocmedvec = list(np.arange(0,12,1))+[13,12]+list(np.arange(15,len(PIEZO1_FILES),1))
"""
Hay que invertir la 12 con la 13, y la 14 es la misma que la 12, por las dudas
"""
#pmlocmedvec = [21]
idxvecdr = [85,185, 185,182,182,232, 162,162,175,177,217, 217,184,184,188,186,186,195,195,198,195,198,196,192,192,192,178,175,175,170,140,140,190,185]
plt.figure()
bkg = np.min(Piezo1Counts[1])
pmdepthsdr=[]
errorpmdepthsdr=[]
Intensity = []
errorIntensity = []
idxtest = 185
print(idxtest)
jj=0
for med in pmlocmedvec:
print(med)
if med == 21 or med == 22 or med == 23:
pmdepthsdr.append(1-(Piezo1Counts[med][1:][idxvecdr[jj]]-bkg)/(np.mean(Piezo1Counts[med][1:][0:20])-bkg))
errorpmdepthsdr.append(ErrorDRdepth(Piezo1Counts[med][1:][idxvecdr[jj]],np.mean(Piezo1Counts[med][1:][0:20]), bkg))
Intens = np.mean(Piezo1Counts[med][1:][0:10])-bkg
Intensity.append(Intens)
errorIntensity.append(2*np.sqrt(np.mean(Piezo1Counts[med][1:][0:10]))+np.sqrt(bkg))
else:
pmdepthsdr.append(1-(Piezo1Counts[med][1:][idxvecdr[jj]]-bkg)/(np.mean(Piezo1Counts[med][1:][0:20])-bkg))
errorpmdepthsdr.append(ErrorDRdepth(Piezo1Counts[med][1:][idxvecdr[jj]],np.mean(Piezo1Counts[med][1:][0:20]), bkg))
Intens = np.mean(Piezo1Counts[med][1:][0:20])-bkg
Intensity.append(Intens)
errorIntensity.append(2*np.sqrt(np.mean(Piezo1Counts[med][1:][0:20]))+np.sqrt(bkg))
plt.plot([2*f*1e-6 for f in Piezo1Frequencies[med][1:]], [c for c in Piezo1Counts[med][1:]], '-o', markersize=2, alpha=0.7)
plt.plot([2*f*1e-6 for f in Piezo1Frequencies[med][1:]][idxtest], [c for c in Piezo1Counts[med][1:]][idxtest], 'o',markersize=14)
jj=jj+1
# plt.xlabel('Frecuencia (MHz)')
# plt.ylabel('Counts')
# plt.xlim(432, 446.5)
# plt.grid()
# plt.legend()
# #plt.title('Espectros para distintas geometrías')
plt.figure()
plt.errorbar(np.arange(0,len(Intensity),1), [i/np.max(Intensity) for i in Intensity], yerr=[i/np.max(Intensity) for i in errorIntensity], fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity),1), [p for p in pmdepthsdr], yerr=errorpmdepthsdr, fmt='o',capsize=2, markersize=8)
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.xticks([1,2,3,4,5])
#plt.xlim(200,3200)
plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()
#%%
import seaborn as sns
"""
Resonancias DD configuracion +2/+2 colineal variando la ubicacion del ion en los haces
TODO EL OAM
"""
palette = sns.color_palette("tab10")
mmlocmedvec = list(np.arange(0,10,1))+[11,10]+list(np.arange(12,len(PIEZO2_FILES),1))
"""
s
"""
#mmlocmedvec = [18]
#idxvecdr = [126,129,129,129,128,128,216,215,215,215,242,215,162,162,162,162,162,132,138,129,242,182,182,208,205,205,205,181,140,170,140,140,140,138,128,128,126,]
idxvecdr = [126,129,129,129,128,128,216,215,215,215,242,215,162,162,162,162,162,132,132,129,242,182,182,208,205,205,205,181,140,170,140,140,140,138,128,128,126,]
#idxtest = idxvecdr[18]
plt.figure()
bkg = np.min(Piezo2Counts[1])
mmdepthsdr=[]
errormmdepthsdr=[]
Intensity2 = []
errorIntensity2 = []
print(idxtest)
jj=0
for med in mmlocmedvec:
print(med)
mmdepthsdr.append(1-(Piezo2Counts[med][1:][idxvecdr[jj]]-bkg)/(np.mean(Piezo2Counts[med][1:][0:20])-bkg))
errormmdepthsdr.append(ErrorDRdepth(Piezo2Counts[med][1:][idxvecdr[jj]],np.mean(Piezo2Counts[med][1:][0:20]), bkg))
Intens = np.mean(Piezo2Counts[med][1:][0:20])-bkg
Intensity2.append(Intens)
errorIntensity2.append(2*np.sqrt(np.mean(Piezo2Counts[med][1:][0:20]))+np.sqrt(bkg))
plt.plot([2*f*1e-6 for f in Piezo2Frequencies[med][1:]], [c for c in Piezo2Counts[med][1:]], '-o', markersize=2, alpha=0.7)
plt.plot([2*f*1e-6 for f in Piezo2Frequencies[med][1:]][idxtest], [c for c in Piezo2Counts[med][1:]][idxtest], 'o',markersize=14)
jj=jj+1
# plt.xlabel('Frecuencia (MHz)')
# plt.ylabel('Counts')
# plt.xlim(432, 446.5)
# plt.grid()
# plt.legend()
# #plt.title('Espectros para distintas geometrías')
plt.figure()
#plt.errorbar(np.arange(0,len(Intensity2),1), [i/np.max(Intensity2) for i in Intensity2], yerr=[i/np.max(Intensity2) for i in errorIntensity2], fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity2),1), [p for p in mmdepthsdr], yerr=errormmdepthsdr, fmt='o',capsize=2, markersize=8)
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.xticks([1,2,3,4,5])
#plt.xlim(200,3200)
plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()
#%%
plt.figure()
plt.errorbar(np.arange(0,len(Intensity),1), [i/np.max(Intensity) for i in Intensity], yerr=[i/np.max(Intensity) for i in errorIntensity], fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity),1), [p for p in pmdepthsdr], yerr=errorpmdepthsdr, fmt='o',capsize=2, markersize=8)
plt.errorbar(np.arange(0,len(Intensity),1), [p for p in mmdepthsdr[3:]], yerr=[0.5*m for m in errormmdepthsdr[3:]], fmt='o',capsize=2, markersize=8)
plt.xlabel('Ion position')
plt.ylabel('Intensity / DR Relative depth')
#plt.xticks([1,2,3,4,5])
#plt.xlim(200,3200)
plt.ylim(-0.1,1.1)
plt.grid()
#plt.axvline(3, color='salmon')
plt.legend()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment