Integration and utilities reference

Integrate module

Most functions in the mwave.integrate module use numba to increase performance.

mwave.integrate.bloch(kvec, phi0, tfinal, delta, omega, omega_args, phase, phase_args, t0=0, method='DOP853', atol=1e-10, rtol=1e-10, dense=False, max_step=0.1, transformed=False, Gamma_sps=None)

Evolves the provided wavefunction under the Bloch Hamiltonian. This function internally uses mwave.integrate.bloch_rhs() if Gamma_sps=None and mwave.integrate.bloch_density_rhs() if Gamma_sps is not None.

Parameters:

kvec – The vector of momentum states to simulate.
phi0 – The initial value of phi.
tfinal – The final time to integrate to.
delta – The two-photon detuning.
omega – The function that returns the value of the effective Rabi frequency $\Omega(t, a)$ at an arbitrary time. The function must take two arguments, t and omega_args. The argument t specifies the time at which to evaluate the effective Rabi frequency and the argument omega_args can be used to pass in additional parameters.
omega_args – A tuple of arguments to pass to the function defined by omega.
phase – The function that returns the phase of two photon detuning at an arbitrary time. The function must take two arguments, t and phase_args. The argument t specifies the time at which to evaluate the phase and the argument phase_args can be used to pass in additional parameters. This function can be set to a constant value if the user does not want to simulate a frequency swept process.
phase_args – A tuple of arguments to pass to the function defined by phase.
t0 – The initial time to start the simulation. Defaults to 0. This is useful if you wish to chain the output of one call to bragg with another call to bragg that occurs at a later time.
method – The integration method to call scipy.integrate.solve_ivp with. Defaults to 'DOP853'.
atol – The absolute tolerance given to scipy.integrate.solve_ivp.
atol – The relative tolerance given to scipy.integrate.solve_ivp.
dense – If true dense output is returned (i.e. the integration result can be queried for any intermediate time).
max_step – The max step size to use during the integration.
transformed – See the description of mwave.integrate.bloch_rhs() for details.
Gamma_sps – The rate of single photon scattering. Single photon scattering is only applied if this parameter is provided. When provided the mwave.integrate.bragg() function converts the provided wavefunction into a density matrix and then evolves the density matrix under the Bragg Hamiltonian with decoherence included. As such the returned solution contains the density matrix instead of a wavefunction.

Returns:

The solution object output from scipy.integrate.solve_ivp().

mwave.integrate.bloch_density_rhs(t, rho, nstates, hkvec, vkvec, loss_mat, delta, omega, omega_args, phase, phase_args)

Evaluates the right hand side of the Von Neumann evolution equation for the Bloch Hamiltonian (i.e. $[H,\rho]$) where

\[H=-\hbar\sum_{k}\left[\frac{\Omega_\text{eff}(t,a)}{2}e^{i(\delta t+\theta(t,b))}e^{i\omega_\text{r}(-4k-4)}|k\rangle\langle k+2|+\frac{\Omega_\text{eff}(t,a)^*}{2}e^{-i(\delta t+\theta(t,b))}e^{i\omega_\text{r}(4k-4)}|k\rangle\langle k-2|\right]\]

where $\hbar=1$ and the sum over $k$ is limited to the values of $k$ defined by hkvec and vkvec.

The parameter rho is supplied as a vector (this makes it compatible with scipy.integrate.solve_ivp). This is then converted to a matrix via np.reshape(rho, (len(kvec), len(kvec))) internally. The matrices hkvec and vkvec are composed of horizontal or vertical vectors of the momentum state grid stacked togeather.

The parameter loss_mat is the loss matrix.

The remaining parameters (delta, omega, omega_args, phase, phase_args) are equivalent to those used in the mwave.integrate.bloch_rhs() function.

Parameters:

