Return noise in stochastic solvers

[gautierronan]

When solving Stochastic Schrödinger or Master equations, it is often interesting to know the time-dependent noise that was generated by the solver. For instance when simulating homodyne or heterodyne detection, the measurement signal for the $n$-th stochastic jump operator $C_n$ reads
$$\mathrm{d}Y_n = \sqrt{\eta}\langle C_n + C_n^\dagger\rangle \mathrm{d}t + \xi_n(t)$$
(Using the same notations as the documentation.)
While it is possible to get $\langle C_n + C_n^\dagger\rangle$ from the stochastic.master and stochastic.schrodinger solvers, I was not able to return $\xi_n(t)$. Is this possible to implement this feature?

[gautierronan]

Below is a simple workaround for anyone interested.

using DifferentialEquations
using QuantumOptics
using PyPlot

# Predefine noise
W = RealWienerProcess(0.0, 0.0, save_everystep=true)
​
# Run simulation
tout, rhot = stochastic.master(timespan, rho0, H, J, C; dt=dt, noise=W)
exp = expect(M, rhot)
​
# Plot expectation value and noise against time
figure()
plot(tout, exp)
plot(W.t, W)
display(gcf())