I want the while loop at the end of my code to perform the calculation for each component of delta_omega, and print the result before moving on to the next number in the delta_omega array. When I run the code python outputs 'The rate of capture at a detuning of -5000000 is 3.0E+13' an infinite number of times before I manually stop it. I am unsure of why this is happening and how to fix it? I tried using break at the end of the loop, but this just performed the calculation once for the first component of the delta_omega array.
import numpy as np
from scipy.integrate import odeint
#import matplotlib.pyplot as plt
# Constants
m_Rb = 1.443*10**-25 #mass of rubidium 87
k_b = 1.38*10**-23
h = 6.63*10**-34
hbar = 1.05*10**-34
L = 38.116*10**6 #natural linewidth
epsilon_0 = 8.85418782*10**-12 #permittivity of free space
# Changable paramaters
lmbda = 780*10**-9 #wavelength of laser light
k = (2*np.pi)/lmbda #wavevector of laser light
B = 5*10**-4 #magnetic field strength
# D2 effective magetic moment
gj_gnd = 1 + (0.5*(0.5+1) + 0.5*(0.5+1) - 0*(0+1))/(2*0.5*(0.5+1))
mj_gnd = 0.5
gj_ex = 1 + (1.5*(1.5+1) + 0.5*(0.5+1) - 1*(1+1))/(2*1.5*(1.5+1))
mj_ex = 1.5
Bohr = 9.274*10**-24 #Bohr magneton value
mu_eff = Bohr*(gj_ex*mj_ex - gj_gnd*mj_gnd)
# -------- Before slower --------
T = 700 #temperature of oven
vp = ((2*k_b*T)/m_Rb)**0.5 #mean velocity of particles coming out of the oven
x_os = 0.1
a = (hbar*L*k)/(2*m_Rb) #max decelleration of atoms
vf_oven = (vp**2 + (2*a*x_os))**0.5
t_b = (2*x_os)/(vp + vf_oven) #time taken from oven to start of slower
# -------- During slower --------
length_slow = 0.5
vz_max = (vf_oven**2 + 2*a*length_slow)**0.5
Z = 0.7
#Z = np.linspace(0, length_slow, 100)
P = 10**(4.312-(4040/T)) #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) #atomic number density
I = 1*10**5 #intensity
n0 = 1 #refraction constant for medium
E_0 = ((2*I)/(3*10**8*n0*epsilon_0))**0.5
Rabi = (E_0*3.5844*10**-29)/hbar
II_sat = (2*Rabi**2)/L**2
delta_omega = np.array([-5*10**6, -10*10**6, -30*10**6]) #range of frequencies
i = 0
while i<len(delta_omega):
B_p = (h/mu_eff) * (delta_omega[i] + (1/lmbda)*(vf_oven**2 - (2*a*length_slow))**0.5)
B_n = (h/mu_eff) * (delta_omega[i] - (1/lmbda)*(vf_oven**2 - (2*a*length_slow))**0.5)
delta_n = delta_omega[i] + (k*vf_oven) - (mu_eff*B_n)/hbar
delta_p = delta_omega[i] - (k*vf_oven) + (mu_eff*B_p)/hbar
F = (hbar*k*L)/2 * ((II_sat/(1+II_sat+(2*delta_n/L)**2))
- (II_sat/(1+II_sat+(2*delta_p/L)**2)))
accn = abs(F/m_Rb)
vf_slower = (vf_oven**2 - (2*accn*length_slow))**0.5
t_d = 1/accn * (vf_oven - vf_slower) #time taken during slower
# -------- After slower --------
da = 0.1 #distance from end of slower to the middle of the MOT
vf_MOT = (vf_slower**2 - (2*accn*da))**0.5
t_a = da/vf_MOT #time taken after slower
r0 = 0.01 #MOT capture radius
vr_max = r0/(t_b+t_d+t_a)
# -------- Flux of atoms captured --------
f_oven = ((n*A)/4) * (2/(np.pi)**0.5) * ((2*k_b*T)/m_Rb)**0.5
f = f_oven * (1 - np.exp(-vr_max**2/vp**2))*(1 - np.exp(-vz_max**2/vp**2))
print('The rate of capture at a detuning of', delta_omega[i], 'is', format(f, '.1E'))