t – The time at which to evaluate the right hand side.
rho – The value of rho at which to evaluate the right hand side.
nstates – The number of states in rho, used to properly reshape the density matrix.
hkvec – The momentum state values at which rho is defined along the horizontal axis.
vkvec – The momentum state values at which rho is defined along the vertical axis.
loss_mat – The loss matrix to use.
delta – The value of $\delta$ (the two-photon detuning).
omega – The function that returns the value of the effective Rabi frequency $\Omega(t, a)$ at an arbitrary time. The function must take two arguments, t and omega_args. The argument t specifies the time at which to evaluate the effective Rabi frequency and the argument omega_args can be used to pass in additional parameters.
omega_args – A tuple of arguments to pass to the function defined by omega.
phase – The function that returns the phase of two photon detuning at an arbitrary time. The function must take two arguments, t and phase_args. The argument t specifies the time at which to evaluate the phase and the argument phase_args can be used to pass in additional parameters. This function can be set to a constant value if the user does not want to simulate a frequency swept process.
phase_args – A tuple of arguments to pass to the function defined by phase.

Returns:

A vector containing the evaluated right hand side values.

mwave.integrate.bloch_rhs(t, phi, kvec, delta, omega, omega_args, phase, phase_args, transformed=False)

Evaluates the right hand side of the Schrodinger equation for the Bloch Hamiltonian. The function returns a vector, one for each state included in the Hamiltonian.

The right hand side is defined in a general way so that a time-dependent field intensity and phase can be computed.

\[\text{returned vector}=i\frac{\Omega(t, a)}{2}\left[e^{i(\delta t+\theta(t,b))}e^{i(-4k-4)t}\lvert k\rangle\langle k + 2\rvert + e^{-i(\delta t+\theta(t,b))}e^{i(4k-4)t}\lvert k\rangle\langle k-2\rvert\right]\rvert\phi\rangle\]

where $k$ indexes momentum states spaced by two photon recoils. The time $t$ is evaluated at t, and the state $\lvert\phi\rangle$ is specified by phi.

The states $k$ included in the calculation are specified by kvec. $\delta$ is specified by delta. $\Omega(t, a)$ is specified by omega, which must a function which takes arguments t and omega_args. $\theta(t, b)$ is specified by phase, which must a function which takes arguments t and phase_args.

If transformed is True then the right hand side is evaluated in the following frame:

\[\text{returned vector}=-ik^2\lvert k\rangle\langle k\rvert\phi\rangle + i\frac{\Omega(t, a)}{2}\left[e^{i(\delta t+\theta(t,b))}\lvert k\rangle\langle k + 2\rvert + e^{-i(\delta t+\theta(t,b))}\lvert k\rangle\langle k-2\rvert\right]\rvert\phi\rangle\]

To solve the Bloch Hamiltonian in time the mwave.integrate.bloch_rhs() function can be integrated using scipy.integrate.solve_ivp().

Parameters:

t – The time at which to evaluate the right hand side.
phi – The value of phi at which to evaluate the right hand side.
kvec – The momentum state values at which phi is defined.
delta – The value of $\delta$ (the two-photon detuning).
omega – The function that returns the value of the effective Rabi frequency $\Omega(t, a)$ at an arbitrary time. The function must take two arguments, t and omega_args. The argument t specifies the time at which to evaluate the effective Rabi frequency and the argument omega_args can be used to pass in additional parameters.
omega_args – A tuple of arguments to pass to the function defined by omega.
phase – The function that returns the phase of two photon detuning at an arbitrary time. The function must take two arguments, t and phase_args. The argument t specifies the time at which to evaluate the phase and the argument phase_args can be used to pass in additional parameters. This function can be set to a constant value if the user does not want to simulate a frequency swept process.
phase_args – A tuple of arguments to pass to the function defined by phase.
transformed – See the function description above.

Returns:

A vector containing the evaluated right hand side values.

mwave.integrate.bloch_rhs_gaussian(t, phi, kvec, delta, omega, sigma, t0)

Evaluates the right hand side of the Schrodinger equation for the Bloch Hamiltonian in the case of a Gaussian pulse with constant phase. The function returns a vector, one for each state included in the Hamiltonian.

