Calculating the differential at multiple frequencies

Viewed 36
# Libraries
import numpy as np
from scipy.integrate import quad

# Tow counter propagating beams in x-axis
laser_x = np.linspace(0, 100, 100)
laser_negx= np.linspace(0, -100, 100)

# Constant parameters
m_Rb = 1.443*10**-25 #mass of rubidium 87
k_b = 1.38*10**-23
hbar = 1.05*10**-34
Rabi = 46.567*10**6 #Rabi frequency
L = 38.116*10**6 #spontaneous decay rate

# Changable paramaters
#delta_omega = -20*10**-6
delta_omega = np.array([-20*10**6, -15*10**6, -10*10**6, -5*10**6]) #difference in the laser frequency and the atomic resonance frequency
lmbda = 700*10**-9 #wavelength of laser light
k = (2*np.pi)/lmbda #wavevector of laser light
V = 1.25*10**-4 #volume of MOT space
length = 5*10**-2 #length of MOT
Bohr = 9.274*10**-24
B = 1

# Maxwell Boltzmann distribution variables
T = 300


# Number of particles emittied
T_oven = 300 #oven temperature
P = 10**(4.312-(4040/T_oven)) #vapour pressure for liquid phase (use 4.857 for solid phase)
A = 5*10**-4 #area of the oven aperture

n = P/(k_b*T_oven) #atomic number density
I_oven = ((n*A)/4) * (2/(np.pi)**0.5) * ((2*k_b*T_oven)/m_Rb)**0.5
print("The flux of atoms from the oven is", format(I_oven, '.1E'))

# Calculating the scattering force and the distance travelled by scattered particles
ds = []

n = 800
dt = 0.001
t = np.zeros(n)
v = np.linspace(-600, 600, n)
x = np.zeros(n)
t[0] = 0
x[0] = 0
w=0
# solving the differential
def f(v):
    return 4*np.pi* (m_Rb/(2*np.pi*k_b*T))**(1.5) * v**3 * np.exp((-m_Rb*v**2)/(2*k_b*T))

while w < len(delta_omega):
    for i in np.arange(1, n):
        delta = delta_omega[w] - (k*v[i]) + ((Bohr/hbar)*B) #result of detuning the laser
        F = (hbar * k * (L/2) * ((Rabi**2/2)/(delta**2+(Rabi**2/2)+(L**2/4))))/m_Rb
        t[i]=dt*i
        v[i] = v[i-1] - dt*F*x[i-1]
        x[i] = x[i-1] + dt*v[i-1] 
        for s in x:
            if s<length:
                ds.append(s)
    v_capture = abs(np.min(ds)) 
    Rate, error = quad(f, 0, v_capture)
    print("The rate of capture at", format(delta_omega[w], '.1E'), 'Hz is', format(Rate, '.1E'))
    w=w+1

I want to calculate the value of Rate at each of the components of delta_omega and then print that result, before repeating the calculation with the next value in the delta_omega array. Here is the result I am getting:

runcell(0, 'C:/Users/chris/Desktop/Internship/MOT Capture Multiple Frequencies.py')
The flux of atoms from the oven is 5.7E+09
The rate of capture at -2.0E+07 Hz is 2.5E+02
The rate of capture at -1.5E+07 Hz is 2.5E+02
The rate of capture at -1.0E+07 Hz is 2.5E+02
The rate of capture at -5.0E+06 Hz is 2.5E+02

I can't tell if this is just how the calculation is coming out or if something is wrong with my loop?

0 Answers
Related