using SLHQuantumSystems
using SecondQuantizedAlgebra
using Symbolics
using ControlSystems
using GLMakie

mode = MechanicalMode("")

@variables Ω m Γ
b = Destroy(FockSpace(:mass), operatornames(mode)[1])
paramdict = Dict(zip(nameof.([Ω, m, Γ]), [Ω, m, Γ]))

opdict = Dict(getfield(b, :name) => b)

slh = SLH("mass", [mode], paramdict, opdict, ["in"], ["out"], [1], [Γ * b], Ω * b' * b)

ss = QuantumStateSpace(slh)
qss = toquadrature(ss)

paramdict = Dict(
    Ω => 10, #10Hz resonance
    m => 4, #kg
    Γ => 1.0e-23
)

numeric = substitute(qss, paramdict)


freq_hz = collect(logrange(0.1, 20_000.0, 500))
freq = 2π .* freq_hz
500-element Vector{Float64}:
      0.6283185307179586
      0.6438773900282578
      0.6598215286056845
      0.6761604870027125
      0.6929040420215935
      0.7100622125645388
      0.7276452656287665
      0.7456637224500013
      0.7641283647981058
      0.7830502414286066
      ⋮
 103329.28189307102
 105887.99324887247
 108510.0651901959
 111197.066695812
 113950.60559663222
 116772.32953779229
 119663.92696455984
 122627.12813265534
 125663.70614359173

output noise spectral density

G = freqresp(numeric, freq)

sd = spectral_density(numeric, freq)
SpectralDensityMatrix: 2×2 over 500 frequencies (0.1 – 20000.0 Hz)
  Channels: out_x, out_p

Plotting

fig = Figure()
ax_tf = Axis(
    fig[1, 1];
    xscale = log10, yscale = log10,
    xlabel = "Frequency [Hz]",
    title = "transfer functions (unitless)"
)

for (ii, iname) in enumerate(sd.names)
    for (jj, jname) in enumerate(sd.names)
        lines!(ax_tf, freq_hz, abs.(G[ii, jj, :]), label = "$iname, $jname")
    end
end
fig[1, 2] = Legend(fig, ax_tf)


ax_noise = Axis(
    fig[2, 1];
    xscale = log10, yscale = log10,
    limits = (nothing, nothing, 10.0e-3, 10.0e3),
    xlabel = "Frequency [Hz]",
    title = "spectral densities (quanta per root hertz)"
)

for ii in sd.names
    for jj in sd.names
        lines!(ax_noise, freq_hz, abs.(sd[ii, jj]), label = "$ii, $jj")
    end
end

fig[2, 2] = Legend(fig, ax_noise)

fig
Example block output

This page was generated using Literate.jl.