The right hand side is given by

\[\text{returned vector}=i\frac{\Omega(t)}{2}\left[e^{i\delta t}e^{i(-4k-4)t}\lvert k\rangle\langle k + 2\rvert + e^{-i\delta t}e^{i(4k-4)t}\lvert k\rangle\langle k-2\rvert\right]\rvert\phi\rangle\]

where $k$ indexes momentum states spaced by two photon recoils. The time $t$ is evaluated at t, and the state $\lvert\phi\rangle$ is specified by phi.

The states $k$ included in the calculation are specified by kvec. $\delta$ is specified by delta. $\Omega(t, a)$ is specified as follows

\[\Omega(t)=\Omega\exp\left(-\frac{(t-t_0)^2}{2\sigma^2}\right)\]

where $\Omega$ is given by omega, $\sigma$ is given by sigma, and $t_0$ is given by t0.

The mwave.integrate.bloch_rhs_gaussian() function can be integrated in time using scipy.integrate.solve_ivp().

Parameters:

t – The time at which to evaluate the right hand side.
phi – The value of phi at which to evaluate the right hand side.
kvec – The momentum state values at which phi is defined.
delta – The value of $\delta$ (the two-photon detuning).
omega – The peak effective Rabi frequency.
sigma – The Gaussian width of the Rabi frequency in time.
t0 – The center time of the Gaussian.

Returns:

A vector containing the evaluated right hand side values.

mwave.integrate.bloch_rhs_multifreq_gaussian(t, phi, kvec, delta, omega, sigma, t0, omega_mod)

Evaluates the right hand side of the Schrodinger equation for the multifrequency Bloch Hamiltonian in the case of a Gaussian pulse with constant phase. The function returns a vector, one for each state included in the Hamiltonian.

The right hand side is given by

\[\text{returned vector}=i\frac{\Omega(t)}{2}\left[e^{i\delta t}e^{i(-4k-4)t}\lvert k\rangle\langle k + 2\rvert + e^{-i\delta t}e^{i(4k-4)t}\lvert k\rangle\langle k-2\rvert\right]\rvert\phi\rangle\]

where $k$ indexes momentum states spaced by two photon recoils. The time $t$ is evaluated at t, and the state $\lvert\phi\rangle$ is specified by phi.

The states $k$ included in the calculation are specified by kvec. $\delta$ is specified by delta. $\Omega(t)$ is specified as follows

\[\Omega(t)=2\Omega\cos(\omega_\text{mod}t)\exp\left(-\frac{(t-t_0)^2}{2\sigma^2}\right)\]

where $\Omega$ is given by omega, $\sigma$ is given by sigma, $t_0$ is given by t0, and $\omega_\text{mod}$ is given by omega_mod.

The mwave.integrate.bloch_rhs_multifreq_gaussian() function can be integrated in time using scipy.integrate.solve_ivp().

Parameters:

t – The time at which to evaluate the right hand side.
phi – The value of phi at which to evaluate the right hand side.
kvec – The momentum state values at which phi is defined.
delta – The value of $\delta$ (the two-photon detuning).
omega – The peak effective Rabi frequency.
sigma – The Gaussian width of the Rabi frequency in time.
t0 – The center time of the Gaussian.
omega_mod – The modulation frequency.

Returns:

A vector containing the evaluated right hand side values.

mwave.integrate.gbragg(kvec, phi0, tfinal, delta, omega, sigma, omega_mod=None, method='DOP853', atol=1e-10, rtol=1e-10, dense=False, max_step=0.1)

Performs Bragg diffraction with a Gaussian profile and a constant phase. This function internally uses mwave.integrate.bloch_rhs_gaussian() if omega_mod=None and mwave.integrate.bloch_rhs_multifreq_gaussian() if omega_mod is not None. The Gaussian center is automatically placed at tfinal/2.

Parameters:

