Add logging statements for classification and ordering

[oameye]

using HarmonicBalance

using LinearAlgebra, StatsBase
omega1 = 2π * 0.576635 * 1e6 # 2π*0.5763*1e6 # Hz
omega2 = 2π * 1.744435 * 1e6 # 2π*1.74107*1e6 # Hz
gamma1 = 89.22# 89.22 # Hz
gamma2 = 537 # Hz
alpha1111 = 1.40e24 # 1.26e24 # m² x Hz²
alpha2222 = 2.8e25 # 3.25e25 # m² x Hz² # uncertain
alpha2111 = -2e22 # m² x Hz² # uncertain
alpha1122 = 2.92e23 # m² x Hz²
eta1 = 8e14 # m² x Hz
eta2 = 2e16 # 3.97e16 # m² x Hz # uncertain
F = 500 # m x Hz²

omega = mean( 2*pi*range(581.2,582.1,30)*1e3)

@variables ω₁ ω₂ ω γ₁ γ₂ G t x₁(t) x₂(t)
@variables A₁ A₂ A₃ A₄ ζ₁ ζ₂

# eq1_no_coupling = d(d(x₁, t),t) + (ω₁^2/ω^2)*x₁ + x₁^3 + (γ₁/ω)*d(x₁,t) + ζ₁*ω*x₁^2*d(x₁,t)
# eq2_no_coupling = d(d(x₂, t),t) + (ω₂^2/ω^2)*x₂ + x₂^3 + (γ₂/ω)*d(x₂,t) + ζ₂*ω*x₂^2*d(x₂,t)
eq1_coupling = 3*A₁*x₁^2*x₂ + A₂*x₂^2*x₁
eq2_coupling = A₃*x₁^3 + A₄*x₁^2*x₂

eq1 = d(d(x₁, t),t) ~ - (ω₁^2/ω^2)*x₁ - x₁^3 - (γ₁/ω)*d(x₁,t) - ζ₁*ω*x₁^2*d(x₁,t) - eq1_coupling + (G/ω^3)*cos(t)
eq2 = d(d(x₂, t),t) ~ - (ω₂^2/ω^2)*x₂ - x₂^3 - (γ₂/ω)*d(x₂,t) - ζ₂*ω*x₂^2*d(x₂,t) - eq2_coupling
eq_HB = [eq1, eq2]

fixed = Dict(
    ω₁ => omega1,
    ω₂ => omega2,
    γ₁ => gamma1,
    γ₂ => gamma2,
    A₁ => alpha2111/√(alpha1111*alpha2222),
    A₂ => alpha1122/alpha2222,
    A₃ => alpha2111*√alpha2222/(alpha1111*√alpha1111),
    A₄ => alpha1122/alpha1111,
    ζ₁ => eta1/alpha1111,
    ζ₂ => eta2/alpha2222,
    G => F*√alpha1111
)


system_HB = DifferentialEquation(eq_HB, [x₁,x₂]);
add_harmonic!(system_HB, x₁, 1) # x will rotate at ω
add_harmonic!(system_HB, x₁, 3) # x will rotate at ω
add_harmonic!(system_HB, x₂, 1) # x will rotate at ω
add_harmonic!(system_HB, x₂, 3); # x₂ will rotate at 3*ω

harmonics_HB = get_harmonic_equations(system_HB)

res = 500
varied1 = ω => 2*pi*range(581,582.1,res)*1e3

result_HB = get_steady_states(harmonics_HB, varied1, fixed, threading=true);