Here is my code:
# Libraries
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import quad
# Constant parameters
m_Rb = 1.443*10**-25 #mass of rubidium 87
k_b = 1.38*10**-23
hbar = 1.05*10**-34
L = 38.116*10**6 #spontaneous decay rate
epsilon_0 = 8.85418782*10**-12 #permittivity of free space
# Changable paramaters
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 #Bohr magneton value
B = 5*10**-4 #magnetic field strength
# Scattering force
I_sat = 1 #intensity
n0 = 1 #refraction constant for medium
E_0 = ((2*I_sat)/(3*10**8*n0*epsilon_0))
Rabi = (E_0/hbar) * (4.2275*1.6*10**-19*0.529*10**-10)
n = 20
delta_omega = np.array([-20*10**6, -15*10**6, -10*10**6, -5*10**6])
def F(v):
i = 0
R_i = []
while i<len(delta_omega):
delta = delta_omega[i] + (k*v)
R_scat = L/2 * (Rabi**2/2)/(delta**2+(Rabi**2/2)+(L**2/4)) # scattering rate
R_i.append(R_scat)
i = i+1
R = np.sum(R_i)
R_total = (L*R)/(2*R + L)
def dXdt(t, X):
G, E = X
dGdt = E*(L+R_total) - G*R_total
dEdt = G*R_total - E*(L+R_total)
return [dGdt, dEdt]
t = np.linspace(0, 5, n)
solve = odeint(dXdt, [1, 0], t, tfirst=True)
G = solve.T[0]
E = solve.T[1]
F = hbar*k*(G-E)*R_total
return F
#stopping distance
ds = []
dt = 0.00001
x = np.zeros(n)
time = np.zeros(n)
#v = np.linspace(-700, 700, n)
v = np.zeros(n)
time[0] = 0
x[0] = 0
v[0] = 800
a = 0
print(n)
for a in np.arange(1, n):
time[a]=dt*a
solve = F(v=a)
F_int=solve[1]
v[a] = v[a-1] - dt*(F_int/m_Rb)*x[a-1]
x[a] = x[a-1] + dt*v[a-1]
for s in x:
if abs(s)<=length and abs(s)>0:
ds.append(abs(s))
# Capture velocity
dss = np.max(ds)
I want to find the value of velocity for which ds is at its maximum value (i.e. what's the maximum velocity that a particle can have and not exceed 0.05) but I don't know how to do this. I have calculated the distances travelled by a number of particles, and then made an array (ds) containing the particles which stop within the value of length. dss gives the maximum distance a particle has travelled, but it is the corresponding velocity for this value that I want to find.