kvec – The vector of momentum states to simulate.
phi0 – The initial value of phi.
tfinal – The final time to integrate to.
delta – The two-photon detuning.
omega – The peak effective Rabi frequency.
sigma – The Gaussian width of the Rabi frequency in time.
omega_mod – The modulation frequency.
method – The integration method to call scipy.integrate.solve_ivp() with. Defaults to 'DOP853'.
atol – The absolute tolerance given to scipy.integrate.solve_ivp().
atol – The relative tolerance given to scipy.integrate.solve_ivp().
dense – If true dense output is returned (i.e. the integration result can be queried for any intermediate time).
max_step – The max step size to use during the integration.

Returns:

The solution object output from scipy.integrate.solve_ivp().

>>> from mwave.integrate import make_kvec, make_phi, gbragg, pops_vs_time
>>> n0, nf = 0, 5
>>> sigma = 0.188
>>> omega= 30
>>> kvec, n0_idx, nf_idx = make_kvec(n0,nf)
>>> sol = gbragg(kvec, make_phi(kvec, n0), 6*sigma, 4*(n0+nf), omega, sigma)
>>> pops_vs_time(kvec, sol.t, sol.y.T)
>>> sol.y[:,-1]
array([ 3.42308432e-15+5.79652561e-15j,  1.16013222e-14+1.77518885e-15j,
        1.02665571e-14-1.09796741e-14j,  3.12034740e-14-4.32155254e-14j,
        4.39906564e-14-1.12239410e-13j, -4.03937279e-13+4.94920498e-14j,
       -2.24946429e-10+6.32936084e-11j,  1.72491002e-07-9.66758213e-08j,
        3.59608734e-05-9.70371487e-05j,  1.82135333e-02+1.69593062e-02j,
       -2.89021471e-01-6.78069289e-01j, -4.69228334e-02+1.65328320e-01j,
       -3.73301726e-02-7.29203499e-02j,  2.17166733e-02+5.01757478e-02j,
        8.47644900e-03-2.22739506e-02j,  5.87939231e-01-2.64231147e-01j,
       -1.86014548e-02-2.07189313e-02j, -3.36140745e-05+1.26132182e-04j,
       -2.13202778e-07+1.69602467e-07j,  3.21925038e-10-1.32610954e-10j,
        2.75292990e-13-1.35756738e-13j,  3.65468413e-15-6.73688241e-15j,
        4.95541492e-15+3.82119313e-15j, -1.97357818e-15-3.62804364e-15j,
       -2.39926148e-15+2.33858146e-15j, -1.94677429e-15+1.34722165e-15j])

mwave.integrate.integrate_continuous_to_discrete(kvec, psif, n2hk)

DEPRECATED

Integrates psif over each spacing of $2\hbar k$ on the momentum grid specified by kvec. Returns the mean value of kvec over each $2\hbar k$ and then the integral of np.abs(psif)**2 over each $2\hbar k$. For example usage see Simulating a momentum distribution.

mwave.integrate.kbragg(n0, nf, tfinal, delta, omega, omega_args, phase, phase_args, npad=10, method='DOP853', atol=1e-10, rtol=1e-10, dense=False, max_step=0.1): DEPRECATED

mwave.integrate.make_continuous_kvec(n0, nf, dk, npad=10)

DEPRECATED

Creates a continuous momentum grid with spacing dk. Note that the usual states in the Bragg Hamiltonian are spaced by $2\hbar k$. Here it is assumed that $\hbar k=1$, so the Bragg couples states distance $2$ apart. For example usage see Simulating a momentum distribution.

mwave.integrate.make_kvec(n0, nf, npad=10)

Generates a vector of $k$-states. Note that neighboring $k$-states are spaced by 2 photon recoils.

Parameters:

n0 – The initial momentum state to include in the space.
nf – The final momentums tate to include in the space.
npad – The padding to include on each side of the initial and final momentum states.

Returns:

A tuple containing a vector of momentum states spaced by $2\hbar k$, the index of n0 in the vector, and the index of nf in the vector.

Example

