Constrain to positive angular velocity

tkf · Jun 25, 2018 · 08ecf96 · 08ecf96
1 parent 11d65cc
commit 08ecf96
Showing 1 changed file with 13 additions and 1 deletion.
diff --git a/src/codim2/diffeq.jl b/src/codim2/diffeq.jl
@@ -111,12 +111,24 @@ imag_eigval(prob) = _imag_eigval(eltype(prob.v0), prob.w0)
 _imag_eigval(::Type{T}, ::Real) where {T <: Real} = SVector{0, T}()
 _imag_eigval(::Type{<: Complex}, w0::Real) = SVector(w0)
 
-function isindomain(u, cache::DiffEqCodim2BifurcationCache)
+isindomain(u, cache::DiffEqCodim2BifurcationCache) =
+    _isindomain(u, cache, eltype(cache.prob.v0))
+
+function _isindomain(u, cache, ::Type{<: Real})
     t = SVector(u[end-1], u[end])
     tmin, tmax = cache.prob.t_domain
     return all(tmin .<= t .<= tmax)
 end
 
+function _isindomain(u, cache, ::Type{<: Complex})
+    # Constrain to positive angular velocity (w = u[end-2])
+    p = SVector(u[end-2], u[end-1], u[end])
+    tmin, tmax = cache.prob.t_domain
+    pmin = SVector(zero(tmin[1]), tmin[1], tmin[2])
+    pmax = SVector(typemax(tmax[1]), tmax[1], tmax[2])
+    return all(pmin .<= p .<= pmax)
+end
+
 residual!(H, u, cache::DiffEqCodim2BifurcationCache) =
     _residual!(H, u, cache.prob, cache.augsys_cache, statekind(cache.prob))