# differential equation
function swing_dynamics!(dₜx, x, p, t)

    # unpack
    G,P,K,D = p.grid, p.power, p.coupling, p.damping 
    N = nv(G); ϕ = x[1:N]; ω = x[N+1:2N]

    # second-order, frequency is derivative of phase
    dₜx[1:N] .= ω;

    # swing equation
    L = laplacian_matrix(G)
    s = sin.( L * ϕ )
    dₜx[N+1:2N] .= @. P - D * ω - K * s 
end


function sync_rhs(x,p)
    N = length(x)÷2
    
    # time derivative at t = 0
    dₜx = similar(x)
    swing_dynamics!(dₜx, x, p, 0.0)
    
    # ignore constant angle shift 
    dₜx[1:N] .-= mean(dₜx[1:N])

    return dₜx[2:2N]
end

# steady state 
function synchronous_state(p)
    nlp = NonlinearProblem(sync_rhs, zeros(2*nv(p.grid)), p)
    sync = solve(nlp, NewtonRaphson(), reltol=1e-6).u
    return sync
end