>>> from mwave.integrate import make_kvec
>>> n0, nf = 0, 5
>>> make_kvec(n0,nf)
(array([-20., -18., -16., -14., -12., -10.,  -8.,  -6.,  -4.,  -2.,   0.,
         2.,   4.,   6.,   8.,  10.,  12.,  14.,  16.,  18.,  20.,  22.,
        24.,  26.,  28.,  30.]), 10, 15)

mwave.integrate.make_phi(kvec, n0)

Creates a vector the same length as kvec filled with all zeros aside from in the index where kvec==2*n0, which is set to 1.

Parameters:

kvec – The momentum state values at which phi is defined.
n0 – The state to initialize all amplitude in.

mwave.integrate.multi_omega_fnc(t, args)

Function defining a multifrequency Gaussian pulse profile in time, i.e.

\[\Omega(t)=2\Omega\cos(\omega_\text{mod}t)\exp\left(-\frac{(t-t_0)^2}{2\sigma^2}\right)\]

where $\Omega$, $\sigma$, $t_0$, and $\omega_\text{mod}$ are given by args[0], args[1], args[2], and args[3], respectively.

Parameters:

t – The time at which to evalute the Gaussian.
args – A tuple of four parameters defining $\Omega$, $\sigma$, $t_0$, and omega_text{mod}.

Returns:

The function value at the provided time.

mwave.integrate.omega_fnc_gaussian(t, args)

Function defining a Gaussian pulse profile in time, i.e.

\[\Omega(t)=\Omega\exp\left(-\frac{(t-t_0)^2}{2\sigma^2}\right)\]

where $\Omega$, $\sigma$, and $t_0$ are given by args[0], args[1], and args[2], respectively.

Parameters:

t – The time at which to evalute the Gaussian.
args – A tuple of three parameters defining $\Omega$, $\sigma$, and $t_0$.

Returns:

The function value at the provided time.

mwave.integrate.opt_states(optfnc, arg_guess, n0_idx, nf_idx, pi2_weight=1, pi_weight=0, nontarget_weight=0)

DEPRECATED

Optimizes a pulse for a particular set of final states. mwave.integrate.optfunc() should return an array of amplitudes. The amplitudes at n0_idx and nf_idx are specified to be the original and target states. Different combinations of the final arguments optimize for different ratios of these two states.

>>> from mwave.integrate import gbragg, make_kvec, make_phi, opt_states, pops_vs_time
>>> n0, nf = 0, 5
>>> kvec, n0_idx, nf_idx = make_kvec(n0, nf)
>>> opt = opt_states(lambda x: gbragg(kvec, make_phi(kvec, n0), 6*x[1], 4*(n0+nf), x[0], x[1]), [30, 0.188], n0_idx, nf_idx, pi2_weight=0, pi_weight=1)
>>> sol = gbragg(kvec, make_phi(kvec, n0), 6*opt.x[1], 4*(n0+nf), opt.x[0], opt.x[1])
>>> pops_vs_time(kvec, sol.t, sol.y.T)
>>> sol.y[:,-1]

mwave.integrate.phase_fnc_constant(t, args)

Function defining a constant phase as a function of time.

Parameters:

t – The time at which to evalute the phase. Since the phase is constant this parameter has no effect.
args – A tuple of one parameters defining the value of the constant phase.

Returns:

The phase value at the provided time.

mwave.integrate.pops_vs_time(kvec, t, phi, ax=None, legend=False)

Plots the population of each momentum state as a function of time.

Parameters:

kvec – The momentum state vector.
t – The time vector.
phi – The wavefunction array with shape (len(t), len(kvec)).
ax – The matplotlib axis to plot on. If None, a new figure is created.
legend – If True, a legend is added to the plot.

Simulation utils module

mwave.simulation_utils.cloud_init(natoms, sigma_cloud, sigma_transverse_v, sigma_vertical_v, x_offset=0, y_offset=0, z_offset=0, vx_offset=0, vy_offset=0, vz_offset=0, seed=None)

Basic cloud initialization. For more complicated atom clouds (i.e. clouds where positions/velocities are correlated with each other), the user should write their own initialization function.

