lgssystem

Define and solve models for LGS systems

source

LGSSystem

def LGSSystem(
    system:str, # Name of the atomic system to be loaded from file
    fixed_params:dict={}, # Values for system parameters that will be held fixed
):

LGSSystem creates numerical models for the chosen laser guide star atomic system.

source

LGSSystem.builtins

def builtins():

List the available LGS systems:

LGSSystem.builtins()
['Na330',
 'Na589+569',
 'Na589+569+repump',
 'NaD1',
 'NaD1_Toy',
 'NaD2',
 'NaD2_Repump']

source

LGSSystem.diagram

def diagram(
    name:str | pathlib._local.Path, # Name of built-in LGS system definition, or directory of supplied system
    kind:str='Toy', # "Toy" to show level diagram and pumped transitions, "ToScale" to show hyperfine structure to scale, "NotToScale" to show unscaled hf structure
):

Draw a level diagram for the system of the specified kind.

LGSSystem.diagram("Na330", "Toy")

LGSSystem.diagram("Na330", "ToScale")

LGSSystem.diagram("Na330", "NotToScale")

source

LGSSystem.info

def info(
    name:str | pathlib._local.Path, # Name of built-in LGS system definition, or directory of supplied system
):

Display information about LGS system name.

3S

LGSSystem.info("NaD2_Repump")

Atomic levels

\text{3S}_{\frac{1}{2}}, \text{3P}_{\frac{3}{2}}

Pump transitions

  1. \text{3S}_{\frac{1}{2}}\land F=1\to \text{3P}_{\frac{3}{2}}
  2. \text{3S}_{\frac{1}{2}}\land F=2\to \text{3P}_{\frac{3}{2}}

Transition wavelengths

  • \text{3P}_{\frac{3}{2}}\to \text{3S}_{\frac{1}{2}}: 589.158 nm

Substructure

  • Hyperfine structure included
  • Zeeman structure included
  • 24 total sublevels

Density matrix elements

  • All populations included
  • All Zeeman coherences (between same level and same F) included
  • All hyperfine coherences (between same level and different F) neglected
  • Optical coherences (between different levels) included for pump transitions only
  • 374 total density matrix elements

Input parameters

  • BeamTransitRatePerS
  • BFieldG
  • MagneticAzimuthDegrees
  • MagneticZenithDegrees
  • RecoilParameter
  • SDampingCollisionRatePerS
  • TemperatureK
  • VccRatePerS
  • DetuningHz1
  • DetuningHz2
  • EllipticityDegrees1
  • EllipticityDegrees2
  • IntensitySI1
  • IntensitySI2
  • LaserWidthHz1
  • LaserWidthHz2
  • PolarizationAngleDegrees1
  • PolarizationAngleDegrees2

Load a “toy” (no angular momentum) model LGS system and fix values for most of the parameters:

lgs = LGSSystem(
    'NaD1_Toy',
    fixed_params={
      'DetuningHz1': 0,
      'LaserWidthHz1': 10.0e6,
      'BFieldG': 0.5,
      'BeamTransitRatePerS': 1e2,
      'VccRatePerS': 3e4,
      'TemperatureK': 2e2,
      'RecoilParameter': 1,
    }
)

The only remaining parameter is the laser intensity. Here the value 1 indicates that the parameter is a scalar quantity:

lgs.parameters
Parameters({IntensitySI1: 1})

Define a value for the intensity to use for the following examples:

mu = Mu(IntensitySI1=46.)

source

LGSSystem.operator

def operator(
    vg:pylgs.velocitygroups.VelocityGroups, # `VelocityGroups` specification
)->pymor.operators.interface.Operator:

The evolution operatior for the system, given the VelocityGroups vg.

Form the operator with the example VelocityGroups:

op = lgs.operator(vg)

The operator operates on vectors from the 2\times4 (velocity \times density matrix) product vector space:

