Dynamics.py

import math
import numpy as np
import scipy.constants as const

data_folder = './Results/'
pos_fname = 'position.npz'
vel_fname = 'velocity.npz'
time_fname = 'time.npz'

# Physical constants
c = const.c                     # light vel in vacuum (m/s)
kB = const.Boltzmann            # Boltzmann constant (J/K)
g = const.g                     # Grav. Acc. (m/s^2)
hbar = const.hbar               # reduced Planck constant (J*s)

#Trap                           # power beam (W)
R_trap = 30E-6                  # Fiber radius (m)

# Rb
m_Rb = 87 * const.physical_constants['atomic mass constant'][0] # Mass Rb 87
Gamma_1 = 2 * np.pi * 6.065e6   # Natural linewidth in Hz
Gamma_2 = 2 * np.pi * 5.746e6
omega_1 = 2 * np.pi * 377.1075 * 1e12 # Transition frequency 5S1/2 -> 5P3/2 in Hz
omega_2 = 2 * np.pi * 384.2305 * 1e12 # Transition frequency 5S1/2 -> 5S3/2  in Hz

def print_physical_constants():

    """
    Print all physical constants, trap parameters, Rb constants,
    and derived scales used in the MOT simulation.
    """

    print("\n=== PHYSICAL CONSTANTS ===")
    print(f"Speed of light (c): {c:.2e} m/s")
    print(f"Boltzmann constant (kB): {kB:.2e} J/K")
    print(f"Gravity (g): {g:.2e} m/s²")

    print("\n=== TRAP PARAMETERS ===")
    print(f"Fiber dimension (R_trap): {R_trap:.2e} m")
    print(f"Beam waist (w0): {w0:.2e} m (length scale along radial coordinate)")

    print("\n=== RUBIDIUM CONSTANTS ===")
    print(f"Rb-87 mass (m_Rb): {m_Rb:.2e} kg")
    print(f"Gamma_1: {Gamma_1:.2e} Hz")
    print(f"Gamma_2: {Gamma_2:.2e} Hz")
    print(f"Transition frequency omega_1: {omega_1:.2e} Hz")
    print(f"Transition frequency omega_2: {omega_2:.2e} Hz")

    print("\n==============================\n")

def alpha_func(lambda_b):
    omega = 2 * np.pi * c / lambda_b
    Delta_1 = omega - omega_1  # Detuning from the transition frequency
    Delta_2 = omega - omega_2  # Detuning from the transition frequency
    alpha = - np.pi * c**2 * (Gamma_1/(Delta_1*omega_1**3) + 2 * Gamma_2/(Delta_2*omega_2**3))
    return alpha

w0 = 19E-6                      # length scale along rho ( Beam waist at fiber tip )
def zR_func(lambda_b):
    # scales
    zR = np.pi * w0**2/lambda_b   # length scale along z ( Rayleigh length )
    return zR

def w(z: float, zR: float):
    """
    Gaussian beam radius.

    Parameters
    ----------
    z : float
        Axial coordinate.

    Returns
    -------
    beta : float
        Gaussian beam radius at z:
            w(z) = w0 * np.sqrt(1 + (z/zR)**2).
    """
    return w0 * np.sqrt(1 + (z/zR)**2)

if __name__ == '__main__':
    print_physical_constants()