Parameters:

natoms – The number of atoms to initialize.
sigma_cloud – The standard deviation of atom position within the cloud. The user can either provide a single number (in which case the standard deviation is the same in all directions), or three numbers (in which case the standard deviation is defined as sigma_cloud=(sigma_x, sigma_y, sigma_z)).
sigma_transverse_v – The standard deviation of atom velocity within the cloud. The user can either provide a single number (in which case the standard deviation is same in both transverse directions), or two numbers (in which case the standard deviation is defined as sigma_transverse_v=(sigma_x, sigma_y)).
sigma_vertical_v – The standard deviation of atom velocity with the vertical direction of the cloud.
x_offset – The offset of the cloud from $x=0$.
y_offset – The offset of the cloud from $y=0$.
z_offset – The offset of the cloud from $z=0$.
vx_offset – The offset of the cloud velocity from $v_x=0$.
vy_offset – The offset of the cloud velocity from $v_y=0$.
vz_offset – The offset of the cloud velocity from $v_z=0$.
seed – The seed to use when drawing the positions and velocities. If this function is called within an optimization loop the seed should remain the same across function calls!

Returns:

A tuple of (x0, y0, z0, vz, vx, vy), where each element of the tuple is a numpy array of length natoms.

mwave.simulation_utils.cloud_to_scrbi_ellipse_xy(x0, y0, z0, vx, vy, vz, T, Tp, n, N, phi_c, phi_d, bragglookup, omegalookup, lplookup, deltalookup, incl_junk=False)

Computes the $x$ and $y$ points of a simultaneous conjugate Ramsey-Borde interferometer (SCRBI) ellipse from the provided vectors of atom positions and velocities. In addition to the positions x0,y0,z0 and velocities vx,vy,vz the user must provide big $T$ as T (in seconds), $T'$ as Tp (in seconds), $n_\mathrm{Bragg}$ as n, and $N_\mathrm{Bloch}$ as N. The average population at each output port is then used to compute x and y.

In order to generate a single ellipse from a single function call the user can pass a vector of common mode phases phi_c. This common mode phase is then added to the appropriate wavefunctions before the calculation of the output port populations to produce an ellipse. The length of the returned x and y is given by the length of phi_c. The user may also pass in a value for phi_d, which sets the opening of the ellipse.

The user must provide the functions bragglookup, omegalookup, and lplookup. See below for more information on each of these functions.

Currently the function does not account for gravity in the calculation of z positions at each beamsplitter.

Currently the function accounts for the phase produced from the Bragg beamsplitter process (determined by calling bragglookup), and the local wavefront phase (determined by calling lplookup). Unless the user implements it manually in bragglookup or lplookup the global laser phase (i.e. $\omega t-kz$) will be ignored. Other ignored sources of phase include the free evolution phase, gravity gradient phase, and separation phases.

As it is difficult to program a completely general function for all of the scenarios that we might be interested in, it might be a good idea to open up the source code of this function and write a custom version for your use case!

Parameters:

x0 – Vector of initial atom x positions.
y0 – Vector of initial atom y positions.
z0 – Vector of initial atom z positions.
vx – Vector of initial atom x velocities.
vy – Vector of initial atom y velocities.
vz – Vector of initial atom z velocities.
T – Big $T$ in seconds.
Tp – $T'$ in seconds.
n – $n_\mathrm{Bragg}$.
N – $N_\mathrm{Bloch}$.
phi_c – Vector of common mode phases.
phi_d – Vector of differential mode phases.
bragglookup – Function that returns the phase difference between the provided initial and final states after undergoing a Bragg pulse. Must accept arguments ni, nf, omega, delta, where omega and delta can be vectors.
omegalookup – Function that returns the effective Rabi frequency at the specified location. Must accept arguments x, y, z, where all arguments are vectors of the same length.
lplookup – Function that returns the local laser phase at the specified location. Must accept arguments x, y, z, where all arguments are vectors of the same length. Note that this should return the kvector, and not the effective kvector as the mwave.simulation_utils.cloud_to_scrbi_ellipse_xy() function will apply a factor of two the the phase returned by lplookup.
deltalookup – Function that returns the detuning specified velocity. Must accept arguments v in whatever units are used by the user to define the cloud velocity spread. The returned detuning must be in units of the recoil frequency.
incl_junk – If True the junk ports are included in the calculation of x and y (i.e. x is calculated via (pA - pB)/(pA + pB) where pA is the population in output port A plus the population in the junk ports that exit with the same momentum as port A, same for pB, etc.). Defaults to False.

