Squeezing example

[oameye]

Squeezing can be measured in terms of eigenvectors. Following, Huber et al. it would be $\frac{(1-\alpha)}{(1+\alpha)}$ where $\alpha_{i,k} = 2 \left( \text{Im}[\tilde{c}^{u_i}_k] \text{Re}[\tilde{c}^{v_i}_k] - \text{Re}[\tilde{c}^{u_i}_k] \text{Im}[\tilde{c}^{v_i}_k] \right)$.

@variables α, ω, ω0, F, γ, t, x(t); # declare constant variables and a function x(t)

# define ODE
diff_eq = DifferentialEquation(d(x, t, 2) + ω0 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x)

# specify the ansatz x = u(T) cos(ω*t) + v(T) sin(ω*t)
add_harmonic!(diff_eq, x, ω)

# implement ansatz to get harmonic equations
harmonic_eq = get_harmonic_equations(diff_eq)

fixed = (α => 1, ω0 => 1.0, γ => 0.005, F => 0.002)   # fixed parameters
varied = ω => range(0.95, 1.05, 100)           # range of parameter values
result = get_steady_states(harmonic_eq, varied, fixed)

plot(result, x="ω", y="sqrt(u1^2 + v1^2)")

plot(
    plot_eigenvalues(result, branch=1),
    plot_eigenvalues(result, branch=1, type=:real, ylims=(-0.003, 0)),
    plot_eigenvalues(result, branch=2),
    plot_eigenvalues(result, branch=2, type=:real, ylims=(-0.003, 0)),
)

using LinearAlgebra
res = result
class = "stable"
branch = 1

filter = HB._get_mask(res, class)
filter_branch = map(x -> getindex(x, branch), replace.(filter, 0 => NaN))

x = string(first(keys(res.swept_parameters)))
varied = Vector{Float64}(collect(first(values(res.swept_parameters))))

eigenvalues = [
    eigvals(res.jacobian(get_single_solution(res; branch=branch, index=i))) for
    i in eachindex(varied)
]
eigenvalues_filtered = map(.*, eigenvalues, filter_branch)

eigenvectors = [
    eigvecs(res.jacobian(get_single_solution(res; branch=branch, index=i))) for
    i in eachindex(varied)
]
eigvecs_filtered = map(.*, eigenvectors, filter_branch);

vecs = eachcol(eigvecs_filtered[100])
function squeeze(v)
    alpha = 2 * (imag(v[1]) * real(v[2]) - real(v[1]) * imag(v[2]))
    alpha > 0 ? (1 - alpha) / (1 + alpha) : (1 + alpha) / (1 - alpha)
end
tmp = [squeeze.(eachcol(mat))[1] for mat in eigvecs_filtered]
plot(tmp)