Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from OBE3levels_functions import *
from scipy.optimize import curve_fit
from scipy.stats import poisson, norm
import os
import sys
plt.ion()
"""
Este código simula la dinámica temporal de un ion de calcio sometido a un láser UV y dos láseres IR.
"""
def expo(T, tau, N0, C):
#global T0
return N0*np.exp(-(T-0e-6)/tau) + C
ADD_NOISE = True
NOISE_VAR = 30 # Varianza del ruido
NOISE_ADJ = 78 # Escaleo del ruido versus el maximo
NOISE_MEAN = 0 # Media del ruido
NRO_REPS = 100 # Veces a promediar
species = 'Calcium'
#Parámetros de saturacion de los láseres rebombeo (IR), que son proporcionales a la intensidad de los mismos
#Para reproducir SP, hay que poner 0 a los dos
SatParRep1, SatParRep2 = 0, 0
sat_intensity = 43.3e-5 # uW um-2
g21 = 135591138.92893547
#El número de los láseres IR aplicados en el experimento. Para el experimento de la transicion SP, es 0.
Number_IR_lasers = 0
#Detunings de los láseres en MHz. Si es 0, están en resonancia. En principio
#podría ser un parámetro libre a ajustar porque influye mucho
DetDoppler = -int(sys.argv[1])
DetRep2 = 0
DetRep1 = 0
#Anchos de línea de los láseres, en MHz. No son tan cruciales, están entre 0.1 y 0.5 MHz
DopplerLaserLinewidth, RepumpLaserLinewidth = 1e3, 1
# Para trabajar con las curvas teóricas:
# range_inicio = np.arange(sat_intensity/10, sat_intensity, sat_intensity/50)
# range_fin = np.arange(sat_intensity, 30*sat_intensity, sat_intensity/5)
# laser_intensities = np.append(range_inicio, range_fin)
# Para trabajar con nuestros datos, saco las intensidades, a un radio de 100um:
laser_intensities = np.array([0.00662085, 0.00582507, 0.0049338, 0.00397887, 0.00305577, 0.00226, 0.00155972, 0.00101859, 0.00060479, 0.00031831])
OmegaRabi2 = (g21**2) * (laser_intensities / sat_intensity)
SatParDopplerVec = OmegaRabi2 / ((DetDoppler*2*np.pi*1e6)**2 + g21**2)
#Parámetros de la simulación, en us
tfinOBE = 1
pasoOBE = 1e-5
#Estado electrónico inicial de la simulacion. Posibilidades: |S> = 1, |P> = 2, |D> = 3
initial_state = 1
#Corre y simula la dinámica. Devuelve elementos de matriz densidad
# fig1, ax1 = plt.subplots()
alltaus = list()
allN0 = list()
allmeans = list()
allstds = list()
temporal_taus = list()
temporal_N0 = list()
for i, SatParDoppler in enumerate(SatParDopplerVec):
t, RealRho11, RealRho22, RealRho33, AbsRho12, AbsRho13, AbsRho23 = Solve3LevelBloch(EquationSystemModel, species, SatParDoppler, SatParRep1, SatParRep2, DetDoppler, DetRep1, DetRep2, DopplerLaserLinewidth, RepumpLaserLinewidth, no_IR_lasers=Number_IR_lasers, initcond=initial_state, tfin=tfinOBE, paso=pasoOBE, printprogress=False)
k = int(0.1*len(t)) #parametro para que fitee solo la parte final
k = np.argmin(np.abs(t - 0.2e-6))
temporal_taus = list()
incert_estad = 0
for nro_rep in range(NRO_REPS):
poblacion_fit = RealRho22[k:].copy()
if ADD_NOISE:
# Genero ruido similar al medido
ruido_poisson = poisson.rvs(NOISE_VAR, size=len(poblacion_fit))
# Reescaleo el ruido para que sea compatible con las simulaciones
ruido_poisson *= np.max(poblacion_fit)/NOISE_ADJ
poblacion_fit += ruido_poisson
# poblacion_fit += norm.rvs(0, NOISE_VAR, size=len(poblacion_fit))*np.max(poblacion_fit)/NOISE_ADJ
poblacion_fit[np.where(poblacion_fit < 0 )] = 0 # Limpio posibles valores negativos
t_fit = t[k+1:]
popt, pcov = curve_fit(expo, t_fit, poblacion_fit, p0=(1e-6, 1e-2, 1e-2))
temporal_taus.append(popt[0])
incert_estad += np.sqrt(pcov[0,0])
# temporal_N0.append(popt[1])
if not ADD_NOISE: # si no quiere que agregue ruido, entonces solo hace el for una vez
break
print(f" {nro_rep} ", end='\r')
if ADD_NOISE:
allmeans.append(np.mean(temporal_taus))
allstds.append(np.std(temporal_taus) + incert_estad/NRO_REPS)
else:
allmeans.append(np.mean(temporal_taus))
#ploteo en funcion de 100*Rho22 para que me devuelva en porcentaje la población del excitado que es proporcional a la fluorescencia detectada
# lbl='Detuning Rep1:' + str(DetRep1) + ' MHz'
# ax1.plot(t[:-1]*1e6, RealRho22, '.', label=lbl, markersize=1)
# ax1.plot(t[:-1]*1e6, 100*RealRho22, '.', label=lbl, markersize=1)
# ax1.plot(t_fit*1e6, 100*expo(t_fit, *popt))
print(f'({i:02d}/{len(SatParDopplerVec)}) Done {SatParDoppler}'.ljust(50, ' '))
# ax1.legend(SatParDopplerVec)
# ax1.set_xlabel('Time (us)')
# ax1.set_ylabel('100*Rho22')
#%% #####################################################################################
### Guardo valores al CSV para analizar luego
# nDet=f"Det{np.abs(DetDoppler)}"
# CSV_FNAME = f"error_teorico_tau_{nDet}.csv"
# try:
# import pandas as pd
# data = pd.read_csv(CSV_FNAME)
# data["meanval"] = allmeans
# data["stdval"] = allstds
# data.to_csv(CSV_FNAME, index=False)
# print(f"SAVED {nDet}".ljust(50, ' '))
# except FileNotFoundError:
# data = pd.DataFrame({'I': laser_intensities, "meanval": allmeans, "stdval": allstds})
# data.to_csv(CSV_FNAME, index=False)
# print(f"CREATED CSV FILE".ljust(50, ' '))
# print(f"SAVED DETUNING {-DetDoppler}".ljust(50, ' '))
#%% ####################################################################################
### Plots de los valores de intensidades medidas para distintos parametros
### Cada axes arma un eje horizontal distinto, es para comparar las posibilidades
# fig4, ax4 = plt.subplots()
# allmeans = [t*1e6 for t in allmeans]
# ax4b = ax4.twinx()
# ax4.plot(SatParDopplerVec, allmeans, '-o', lw=0.4)
# ax4b.plot(SatParDopplerVec, allN0, 'k-^', lw=0.4)
# ax4.set_xlabel("Parametro de saturación")
# ax4.set_ylabel("Tau (circulo)")
# ax4b.set_ylabel("Alturas (triang)")
# ax4.set_title(f"SP; IR_Lasers=0; DetuningDoppler={DetDoppler}; LineWidth={DopplerLaserLinewidth}")
# ax4by = ax4.twiny()
# ax4by.plot(laser_intensities, allmeans, '-', lw=0)
# ax4by.xaxis.set_label_position('bottom')
# ax4by.xaxis.tick_bottom()
# ax4by.spines['bottom'].set_position(("axes", -0.27))
# ax4by.set_xlabel(r"Intensidad [ $\mu$W/$\mu$m$^2$ ]")
# ax4cy = ax4.twiny()
# ax4cy.plot(laser_intensities*np.pi*(35**2), allmeans, '-', lw=0)
# ax4cy.xaxis.set_label_position('bottom')
# ax4cy.xaxis.tick_bottom()
# ax4cy.spines['bottom'].set_position(("axes", -0.58))
# ax4cy.set_xlabel(r"Potencia (r=35$\mu$m) [ $\mu$W ]")
# ax4dy = ax4.twiny()
# ax4dy.plot(laser_intensities*np.pi*(50**2), allmeans, '-', lw=0)
# ax4dy.xaxis.set_label_position('bottom')
# ax4dy.xaxis.tick_bottom()
# ax4dy.spines['bottom'].set_position(("axes", -0.9))
# ax4dy.set_xlabel(r"Potencia (r=50$\mu$m) [ $\mu$W ]")