Returns:

The x and y values of the ellipse. These will have the same length as the provided vector phi_c.

import numpy as np
from mwave.precompute import load_fast_bragg_evaluator
from mwave.simulation_utils import cloud_to_scrbi_ellipse_xy
from alphautil.analysis import fit_ellipse_coeff
from alphautil.ellipse import get_sci_params
from matplotlib import pyplot as plt

np.random.seed(13703599)

npoints = 100000
x0 = np.random.randn(npoints)*0.75e-3
y0 = np.random.randn(npoints)*0.75e-3
z0 = np.random.randn(npoints)*0.75e-3
v0 = np.random.randn(npoints)*0.35e-3
vx = np.random.randn(npoints)*3.5e-3*1.8
vy = np.random.randn(npoints)*3.5e-3*1.8
T = 100e-3
Tp = 10e-3
phi_c = np.linspace(0,2*np.pi)
phi_d = np.pi/4
Omega0 = 32
w0 = 10e-3

n_init = 0
n_bragg = 5
N_bloch = 100

bragglookup = load_fast_bragg_evaluator('sig0.260.h5', n_init, n_bragg, N_bloch)

def omegalookup(x, y, z):
    return Omega0*np.exp(-2*(x**2 + y**2)/(w0**2))

# (l)ocal (p)hase (lookup)
def lplookup(x, y, z):
    wavelen = 852e-9
    zR = np.pi*w0**2/wavelen
    kk = 2*np.pi/wavelen
    return kk*(x**2 + y**2)/(2*zR)

def deltalookup(v, n_bragg):
    return 4*n_bragg + 4*(v/0.0035) # The modification to delta is 4 times the velocity defined in units of recoil velocities

x, y = cloud_to_scrbi_ellipse_xy(x0, y0, z0, vx, vy, v0, T, Tp, n_bragg, N_bloch, np.exp(1j*phi_c), np.exp(1j*phi_d), bragglookup, omegalookup, lplookup, deltalookup)

coeff = fit_ellipse_coeff(x, y)
bx, by, Ax, Ay, phi_d_fit = get_sci_params(coeff)

plt.scatter(x, y)
plt.ylim([-1,1])
plt.xlim([-1,1])
plt.gca().set_aspect('equal')
plt.show()

print('differential phase fit error=%0.3f mRad' % ((phi_d - phi_d_fit)*1e3))
print(Ax)
print(Ay)

Precompute module

exception mwave.precompute.NotPrecomputedError(path)

mwave.precompute.load_fast_bragg_evaluator(fname, n_init, n_bragg, N_bloch)

This function loads a function that quickly evaluates Bragg pulse precompute tables for SCRBI geometries on a grid of inputs using the scipy regular grid interpolator with the cubic method enabled. This is useful for simulating a atom cloud with transverse motion.

Parameters:

fname – The name of the HDF5 precompute table to load. phi Datasets are loaded using the read_bragg_precompute function. It is assumed phi is computed on a grid of (omega,delta). If this is not the case this function will return gibberish!
n_init – The initial momentum state, this is zero for most precompute tables.
n_bragg – The Bragg order.
N_bloch – The Bloch order (used when loading multifrequency pulses).

Returns:

A function fbe that takes parameters n0,nf,omega,delta, which are the initial and target momentum states, the value of omega, and the value of delta. These arguments must be supplied as equal length numpy arrays.

