Sample quantum initial states commonly encountered in quantum phase space simulations of Bose fields, including those encountered in quantum optics and Bose-Einstein condensates.
Wigner and positive-P distributions are available, being the most useful for dynamical simulations.
Available distributions are glauberP, positiveP wigner, positiveW, husimiQ.
julia> ]add PhaseSpaceToolsjulia> using PhaseSpaceTools
help?> positiveP
search: positiveP positiveW
  α,α⁺ = positiveP(state <: State,N)
  Generate N samples from the positive-P (+P) phase-space distribution for state.
  Moments of the +P distribution generate quantum operator averages that are normally ordered.
  In general the two random variates α,α⁺ are statistically independent for the +P distribution. help?> State
search: State state estimate InvalidStateException AbstractSet AbstractVector AbstractVecOrMat stacktrace StackTraces istaskstarted abstract type AbstractRange AbstractPattern
  State
  Abstract supertype for all sampled states: state <: State.
  Examples
  ≡≡≡≡≡≡≡≡≡≡
  Find all states that may be sampled
  =====================================
  julia> subtypes(State)
  
  6-element Vector{Any}:
  Bogoliubov
  Coherent
  Crescent
  Fock
  Squeezed
  Thermal
  Create and sample a particular state (vacuum)
  ===============================================
  julia> s = Fock(0)
  Fock(0)
  julia> wigner(s,100)
  ┌ Warning: Fock state sampling for W is only valid for n ≫ 1.
  Here a warning is generated since Fock sampling is not well defined for small n. 
  A simpler way to sample the vacuum is 
  julia> s = Coherent(0) 
  Coherent(0.0 + 0.0im)  # type conversion to ComplexF64.
  
  julia> wigner(s,100)
  (ComplexF64[0.33820868828162637 + 0.4407579103538181im, 0.057183146091823775 - 0.2772571883006981im, ...
  generating two vectors of sampled points α,α⁺ in the complex plane. In this case, α = conj(α⁺), as we are not working with a doubled phase space.A coherent state |α⟩ is sampled as
α = 1.0+im*2.0 # coherent amplitude
s = Coherent(α) # define state |α⟩
N = 1000 # number of samples
a,a⁺ = positiveP(s,N)This is a special case where the two phase space variables a and a⁺ are complex conjugate, and non-stochastic in the +P representation.
An approximate Fock state sampler in the Wigner representation:
n = 100
s = Fock(n) # define number state |n⟩
N = 1000 # number of samples
a,a⁺ = wigner(s,N)Provides an approximate sampling of W that reproduces operator averages for large n.
See  /examples/PhaseSpaceTools.ipynb for more usage.
Numerical representation of quantum states in the positive-P and Wigner representations, 
M K Olsen, A S Bradley, 
Optics Communications 282, 3924 (2009)