op.source
Atomic velocity(2) ⛒ Density matrix(4) 
<xarray.DataArray 'XarrayVectorSpace' (Atomic velocity: 2, Density matrix: 4)> Size: 64B
array([[    0   ,     0   ,     0   ,     0   ],
       [    0   ,     0   ,     0   ,     0   ]])
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 16B -1.5 1.5
  * Density matrix   (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, ...

And returns a vector from the same space:

op.range
Atomic velocity(2) ⛒ Density matrix(4) 
<xarray.DataArray 'XarrayVectorSpace' (Atomic velocity: 2, Density matrix: 4)> Size: 64B
array([[    0   ,     0   ,     0   ,     0   ],
       [    0   ,     0   ,     0   ,     0   ]])
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 16B -1.5 1.5
  * Density matrix   (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, ...

Substitute the free parameter values into the operator and write it explicitly as a block matrix, with 2\times2 blocks for the velocity space, each block being 4\times4 on the density matrix space:

_mat = to_matrix(op, mu=mu).toarray()/1e6
_mat
array([[  0.015 ,     0   , -3.7e+01, -6.1e+01,  -0.015 ,     0   ,     0   ,     0   ],
       [    0   ,  6.2e+01, -6.1e+03,     0   ,     0   ,     0   ,     0   ,     0   ],
       [ 1.9e+01,  6.1e+03,  6.2e+01, -1.9e+01,     0   ,     0   ,     0   ,     0   ],
       [    0   ,     0   ,  3.7e+01,  6.1e+01,     0   ,     0   ,     0   ,  -0.015 ],
       [ -0.015 ,     0   ,     0   , -0.0016 ,   0.015 ,     0   , -3.7e+01, -6.1e+01],
       [    0   ,     0   ,     0   ,     0   ,     0   ,  6.2e+01,  6.1e+03,     0   ],
       [    0   ,     0   ,     0   ,     0   ,  1.9e+01, -6.1e+03,  6.2e+01, -1.9e+01],
       [    0   ,     0   ,     0   ,  -0.015 ,     0   ,     0   ,  3.7e+01,  6.1e+01]])

source

LGSSystem.rhs

def rhs(
    vg:pylgs.velocitygroups.VelocityGroups, # `VelocityGroups` specification
)->pymor.vectorarrays.interface.VectorArray:

The right-hand side vector for this system given the VelocityGroups vg.

Form the right-hand side VectorArray:

vec = lgs.rhs(vg)
vec
∅ ⛒ Atomic velocity(2) ⛒ Density matrix(4) 
<xarray.DataArray (Atomic velocity: 2, Density matrix: 4)> Size: 64B
array([[  5e+01 ,     0   ,     0   ,     0   ],
       [  5e+01 ,     0   ,     0   ,     0   ]])
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 16B -1.5 1.5
  * Density matrix   (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, ...

The vector is in the 2\times4 (velocity \times density matrix) product vector space:

vec.space
Atomic velocity(2) ⛒ Density matrix(4) 
<xarray.DataArray 'XarrayVectorSpace' (Atomic velocity: 2, Density matrix: 4)> Size: 64B
array([[    0   ,     0   ,     0   ,     0   ],
       [    0   ,     0   ,     0   ,     0   ]])
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 16B -1.5 1.5
  * Density matrix   (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, ...
vec.space.shape
[2, 4]

Write as a flattened 1D vector:

_arr = vec.to_numpy().T
_arr
array([[  5e+01 ,     0   ,     0   ,     0   ,   5e+01 ,     0   ,     0   ,     0   ]])

source

LGSSystem.stationary_model

def stationary_model(
    vg:Union[int, Iterable, pylgs.velocitygroups.VelocityGroups]=6, # Specification of velocity groups
    solver:pymor.solvers.interface.Solver | dict | None=None, # output='total_flux' # Operator to be applied to solution to obtain the output
)->pylgs.pymor.models.StationaryModel: # Model for the steady-state LGS system

Return a pyMOR StationaryModel for the steady state cw LGS system.

mu = Mu(IntensitySI1=46.)

An integer can be supplied for vg to get that many evenly spaced velocity groups:

model = lgs.stationary_model(vg=30)

Solve the model and calculate the return flux by supplying values for all free parameters:

sol = model.solve({'IntensitySI1': 1e3})
sol
∅ ⛒ Atomic velocity(30) ⛒ Density matrix(4) 
<xarray.DataArray (Atomic velocity: 30, Density matrix: 4)> Size: 960B
array([[ 2.6e-05, -1.9e-07,  -1e-09 ,  3.2e-09],
       [  8e-05 , -6.4e-07, -3.6e-09,  1.1e-08],
       [ 0.00022, -1.9e-06, -1.2e-08,  3.6e-08],
       [ 0.00058, -5.4e-06, -3.6e-08,  1.1e-07],
       [ 0.0014 , -1.4e-05,  -1e-07 ,  3.1e-07],
       [ 0.0031 , -3.5e-05, -2.8e-07,  8.3e-07],
       [ 0.0063 ,  -8e-05 , -7.2e-07,  2.1e-06],
       [  0.012 , -0.00017, -1.7e-06,  5.1e-06],
       [  0.021 , -0.00034, -4.1e-06,  1.2e-05],
       [  0.034 , -0.00065, -9.1e-06,  2.6e-05],
       [  0.05  , -0.0012 ,  -2e-05 ,  5.8e-05],
       [  0.068 , -0.0021 , -4.6e-05,  0.00013],
       [  0.087 , -0.0037 , -0.00011,  0.00032],
       [   0.1  , -0.0071 , -0.00036,   0.001 ],
       [  0.097 ,  -0.019 , -0.0029 ,  0.0081 ],
       [   0.1  ,   0.02  ,  -0.003 ,  0.0085 ],
       [  0.11  ,  0.0076 , -0.00039,  0.0011 ],
       [  0.087 ,  0.0037 , -0.00011,  0.00033],
       [  0.069 ,  0.0021 , -4.6e-05,  0.00013],
       [  0.05  ,  0.0012 ,  -2e-05 ,  5.8e-05],
       [  0.034 ,  0.00065, -9.1e-06,  2.6e-05],
       [  0.021 ,  0.00034, -4.1e-06,  1.2e-05],
       [  0.012 ,  0.00017, -1.7e-06,  5.1e-06],
       [ 0.0063 ,   8e-05 , -7.2e-07,  2.1e-06],
       [ 0.0031 ,  3.5e-05, -2.8e-07,  8.3e-07],
       [ 0.0014 ,  1.4e-05,  -1e-07 ,  3.1e-07],
       [ 0.00058,  5.4e-06, -3.6e-08,  1.1e-07],
       [ 0.00022,  1.9e-06, -1.2e-08,  3.6e-08],
       [  8e-05 ,  6.4e-07, -3.6e-09,  1.1e-08],
       [ 2.6e-05,  1.9e-07,  -1e-09 ,  3.2e-09]])
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 240B -2.9 -2.7 -2.5 ... 2.7 2.9
  * Density matrix   (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, ...
_flux = model.total_flux({'IntensitySI1': 1e3}).item()
_flux
96977.58974001158
model.integrated_flux({'IntensitySI1': 1e3})
∅ ⛒ Integral(1) ⛒ Transition(1) 
<xarray.DataArray (Integral: 1, Transition: 1)> Size: 8B
array([[ 9.7e+04]])
Coordinates:
  * Integral    (Integral) int64 8B 1
  * Transition  (Transition) <U33 132B '3P<sub>1/2</sub>→3S<sub>1/2</sub>'
Attributes:
    long_name:  Flux
    units:      photon/atom/s/sr

Solve the model and return the steady-state density matrix:

sol = model.solve({'IntensitySI1': 1e3})

The solution can be visualized in several ways. Plot the flux for each velocity bin:

model.flux_distribution(sol).visualize()

Plot the real and imaginary parts of the density-matrix elements for each velocity bin:

sol.visualize()

Plot just the level populations:

model.level_population(sol).visualize()

Plot the total atomic population:

model.population(sol).visualize()

Note that the above visualizations plot the contribution to the quantity from each velocity group. If the velocity groups are not all of the same width, this will cause the plots to appear distorted. In this case one can instead plot the velocity distribution of the quantities using flux_distribution, population_distribution, etc. To plot the velocity distribution of the density-matrix elements, use velocity_normalize(sol).

Table 1: Operators that can be applied to the solution of an LGS model
Description
population Fraction of total atomic population in each velocity group
total_population Total atomic population (should equal 1)
population_distribution Atomic population per unit velocity interval
level_population Population of each atomic level in each velocity group
total_level_population Total atomic population in each level
level_population_distribution Level population per unit velocity interval
flux Return flux on each transition from each velocity group
total_flux Total return flux
integrated_flux Integral of return flux over velocity groups
flux_distribution Return flux per unit velocity interval
velocity_sum Sum over velocity groups
velocity_normalize Normalize by velocity-group width

Check that the atomic population summed over all velocity groups is approximately 1:

_pop = model.total_population(sol).item()
_pop
0.9933768838135858
type(to_matrix(model.operator.assemble(Mu({'IntensitySI1': 1e3}))))
scipy.sparse._csr.csr_array

source

XarrayVectorArray.interp_like

def interp_like(
    other:pylgs.pymor.vectorarrays.XarrayVectorArray | pylgs.pymor.vectorarrays.XarrayVectorSpace
):

source

refine_solution

def refine_solution(
    old_sol, old_model, new_model
):

Refine the solution from an old model to a new model by interpolating the normalized solution onto the new model’s velocity grid.

source

transfer_ilu

def transfer_ilu(
    old_model, new_model
):

source

create_adaptive_model

def create_adaptive_model(
    model_func, solve_func, flux_func, vg, rtol, atol:int=0, refine_fraction:float=0.2, max_iterations:int=50,
    max_velocity_groups:int=500, use_initial_guess:bool=True, reuse_ilu:bool=True, return_intermediate:bool=False
):

source

LGSSystem.adaptive_stationary_model

def adaptive_stationary_model(
    mu, # Parameter values to use during velocity-group refinement
    vg:int=20, # Initial velocity groups
    solver:NoneType=None, # Solver to use for the linear system
    rtol:float=0.02, # Global error tolerance for the return flux integral
    atol:int=0, # Absolute tolerance for the return flux integral
    refine_fraction:float=0.2, # Fraction of velocity groups to refine at each iteration
    return_intermediate:bool=False, # Whether to return intermediate integration results
    max_iterations:int=50, # Maximum number of refinement iterations
    max_velocity_groups:int=500, # Maximum number of velocity groups allowed during refinement
    use_initial_guess:bool=True, # Whether to use the initial guess when solving the refined model
    reuse_ilu:bool=True, # Whether to reuse the ILU factorization from the previous model
)->pylgs.pymor.models.StationaryModel: # Steady-state LGS model with more velocity groups in regions that produce more return flux

Create a StationaryModel with narrower velocity bins in velocity regions with higher flux.

Create a model with adaptively refined velocity groups using a particular value for the intensity parameter during refinement:

model, sol, intermediate_results = lgs.adaptive_stationary_model({'IntensitySI1': 1e3}, return_intermediate=True, solver={'preconditioner': {'n_workers': 1}})
# import logging

# logger = logging.getLogger("pylgs.pymor.solvers")
# while logger:
#     print(
#         logger.name or "root",
#         [(type(h).__name__, id(h), h.level) for h in logger.handlers],
#         "propagate =", logger.propagate,
#     )
#     logger = logger.parent

Solve and find total flux:

_flux = model.total_flux({'IntensitySI1': 1e3}).item()
_flux
145695.62157053617
_flux = model.integrated_flux({'IntensitySI1': 1e3}).item()
_flux
146630.76017041213

Solve for the density matrix. The solution has density-matrix values for 17 velocity groups:

sol = model.solve({'IntensitySI1': 1e3})
sol
∅ ⛒ Atomic velocity(53) ⛒ Density matrix(4) 
<xarray.DataArray (Atomic velocity: 53, Density matrix: 4)> Size: 2kB
array([[ 5.6e-05, -4.2e-07, -2.3e-09,  7.2e-09],
       [ 0.00028, -2.3e-06, -1.4e-08,  4.4e-08],
       [ 0.0011 , -1.1e-05, -7.4e-08,  2.3e-07],
       [ 0.0039 , -4.3e-05, -3.4e-07,   1e-06 ],
       [  0.011 , -0.00015, -1.4e-06,  4.1e-06],
       [  0.028 , -0.00044,  -5e-06 ,  1.5e-05],
       [  0.056 , -0.0011 , -1.7e-05,  4.8e-05],
       [  0.042 , -0.0011 , -2.1e-05,  5.9e-05],
       [  0.053 , -0.0017 , -3.8e-05,  0.00011],
       [  0.063 , -0.0026 , -7.5e-05,  0.00021],
       [  0.072 , -0.0041 , -0.00017,  0.00048],
       [ 0.0021 , -0.00015, -7.8e-06,  2.2e-05],
       [ 0.0023 , -0.00017, -8.7e-06,  2.5e-05],
       [ 0.0047 , -0.00035, -1.9e-05,  5.3e-05],
       [ 0.0095 , -0.00074, -4.2e-05,  0.00012],
       [  0.019 , -0.0017 , -0.0001 ,  0.0003 ],
       [  0.039 , -0.0043 , -0.00035,  0.00099],
       [ 0.0087 , -0.0013 , -0.00014,  0.00039],
       [ 0.0088 , -0.0014 , -0.00018,  0.00051],
       [  0.017 , -0.0034 , -0.00056,  0.0016 ],
...
       [  0.013 ,   0.002 , -0.00026,  0.00073],
       [  0.012 ,  0.0017 , -0.00018,  0.00052],
       [  0.011 ,  0.0014 , -0.00014,  0.00039],
       [  0.011 ,  0.0012 , -0.00011,  0.0003 ],
       [  0.02  ,  0.0021 , -0.00015,  0.00044],
       [  0.02  ,  0.0017 , -0.00011,  0.00031],
       [  0.02  ,  0.0015 ,  -8e-05 ,  0.00023],
       [  0.037 ,  0.0024 , -0.00011,  0.0003 ],
       [  0.035 ,  0.0018 , -6.8e-05,  0.00019],
       [  0.033 ,  0.0015 , -4.6e-05,  0.00013],
       [  0.031 ,  0.0012 , -3.2e-05,   9e-05 ],
       [  0.053 ,  0.0017 , -3.8e-05,  0.00011],
       [  0.043 ,  0.0011 , -2.1e-05,  5.9e-05],
       [  0.056 ,  0.0012 , -1.7e-05,  4.8e-05],
       [  0.028 ,  0.00044,  -5e-06 ,  1.5e-05],
       [  0.011 ,  0.00015, -1.4e-06,  4.1e-06],
       [ 0.0039 ,  4.3e-05, -3.4e-07,   1e-06 ],
       [ 0.0011 ,  1.1e-05, -7.4e-08,  2.3e-07],
       [ 0.00028,  2.3e-06, -1.4e-08,  4.4e-08],
       [ 5.6e-05,  4.2e-07, -2.3e-09,  7.2e-09]])
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 424B -2.85 -2.55 ... 2.55 2.85
  * Density matrix   (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, ...

Plot the velocity distribution of flux. Each point marks the center of a velocity group:

model.flux_distribution(sol).visualize(markers=True)

Velocity distribution of the density matrix:

model.velocity_normalize(sol).visualize()

Of the total population:

model.population_distribution(sol).visualize()

Of the level populations:

model.level_population_distribution(sol).visualize()

source

LGSSystem.instationary_model

def instationary_model(
    vg:Union[int, Iterable, pylgs.velocitygroups.VelocityGroups]=6, # Specification of velocity groups
    T:float=1e-06, # Final time (seconds)
    num_values:int=100, # Number of time step values to return
    initial_rho:NoneType=None, # Initial density matrix. `None` means to use the thermal distribution
    time_stepper:NoneType=None, # Solver to use. `None` means to use the default `BDFTimeStepper`
    preconditioner:NoneType=None, **kwargs
)->pylgs.pymor.models.InstationaryModel: # Model for the LGS system dynamics

Create an InstationaryModel for the LGS system dynamics.

An instationary model with 30 velocity groups:

vg=30
num_values=100
T=1.e-4
inst_model = lgs.instationary_model(
    vg=vg, 
    T=T, 
    num_values=num_values
)

The BDF solver is efficient but can be very sensitive to the absolute and relative tolerance values. Set appropriate values for this problem:

pymor.basic.set_defaults({
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_atol': 1e-3,
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_rtol': 1e-6
})

Solve for a constant value of intensity:

inst_sol = inst_model.solve({'IntensitySI1': 1e3})

Visualize the total flux as a function of time:

inst_model.total_population(inst_sol).visualize().show()
inst_model.total_flux(inst_sol).visualize().show()

Plot the density matrix elements as a function of atomic velocity and animate as a function of time:

inst_model.velocity_normalize(inst_sol).visualize().update_layout(sliders=[dict(active=4)]).show()

Animate the population distribution as a function of time:

inst_model.population_distribution(inst_sol).visualize().show()

Animate the level population distributions:

inst_model.level_population_distribution(inst_sol).visualize().show()

Set a shorter evolution time for the model:

T=1e-7
num_values=200
inst_model = inst_model.with_(T=T, num_values=num_values)

Solve the model with a modulated intensity parameter:

inst_mu = {'IntensitySI1': "5000.*sin(1.e8*t)**2"}
pymor.basic.set_defaults({
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_atol': 3e-4,
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_rtol': 1e-7
})
inst_sol = inst_model.solve(inst_mu)

Visualize the total flux as a function of time:

inst_plot = inst_model.total_flux(inst_sol).visualize()
inst_plot.show()

Visualize the density matrix elements:

inst_model.velocity_normalize(inst_sol).visualize().show()

Visualize the flux distribution:

inst_model.flux_distribution(inst_sol).visualize().show()
inst_model.population_distribution(inst_sol).visualize().show()

The numpy model is currently much faster for explicit methods than the xarray model. However, it requires Kronecker-expanding the DM x VG product-space operators. It might be a little faster if this wasn’t necessary, but no more than a factor of 2 or 3. It seems like BlockOperator and BlockVectorArray are designed for this, but it looks like they are not optimized for speed. It could be implemented by optimizing XarrayMatrixOperator so that all apply operations can bypass xarray and use numpy instead. Note that sparse DataArrays currently use COO format, would want to switch to CSR, but need to enforce 2D matrices or allow N-D CSR arrays.

Evaluating the coefficients at a given time also could be sped up a little (use lambdify to convert ExpressionFunctions to GenericFunctions), but it is not the dominant source of computation time.

# def to_dense_matrix(op):
#     a = to_matrix(op)
#     try: a = a.todense()
#     except: pass
#     return a

source

LGSSystem.numpy_instationary_model

def numpy_instationary_model(
    vg, *, T:float=1e-06, # Final time (seconds)
    num_values:int=100, # Number of time step values to return
    initial_rho:NoneType=None, # Initial density matrix. `None` means to use the thermal distribution
    time_stepper:NoneType=None, # Solver to use. `None` means to use the default `BDFTimeStepper`
    preconditioner:NoneType=None
): # Model for the LGS system dynamics

Create an InstationaryModel for the LGS system dynamics.

npin_model = lgs.numpy_instationary_model(
    vg=vg, 
    T=T, 
    num_values=num_values,
    time_stepper=AdamsTimeStepper()
)

Solve the model with a modulated intensity parameter:

npin_sol = npin_model.output(inst_mu)
plt.plot(npin_sol.sum(0))

plt.plot(inst_model.total_flux(inst_sol).to_numpy()[0])

plt.plot(npin_sol.T);

plt.plot(npin_sol.sum(1));

source

transpose_vector

def transpose_vector(
    v
):

source

transpose_operator

def transpose_operator(
    op
):

source

transpose_floq_model

def transpose_floq_model(
    floq_model
):

source

LGSSystem.stationary_floquet_model

def stationary_floquet_model(
    mu, # Parameter specification for parameters modulated at a frequency `nu`
    time, # List of time values to use for reconstruction
    n_vector, # Number of terms in the system Fourier series
    nu:str='nu', # Name of the frequency variable
    n_operator:NoneType=None, # Number of terms in the operator Fourier series. `None` defaults to `n_vector` - 1
    vg:Union[int, Iterable, pylgs.velocitygroups.VelocityGroups]=6, # Specification of velocity groups
    solver:NoneType=None
)->pylgs.pymor.models.StationaryFloquetModel: # Model for the Fourier coefficients of the periodic LGS system

Return a StationaryFloquetModel for the Fourier coefficients of the periodic state of a modulated LGS system.

The StationaryFloquetModel solves for the periodic state of a modulated system after the transient dynamics have died out.

Build a stationary Floquet model for six Fourier harmonics of a system with sinusoidally modulated intensity:

floq_model = lgs.stationary_floquet_model(
    mu={'IntensitySI1': "5000.*sin(nu*t)**2"},
    time=np.linspace(0, T, num_values),
    n_vector=10,
    vg=30
)

Define values for the modulation frequency and light detuning:

Solve for the periodic state. The solution contains 7 Fourier coefficients for each density-matrix element in each velocity group:

floq_mu = Mu(nu=1e8)
floq_sol = floq_model.solve(floq_mu)
floq_plot = floq_model.reconstructed_total_flux(floq_sol, floq_mu).real.visualize()
floq_plot
floq_model.reconstructed_flux_distribution(floq_sol, floq_mu).real.visualize()

The real and imaginary parts of the Fourier coefficients of the time-dependent return flux:

floq_model.total_flux(floq_sol).real.visualize()
floq_model.reconstruct(floq_sol, floq_mu).visualize()

The total return flux reconstructed as a function of time:

The real and imaginary parts of the Fourier coefficients of the real and imaginary parts of the density-matrix elements:

floq_model.velocity_normalize(floq_mu).visualize()

The real and imaginary parts of the density-matrix elements reconstructed as a function of time:

floq_model.velocity_normalize(floq_model.reconstruct(floq_sol, floq_mu)).real.visualize()

source

LGSSystem.adaptive_stationary_floquet_model

def adaptive_stationary_floquet_model(
    mu, # Parameter specification for parameters modulated at a frequency `nu`
    time, # List of time values to use for reconstruction
    n_vector, # Number of terms in the system Fourier series
    mu_refine:NoneType=None, # Additional parameter specifications used during velocity-group refinement
    nu:str='nu', # Name of the frequency variable
    n_operator:NoneType=None, # Number of terms in the operator Fourier series. `None` defaults to `n_vector` - 1
    vg:int=20, # Initial (number of) velocity groups
    solver:NoneType=None, # Steady-state solver or options for the default solver
    rtol:float=0.02, # Global error tolerance for the return flux integral
    refine_fraction:float=0.2, # Fraction of velocity groups to refine at each iteration
    max_iterations:int=50, # Maximum number of refinement iterations
    max_velocity_groups:int=500, # Maximum number of velocity groups allowed during refinement
    return_intermediate:bool=False, # Whether to return intermediate results from the refinement process
    atol:int=0, # Absolute error tolerance for the return flux integral
)->pylgs.pymor.models.StationaryFloquetModel: # Steady-state LGS model with more velocity groups in regions that produce more return flux

Create a StationaryModel with narrower velocity bins in velocity regions with higher flux.

Create a model with adaptively refined velocity groups using a particular value for the intensity parameter during refinement:

mu = {'nu': 1e8}
model, sol = lgs.adaptive_stationary_floquet_model(
    mu={'IntensitySI1': "5000.*sin(nu*t)**2"},
    mu_refine=mu,
    time=np.linspace(0, T, num_values),
    n_vector=20,
    rtol=.04,
    solver={'rtol': .01}
)
model.reconstructed_total_flux(sol, mu).real.visualize()
model.reconstructed_flux_distribution(sol, mu).real.visualize(markers=True)