mwave.precompute.load_lookup_table(fname, kvec=None, method='cubic')

Loads a single lookup table where $n_0=0$, returns a function that interpolates this table using the provided method. The lookup table must be a 2D grid over $\Omega$ and $\delta$. If kvec is provided then the function aligns its k-vector with the one defined by kvec.

Also see the mwave.precompute.load_precomputed_gbragg() function.

Parameters:

fname – The filename to load the precompute table from. Internally the read is performed using read_bragg_precompute.
kvec – Optional. If provided only the values in kvec will be returned by the interpolation function.
method – The interpolation method to use. This is directly passed to scipy.interpolate.RegularGridInterpolator().

Returns:

A tuple containing (kvec, interpolation_function), where the interpolation_function is returned by scipy.interpolate.RegularGridInterpolator().

mwave.precompute.load_precomputed_gbragg(single_path, multi_path=None, method='cubic', table_sigma=None, table_modulation_frequency=None, flip_negatives=True)

Returns a function that provides similar functionality to the mwave.integrate.gbragg() function but uses a lookup table for faster evaluation.

For example usage of this function checks against equivalent direct computation methods see Using precompute tables.

Internally the function loads a lookup table where $n_0=0$ using the mwave.precompute.load_lookup_table() function. Then several frame transformations are performed to transform an arbitrary input state to make $n_0 o n_0’=0$. The user can optionally provide a multifrequency lookup table as well, in which case the returned function has the capability to return wavefunctions for both single and multifrequency Bragg pulses.

Parameters:

single_path – The name of the single frequency HDF5 precompute table to load. phi Datasets are loaded using the read_bragg_precompute function. It is assumed phi is computed on a grid of (omega,delta). If this is not the case this function will return gibberish!
multi_path – The name of the multi-frequency HDF5 precompute table to load. phi Datasets are loaded using the read_bragg_precompute function. It is assumed phi is computed on a grid of (omega,delta). If this is not the case this function will return gibberish!
method – The interpolation method to use. This is directly passed to scipy.interpolate.RegularGridInterpolator().
table_sigma – The value of sigma used in the lookup table. This is not directly checked against the precompute table contents, but it is used to check user input to the returned function.
table_modulation_frequency – The value of the modulation frequency used in the lookup table. This is not directly checked against the precompute table contents, but it is used to check user input to the returned function.
flip_negatives – If True the symmetry of the problem is used to interpolate negative provided values of delta to the returned function. This means that the precompute table does not need to include negative values of delta.

Returns:

A function that utilizes the referenced lookup tables to compute the effect of the Bragg pulse that takes arguments kvec, k0, sigma, omega, delta, delta_phase, mod_freq=None, mod_phase=0.0. See Using precompute tables for example usage.

mwave.precompute.read_bragg_precompute(fname, n0, nf, n_bragg, N_bloch=None)

Reads a precomputed Bragg dataset from an HDF5 file.

Parameters:

fname – The name of the HDF5 file to read from.
n0 – The initial momentum state of the Bragg process.
nf – The final momentum state of the Bragg process.
n_bragg – The Bragg order used.
N_bloch – Optional. The Bloch order used in the simulation. This will load a multifrequency simulation.

Returns:

A tuple containing phi, kvec, grid. For a description of grid see write_bragg_precompute.

mwave.precompute.write_bragg_precompute(fname, phi, kvec, grid, n0, nf, n_bragg, N_bloch=None)

Saves the precomputed dataset in an HDF5 file with the given name.

Parameters:

fname – The name of the HDF5 file to write to.
phi – The table of phi values to save.
kvec – The momentum space vector.
grid – A list of tuples describing the parameter grid over which phi was applied. Each tuple contains information about a grid axis over which phi was computed. The first element contains a vector of parameter values, the second element is the name associated with the parameter.
n0 – The initial momentum state.
nf – The final momentum state.
n_bragg – The Bragg order used.
N_bloch – Optional. The Bloch order used in the simulation. This implies a multifrequency simulation was used.