Hello everyone!
I am new to bifurcation analysis and I am trying to work with BifurcationKit.jl for my Master’s project.
The system I’m working on consists of a freely moving excentric mass coupled to a 2-dof aeroelastic structure that undergoes Hopf bifurcation. The eccentric mass can perform full revolutions (unbounded phase portrait but of periodic nature from 0 to 2π, similar to a pendulum). A time integration example is presente using DifferentialEquations.jl for a full revolution regime, in which θ corresponds to the excentric mass angle (Figure 1).
Continuation of equilibria works well for the system, determining the Hopf bifurcation point without issues (Figure 2). However, when I try to continue the periodic orbits, it only retrieves a residual orbit in the order of 1e-14 amplitude. So far I’ve tried PeriodicOrbitTrapProblem() and PeriodicOrbitOCollProblem(), and haven’t really explored Shooting methods.
I suspect the problem lies in the 2π-periodic nature of the system, which might be affecting the continuation. Full revolutions of the excentric mass seems to be periodic, but at f(t + T) = f(t) + k*2π. Has anyone dealt with this type of issue before? I’d appreciate any advice on how to handle it, or maybe it’s something simpler that I’m overlooking?
A working example of the code I’ve been working with is below.
P.S.: The vector field definition is huge due to the highly coupled behavior of the second derivatives, but it works well for time integration.
using Plots
using BifurcationKit
# vector field
function TS!(Y, p)
(;U, m, a, b, ra, xa, ωh, ωα, H, ρ, c0, c1, c2, c3, c4, mr, rr, dr, ln) = p
ξ, ξd, α, αd, θ, θd, xb, xbd = Y
[
ξd
-(U*αd*b^4*ρ^2*π^2 + 8*b*m^2*ωh^2*ra^2*ξ + 8*b*dr^2*m^2*mr*ωh^2*ξ - 8*α*b*m^2*ωα^2*ra^2*xa + 2*U^2*α*b^3*c0*ρ^2*π^2 - 2*U^2*α*b^3*c1*ρ^2*π^2 - 2*U^2*α*b^3*c3*ρ^2*π^2 + π*b^3*m*ωh^2*ρ*ξ + 8*αd^2*b*dr^3*m^2*mr^2*sin(α) + 4*U*a*αd*b^4*ρ^2*π^2 + U*αd*b^4*c0*ρ^2*π^2 - U*αd*b^4*c1*ρ^2*π^2 - U*αd*b^4*c3*ρ^2*π^2 + 2*U*b^4*c0*ρ^2*ξd*π^2 - 2*U*b^4*c1*ρ^2*ξd*π^2 - 2*U*b^4*c3*ρ^2*ξd*π^2 - 8*αd^2*b*dr^3*m^2*mr^2*sin(α - θ)^2*sin(α) - 6*U*a*αd*b^4*c0*ρ^2*π^2 + 6*U*a*αd*b^4*c1*ρ^2*π^2 + 6*U*a*αd*b^4*c3*ρ^2*π^2 - 8*U*a*b^4*c0*ρ^2*ξd*π^2 + 8*U*a*b^4*c1*ρ^2*ξd*π^2 + 8*U*a*b^4*c3*ρ^2*ξd*π^2 + 8*b*dr^2*m^2*mr^2*rr*θd^2*cos(θ) + 8*π*a^2*b^3*m*ωh^2*ρ*ξ - 8*b*dr^2*m^2*mr*ωh^2*ξ*sin(α - θ)^2 - 8*U^2*a*α*b^3*c0*ρ^2*π^2 + 8*U^2*a*α*b^3*c1*ρ^2*π^2 + 8*U^2*a*α*b^3*c3*ρ^2*π^2 + 8*U*a^2*αd*b^4*c0*ρ^2*π^2 - 8*U*a^2*αd*b^4*c1*ρ^2*π^2 - 8*U*a^2*αd*b^4*c3*ρ^2*π^2 - 4*π*U*αd*b^2*m*ρ*xa + 2*U^2*b^3*c1*c2*ρ^2*xbd*π^2 + 2*U^2*b^3*c3*c4*ρ^2*xbd*π^2 - 8*H*α^3*b*m^2*ωα^2*ra^2*xa - 8*αd^2*b*dr^3*m^2*mr^2*sin(θ)*cos(α - θ) + 8*αd^2*b*dr*m^2*mr*ra^2*sin(α) + 8*π*U*αd*b^2*m*ra^2*ρ + 8*b*m^2*mr*ra^2*rr*θd^2*cos(θ) - 8*b*ln*m^2*mr*ra^2*rr*θd*sin(θ) + 8*b*dr*m^2*mr*rr*θd^2*xa*cos(α - θ) - 8*b*dr^2*m^2*mr^2*rr*θd^2*cos(α)*cos(α - θ) + 8*α*b*dr*m^2*mr*ωα^2*ra^2*cos(α) + 8*π*U*αd*b^2*dr^2*m*mr*ρ + 16*π*U^2*α*b*c0*m*ra^2*ρ - 16*π*U^2*α*b*c1*m*ra^2*ρ - 16*π*U^2*α*b*c3*m*ra^2*ρ + 8*π*U*αd*b^2*c0*m*ra^2*ρ - 8*π*U*αd*b^2*c1*m*ra^2*ρ - 8*π*U*αd*b^2*c3*m*ra^2*ρ + 16*π*U*b^2*c0*m*ra^2*ρ*ξd - 16*π*U*b^2*c1*m*ra^2*ρ*ξd - 16*π*U*b^2*c3*m*ra^2*ρ*ξd - 8*b*dr^2*ln*m^2*mr^2*rr*θd*sin(θ) - 8*b*dr^2*m^2*mr^2*rr*θd^2*sin(α - θ)^2*cos(θ) - 8*αd^2*b*dr^3*m^2*mr^2*sin(α - θ)*cos(α)*cos(α - θ) + 8*π*a*α*b^3*m*ωα^2*ra^2*ρ + π*αd^2*b^3*dr*m*mr*ρ*sin(α) - 8*αd^2*b*dr*m^2*mr*ra^2*sin(θ)*cos(α - θ) + π*b^3*m*mr*ρ*rr*θd^2*cos(θ) + 8*π*U*a*αd*b^2*m*ρ*xa + 8*π*U^2*α*b*c0*m*ρ*xa - 8*π*U^2*α*b*c1*m*ρ*xa - 8*π*U^2*α*b*c3*m*ρ*xa + 4*π*U*αd*b^2*c0*m*ρ*xa - 4*π*U*αd*b^2*c1*m*ρ*xa - 4*π*U*αd*b^2*c3*m*ρ*xa - 8*U^2*a*b^3*c1*c2*ρ^2*xbd*π^2 - 8*U^2*a*b^3*c3*c4*ρ^2*xbd*π^2 + 2*U^3*b^2*c1*c2*c4*ρ^2*xb*π^2 + 2*U^3*b^2*c2*c3*c4*ρ^2*xb*π^2 + 8*αd^2*b*dr^2*m^2*mr*xa*sin(α - θ)*cos(α - θ) + 8*π*U*b^2*c0*m*ρ*xa*ξd - 8*π*U*b^2*c1*m*ρ*xa*ξd - 8*π*U*b^2*c3*m*ρ*xa*ξd + 8*b*dr*ln*m^2*mr*rr*θd*xa*sin(α - θ) - 8*b*dr^2*ln*m^2*mr^2*rr*θd*sin(α - θ)*cos(α) + 16*π*U^2*a*α*b*c0*m*ρ*xa - 16*π*U^2*a*α*b*c1*m*ρ*xa - 16*π*U^2*a*α*b*c3*m*ρ*xa - 8*b*dr^2*m^2*mr^2*rr*θd^2*sin(α - θ)*sin(θ)*cos(α - θ) + 8*α*b*dr*m^2*mr*ωα^2*ra^2*sin(α - θ)*sin(θ) - 8*U^3*a*b^2*c1*c2*c4*ρ^2*xb*π^2 - 8*U^3*a*b^2*c2*c3*c4*ρ^2*xb*π^2 + 16*π*U*a*b^2*c0*m*ρ*xa*ξd - 16*π*U*a*b^2*c1*m*ρ*xa*ξd - 16*π*U*a*b^2*c3*m*ρ*xa*ξd + 8*π*H*a*α^3*b^3*m*ωα^2*ra^2*ρ + 8*π*U^2*b*c1*c2*m*ρ*xa*xbd + 8*π*U^2*b*c3*c4*m*ρ*xa*xbd + 8*π*U^3*c1*c2*c4*m*ρ*xa*xb + 8*π*U^3*c2*c3*c4*m*ρ*xa*xb + 8*π*a^2*αd^2*b^3*dr*m*mr*ρ*sin(α) + 16*π*U^2*α*b*c0*dr^2*m*mr*ρ - 16*π*U^2*α*b*c1*dr^2*m*mr*ρ - 16*π*U^2*α*b*c3*dr^2*m*mr*ρ + 8*π*U*αd*b^2*c0*dr^2*m*mr*ρ - 8*π*U*αd*b^2*c1*dr^2*m*mr*ρ - 8*π*U*αd*b^2*c3*dr^2*m*mr*ρ + 8*π*a^2*b^3*m*mr*ρ*rr*θd^2*cos(θ) - 16*π*U*a*αd*b^2*c0*m*ra^2*ρ + 16*π*U*a*αd*b^2*c1*m*ra^2*ρ + 16*π*U*a*αd*b^2*c3*m*ra^2*ρ + 4*π*U*αd*b^2*dr*m*mr*ρ*cos(α) - 8*π*U*αd*b^2*dr^2*m*mr*ρ*sin(α - θ)^2 - 16*π*U*a^2*αd*b^2*c0*m*ρ*xa + 16*π*U*a^2*αd*b^2*c1*m*ρ*xa + 16*π*U*a^2*αd*b^2*c3*m*ρ*xa + 16*π*U*b^2*c0*dr^2*m*mr*ρ*ξd - 16*π*U*b^2*c1*dr^2*m*mr*ρ*ξd - 16*π*U*b^2*c3*dr^2*m*mr*ρ*ξd + 16*π*U^2*b*c1*c2*m*ra^2*ρ*xbd + 16*π*U^2*b*c3*c4*m*ra^2*ρ*xbd + 16*π*U^3*c1*c2*c4*m*ra^2*ρ*xb + 16*π*U^3*c2*c3*c4*m*ra^2*ρ*xb - π*αd^2*b^3*dr*m*mr*ρ*sin(θ)*cos(α - θ) + 8*H*α^3*b*dr*m^2*mr*ωα^2*ra^2*cos(α) - π*b^3*ln*m*mr*ρ*rr*θd*sin(θ) - 8*π*U^2*α*b*c0*dr*m*mr*ρ*cos(α) + 8*π*U^2*α*b*c1*dr*m*mr*ρ*cos(α) + 8*π*U^2*α*b*c3*dr*m*mr*ρ*cos(α) - 4*π*U*αd*b^2*c0*dr*m*mr*ρ*cos(α) + 4*π*U*αd*b^2*c1*dr*m*mr*ρ*cos(α) + 4*π*U*αd*b^2*c3*dr*m*mr*ρ*cos(α) - 16*π*U^2*α*b*c0*dr^2*m*mr*ρ*sin(α - θ)^2 + 16*π*U^2*α*b*c1*dr^2*m*mr*ρ*sin(α - θ)^2 + 16*π*U^2*α*b*c3*dr^2*m*mr*ρ*sin(α - θ)^2 - 8*π*U*αd*b^2*c0*dr^2*m*mr*ρ*sin(α - θ)^2 + 8*π*U*αd*b^2*c1*dr^2*m*mr*ρ*sin(α - θ)^2 + 8*π*U*αd*b^2*c3*dr^2*m*mr*ρ*sin(α - θ)^2 + 16*π*U^2*b*c1*c2*dr^2*m*mr*ρ*xbd + 16*π*U^2*b*c3*c4*dr^2*m*mr*ρ*xbd + 16*π*U^3*c1*c2*c4*dr^2*m*mr*ρ*xb + 16*π*U^3*c2*c3*c4*dr^2*m*mr*ρ*xb - 8*π*U*b^2*c0*dr*m*mr*ρ*ξd*cos(α) + 8*π*U*b^2*c1*dr*m*mr*ρ*ξd*cos(α) + 8*π*U*b^2*c3*dr*m*mr*ρ*ξd*cos(α) - 16*π*U*b^2*c0*dr^2*m*mr*ρ*ξd*sin(α - θ)^2 + 16*π*U*b^2*c1*dr^2*m*mr*ρ*ξd*sin(α - θ)^2 + 16*π*U*b^2*c3*dr^2*m*mr*ρ*ξd*sin(α - θ)^2 - 8*π*a^2*αd^2*b^3*dr*m*mr*ρ*sin(θ)*cos(α - θ) - 8*π*a*αd^2*b^3*dr^2*m*mr*ρ*sin(α - θ)*cos(α - θ) - 8*π*a^2*b^3*ln*m*mr*ρ*rr*θd*sin(θ) + 16*π*U^2*a*b*c1*c2*m*ρ*xa*xbd + 16*π*U^2*a*b*c3*c4*m*ρ*xa*xbd + 16*π*U^3*a*c1*c2*c4*m*ρ*xa*xb + 16*π*U^3*a*c2*c3*c4*m*ρ*xa*xb - 8*π*a*b^3*dr*m*mr*ρ*rr*θd^2*cos(α - θ) + 4*π*U*αd*b^2*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 8*H*α^3*b*dr*m^2*mr*ωα^2*ra^2*sin(α - θ)*sin(θ) - 16*π*U*a*αd*b^2*c0*dr^2*m*mr*ρ + 16*π*U*a*αd*b^2*c1*dr^2*m*mr*ρ + 16*π*U*a*αd*b^2*c3*dr^2*m*mr*ρ - 8*π*U*a*αd*b^2*dr*m*mr*ρ*cos(α) - 8*π*U*a*αd*b^2*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 8*π*U^2*α*b*c0*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 8*π*U^2*α*b*c1*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 8*π*U^2*α*b*c3*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 4*π*U*αd*b^2*c0*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 4*π*U*αd*b^2*c1*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 4*π*U*αd*b^2*c3*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 8*π*U*b^2*c0*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) + 8*π*U*b^2*c1*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) + 8*π*U*b^2*c3*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) - 16*π*U^2*a*α*b*c0*dr*m*mr*ρ*cos(α) + 16*π*U^2*a*α*b*c1*dr*m*mr*ρ*cos(α) + 16*π*U^2*a*α*b*c3*dr*m*mr*ρ*cos(α) + 16*π*U*a*αd*b^2*c0*dr^2*m*mr*ρ*sin(α - θ)^2 - 16*π*U*a*αd*b^2*c1*dr^2*m*mr*ρ*sin(α - θ)^2 - 16*π*U*a*αd*b^2*c3*dr^2*m*mr*ρ*sin(α - θ)^2 - 16*π*U*a*b^2*c0*dr*m*mr*ρ*ξd*cos(α) + 16*π*U*a*b^2*c1*dr*m*mr*ρ*ξd*cos(α) + 16*π*U*a*b^2*c3*dr*m*mr*ρ*ξd*cos(α) - 8*π*U^2*b*c1*c2*dr*m*mr*ρ*xbd*cos(α) - 8*π*U^2*b*c3*c4*dr*m*mr*ρ*xbd*cos(α) - 8*π*U^3*c1*c2*c4*dr*m*mr*ρ*xb*cos(α) - 8*π*U^3*c2*c3*c4*dr*m*mr*ρ*xb*cos(α) - 16*π*U^2*b*c1*c2*dr^2*m*mr*ρ*xbd*sin(α - θ)^2 - 16*π*U^2*b*c3*c4*dr^2*m*mr*ρ*xbd*sin(α - θ)^2 - 16*π*U^3*c1*c2*c4*dr^2*m*mr*ρ*xb*sin(α - θ)^2 - 16*π*U^3*c2*c3*c4*dr^2*m*mr*ρ*xb*sin(α - θ)^2 + 16*π*U*a^2*αd*b^2*c0*dr*m*mr*ρ*cos(α) - 16*π*U*a^2*αd*b^2*c1*dr*m*mr*ρ*cos(α) - 16*π*U*a^2*αd*b^2*c3*dr*m*mr*ρ*cos(α) - 8*π*a*b^3*dr*ln*m*mr*ρ*rr*θd*sin(α - θ) - 16*π*U^2*a*α*b*c0*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 16*π*U^2*a*α*b*c1*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 16*π*U^2*a*α*b*c3*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 16*π*U*a*b^2*c0*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) + 16*π*U*a*b^2*c1*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) + 16*π*U*a*b^2*c3*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) - 8*π*U^2*b*c1*c2*dr*m*mr*ρ*xbd*sin(α - θ)*sin(θ) - 8*π*U^2*b*c3*c4*dr*m*mr*ρ*xbd*sin(α - θ)*sin(θ) - 8*π*U^3*c1*c2*c4*dr*m*mr*ρ*xb*sin(α - θ)*sin(θ) - 8*π*U^3*c2*c3*c4*dr*m*mr*ρ*xb*sin(α - θ)*sin(θ) - 16*π*U^2*a*b*c1*c2*dr*m*mr*ρ*xbd*cos(α) - 16*π*U^2*a*b*c3*c4*dr*m*mr*ρ*xbd*cos(α) - 16*π*U^3*a*c1*c2*c4*dr*m*mr*ρ*xb*cos(α) - 16*π*U^3*a*c2*c3*c4*dr*m*mr*ρ*xb*cos(α) + 16*π*U*a^2*αd*b^2*c0*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 16*π*U*a^2*αd*b^2*c1*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 16*π*U*a^2*αd*b^2*c3*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 16*π*U^2*a*b*c1*c2*dr*m*mr*ρ*xbd*sin(α - θ)*sin(θ) - 16*π*U^2*a*b*c3*c4*dr*m*mr*ρ*xbd*sin(α - θ)*sin(θ) - 16*π*U^3*a*c1*c2*c4*dr*m*mr*ρ*xb*sin(α - θ)*sin(θ) - 16*π*U^3*a*c2*c3*c4*dr*m*mr*ρ*xb*sin(α - θ)*sin(θ))/(b*(8*m^2*ra^2 - 8*m^2*xa^2 + 8*dr^2*m^2*mr + 8*m^2*mr*ra^2 + b^4*ρ^2*π^2 + 8*dr^2*m^2*mr^2 - 8*dr^2*m^2*mr^2*sin(θ)^2 - 8*dr^2*m^2*mr^2*sin(α - θ)^2 + π*b^2*m*ρ - 8*m^2*mr*ra^2*sin(θ)^2 - 8*dr^2*m^2*mr*sin(α - θ)^2 - 8*dr^2*m^2*mr^2*cos(α)^2 + 8*π*a^2*b^2*m*ρ + 8*π*b^2*m*ra^2*ρ + 16*dr*m^2*mr*xa*cos(α) + π*b^2*m*mr*ρ - π*b^2*m*mr*ρ*sin(θ)^2 + 16*π*a*b^2*m*ρ*xa + 16*dr*m^2*mr*xa*sin(α - θ)*sin(θ) - 16*dr^2*m^2*mr^2*sin(α - θ)*cos(α)*sin(θ) + 8*π*a^2*b^2*m*mr*ρ + 8*π*b^2*dr^2*m*mr*ρ - 8*π*a^2*b^2*m*mr*ρ*sin(θ)^2 - 8*π*b^2*dr^2*m*mr*ρ*sin(α - θ)^2 - 16*π*a*b^2*dr*m*mr*ρ*cos(α) - 16*π*a*b^2*dr*m*mr*ρ*sin(α - θ)*sin(θ)))
αd
(4*(2*b*m^2*ωh^2*xa*ξ - U*αd*b^4*ρ^2*π^2 - 2*α*b*m^2*ωα^2*ra^2 - π*U*αd*b^2*m*ρ + 2*U^2*α*b^3*c0*ρ^2*π^2 - 2*U^2*α*b^3*c1*ρ^2*π^2 - 2*U^2*α*b^3*c3*ρ^2*π^2 - 2*H*α^3*b*m^2*ωα^2*ra^2 + U*αd*b^4*c0*ρ^2*π^2 - U*αd*b^4*c1*ρ^2*π^2 - U*αd*b^4*c3*ρ^2*π^2 + 2*U*b^4*c0*ρ^2*ξd*π^2 - 2*U*b^4*c1*ρ^2*ξd*π^2 - 2*U*b^4*c3*ρ^2*ξd*π^2 - 2*α*b*m^2*mr*ωα^2*ra^2 - 2*U*a*αd*b^4*c0*ρ^2*π^2 + 2*U*a*αd*b^4*c1*ρ^2*π^2 + 2*U*a*αd*b^4*c3*ρ^2*π^2 + 2*α*b*m^2*mr*ωα^2*ra^2*sin(θ)^2 - 2*π*α*b^3*m*ωα^2*ra^2*ρ + 2*π*U*a*αd*b^2*m*ρ + 2*π*U^2*α*b*c0*m*ρ - 2*π*U^2*α*b*c1*m*ρ - 2*π*U^2*α*b*c3*m*ρ + π*U*αd*b^2*c0*m*ρ - π*U*αd*b^2*c1*m*ρ - π*U*αd*b^2*c3*m*ρ - 2*b*dr*m^2*mr*ωh^2*ξ*cos(α) + 2*αd^2*b*dr*m^2*mr*xa*sin(α) - π*U*αd*b^2*m*mr*ρ - 2*αd^2*b*dr^2*m^2*mr^2*cos(α)*sin(α) + 2*π*U*αd*b^2*m*ρ*xa + 2*αd^2*b*dr^2*m^2*mr*sin(α - θ)*cos(α - θ) + 2*π*U*b^2*c0*m*ρ*ξd - 2*π*U*b^2*c1*m*ρ*ξd - 2*π*U*b^2*c3*m*ρ*ξd + 2*b*m^2*mr*rr*θd^2*xa*cos(θ) + 2*U^2*b^3*c1*c2*ρ^2*xbd*π^2 + 2*U^2*b^3*c3*c4*ρ^2*xbd*π^2 + 2*b*dr*m^2*mr*rr*θd^2*cos(α - θ) - 2*H*α^3*b*m^2*mr*ωα^2*ra^2 + 2*αd^2*b*dr^2*m^2*mr^2*sin(α - θ)*cos(α - θ) - 2*π*a*b^3*m*ωh^2*ρ*ξ + 2*b*dr*m^2*mr^2*rr*θd^2*cos(α - θ) + 2*b*dr*ln*m^2*mr^2*rr*θd*sin(α - θ) - 4*π*U*a^2*αd*b^2*c0*m*ρ + 4*π*U*a^2*αd*b^2*c1*m*ρ + 4*π*U*a^2*αd*b^2*c3*m*ρ - 2*b*dr*m^2*mr^2*rr*θd^2*sin(θ)^2*cos(α - θ) - 2*αd^2*b*dr*m^2*mr*xa*sin(θ)*cos(α - θ) + 2*αd^2*b*dr^2*m^2*mr^2*cos(α)*sin(θ)*cos(α - θ) - 2*b*dr*m^2*mr*ωh^2*ξ*sin(α - θ)*sin(θ) - 2*b*dr*m^2*mr^2*rr*θd^2*cos(α)*cos(θ) - 2*αd^2*b*dr^2*m^2*mr^2*sin(α - θ)*sin(α)*sin(θ) - 2*b*ln*m^2*mr*rr*θd*xa*sin(θ) + 2*b*dr*ln*m^2*mr*rr*θd*sin(α - θ) + π*U*αd*b^2*m*mr*ρ*sin(θ)^2 + 4*π*U^2*a*α*b*c0*m*ρ - 4*π*U^2*a*α*b*c1*m*ρ - 4*π*U^2*a*α*b*c3*m*ρ + 2*π*U*a*αd*b^2*m*mr*ρ + 2*π*U^2*α*b*c0*m*mr*ρ - 2*π*U^2*α*b*c1*m*mr*ρ - 2*π*U^2*α*b*c3*m*mr*ρ + π*U*αd*b^2*c0*m*mr*ρ - π*U*αd*b^2*c1*m*mr*ρ - π*U*αd*b^2*c3*m*mr*ρ + 4*π*U*a*b^2*c0*m*ρ*ξd - 4*π*U*a*b^2*c1*m*ρ*ξd - 4*π*U*a*b^2*c3*m*ρ*ξd + 4*π*U^2*α*b*c0*m*ρ*xa - 4*π*U^2*α*b*c1*m*ρ*xa - 4*π*U^2*α*b*c3*m*ρ*xa + 2*π*U*αd*b^2*c0*m*ρ*xa - 2*π*U*αd*b^2*c1*m*ρ*xa - 2*π*U*αd*b^2*c3*m*ρ*xa + 2*π*U^2*b*c1*c2*m*ρ*xbd + 2*π*U^2*b*c3*c4*m*ρ*xbd + 2*π*U^3*c1*c2*c4*m*ρ*xb + 2*π*U^3*c2*c3*c4*m*ρ*xb + 2*H*α^3*b*m^2*mr*ωα^2*ra^2*sin(θ)^2 + 2*π*U*b^2*c0*m*mr*ρ*ξd - 2*π*U*b^2*c1*m*mr*ρ*ξd - 2*π*U*b^2*c3*m*mr*ρ*ξd + 2*U^3*b^2*c1*c2*c4*ρ^2*xb*π^2 + 2*U^3*b^2*c2*c3*c4*ρ^2*xb*π^2 + 4*π*U*b^2*c0*m*ρ*xa*ξd - 4*π*U*b^2*c1*m*ρ*xa*ξd - 4*π*U*b^2*c3*m*ρ*xa*ξd - 2*π*H*α^3*b^3*m*ωα^2*ra^2*ρ - 2*π*a*αd^2*b^3*dr*m*mr*ρ*sin(α) - 2*π*U*a*αd*b^2*m*mr*ρ*sin(θ)^2 - 2*π*U^2*α*b*c0*m*mr*ρ*sin(θ)^2 + 2*π*U^2*α*b*c1*m*mr*ρ*sin(θ)^2 + 2*π*U^2*α*b*c3*m*mr*ρ*sin(θ)^2 - π*U*αd*b^2*c0*m*mr*ρ*sin(θ)^2 + π*U*αd*b^2*c1*m*mr*ρ*sin(θ)^2 + π*U*αd*b^2*c3*m*mr*ρ*sin(θ)^2 - 2*π*U*b^2*c0*m*mr*ρ*ξd*sin(θ)^2 + 2*π*U*b^2*c1*m*mr*ρ*ξd*sin(θ)^2 + 2*π*U*b^2*c3*m*mr*ρ*ξd*sin(θ)^2 + 4*π*U^2*a*α*b*c0*m*mr*ρ - 4*π*U^2*a*α*b*c1*m*mr*ρ - 4*π*U^2*a*α*b*c3*m*mr*ρ - 2*π*a*b^3*m*mr*ρ*rr*θd^2*cos(θ) + 2*π*αd^2*b^3*dr^2*m*mr*ρ*sin(α - θ)*cos(α - θ) - 4*π*U*a*αd*b^2*c0*m*ρ*xa + 4*π*U*a*αd*b^2*c1*m*ρ*xa + 4*π*U*a*αd*b^2*c3*m*ρ*xa + 4*π*U^2*a*b*c1*c2*m*ρ*xbd + 4*π*U^2*a*b*c3*c4*m*ρ*xbd + 4*π*U^3*a*c1*c2*c4*m*ρ*xb + 4*π*U^3*a*c2*c3*c4*m*ρ*xb + 4*π*U*a*b^2*c0*m*mr*ρ*ξd - 4*π*U*a*b^2*c1*m*mr*ρ*ξd - 4*π*U*a*b^2*c3*m*mr*ρ*ξd + 2*π*U^2*b*c1*c2*m*mr*ρ*xbd + 2*π*U^2*b*c3*c4*m*mr*ρ*xbd + 2*π*U^3*c1*c2*c4*m*mr*ρ*xb + 2*π*U^3*c2*c3*c4*m*mr*ρ*xb + 4*π*U^2*b*c1*c2*m*ρ*xa*xbd + 4*π*U^2*b*c3*c4*m*ρ*xa*xbd + 4*π*U^3*c1*c2*c4*m*ρ*xa*xb + 4*π*U^3*c2*c3*c4*m*ρ*xa*xb + 2*π*b^3*dr*m*mr*ρ*rr*θd^2*cos(α - θ) - 4*π*U*a^2*αd*b^2*c0*m*mr*ρ + 4*π*U*a^2*αd*b^2*c1*m*mr*ρ + 4*π*U*a^2*αd*b^2*c3*m*mr*ρ - 2*π*U*αd*b^2*dr*m*mr*ρ*cos(α) + 2*b*dr*ln*m^2*mr^2*rr*θd*cos(α)*sin(θ) - 2*b*dr*m^2*mr^2*rr*θd^2*sin(α - θ)*cos(θ)*sin(θ) - 4*π*U^2*α*b*c0*dr*m*mr*ρ*cos(α) + 4*π*U^2*α*b*c1*dr*m*mr*ρ*cos(α) + 4*π*U^2*α*b*c3*dr*m*mr*ρ*cos(α) - 2*π*U*αd*b^2*c0*dr*m*mr*ρ*cos(α) + 2*π*U*αd*b^2*c1*dr*m*mr*ρ*cos(α) + 2*π*U*αd*b^2*c3*dr*m*mr*ρ*cos(α) - 4*π*U*b^2*c0*dr*m*mr*ρ*ξd*cos(α) + 4*π*U*b^2*c1*dr*m*mr*ρ*ξd*cos(α) + 4*π*U*b^2*c3*dr*m*mr*ρ*ξd*cos(α) + 2*π*a*αd^2*b^3*dr*m*mr*ρ*sin(θ)*cos(α - θ) + 2*π*a*b^3*ln*m*mr*ρ*rr*θd*sin(θ) + 2*π*b^3*dr*ln*m*mr*ρ*rr*θd*sin(α - θ) - 4*π*U^2*a*α*b*c0*m*mr*ρ*sin(θ)^2 + 4*π*U^2*a*α*b*c1*m*mr*ρ*sin(θ)^2 + 4*π*U^2*a*α*b*c3*m*mr*ρ*sin(θ)^2 - 4*π*U*a*b^2*c0*m*mr*ρ*ξd*sin(θ)^2 + 4*π*U*a*b^2*c1*m*mr*ρ*ξd*sin(θ)^2 + 4*π*U*a*b^2*c3*m*mr*ρ*ξd*sin(θ)^2 - 2*π*U^2*b*c1*c2*m*mr*ρ*xbd*sin(θ)^2 - 2*π*U^2*b*c3*c4*m*mr*ρ*xbd*sin(θ)^2 - 2*π*U^3*c1*c2*c4*m*mr*ρ*xb*sin(θ)^2 - 2*π*U^3*c2*c3*c4*m*mr*ρ*xb*sin(θ)^2 + 4*π*U^2*a*b*c1*c2*m*mr*ρ*xbd + 4*π*U^2*a*b*c3*c4*m*mr*ρ*xbd + 4*π*U^3*a*c1*c2*c4*m*mr*ρ*xb + 4*π*U^3*a*c2*c3*c4*m*mr*ρ*xb - 2*π*U*αd*b^2*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 4*π*U*a^2*αd*b^2*c0*m*mr*ρ*sin(θ)^2 - 4*π*U*a^2*αd*b^2*c1*m*mr*ρ*sin(θ)^2 - 4*π*U*a^2*αd*b^2*c3*m*mr*ρ*sin(θ)^2 - 4*π*U^2*α*b*c0*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 4*π*U^2*α*b*c1*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 4*π*U^2*α*b*c3*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 2*π*U*αd*b^2*c0*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 2*π*U*αd*b^2*c1*dr*m*mr*ρ*sin(α - θ)*sin(θ) + 2*π*U*αd*b^2*c3*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 4*π*U*b^2*c0*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) + 4*π*U*b^2*c1*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) + 4*π*U*b^2*c3*dr*m*mr*ρ*ξd*sin(α - θ)*sin(θ) + 4*π*U*a*αd*b^2*c0*dr*m*mr*ρ*cos(α) - 4*π*U*a*αd*b^2*c1*dr*m*mr*ρ*cos(α) - 4*π*U*a*αd*b^2*c3*dr*m*mr*ρ*cos(α) - 4*π*U^2*b*c1*c2*dr*m*mr*ρ*xbd*cos(α) - 4*π*U^2*b*c3*c4*dr*m*mr*ρ*xbd*cos(α) - 4*π*U^3*c1*c2*c4*dr*m*mr*ρ*xb*cos(α) - 4*π*U^3*c2*c3*c4*dr*m*mr*ρ*xb*cos(α) - 4*π*U^2*a*b*c1*c2*m*mr*ρ*xbd*sin(θ)^2 - 4*π*U^2*a*b*c3*c4*m*mr*ρ*xbd*sin(θ)^2 - 4*π*U^3*a*c1*c2*c4*m*mr*ρ*xb*sin(θ)^2 - 4*π*U^3*a*c2*c3*c4*m*mr*ρ*xb*sin(θ)^2 + 4*π*U*a*αd*b^2*c0*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 4*π*U*a*αd*b^2*c1*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 4*π*U*a*αd*b^2*c3*dr*m*mr*ρ*sin(α - θ)*sin(θ) - 4*π*U^2*b*c1*c2*dr*m*mr*ρ*xbd*sin(α - θ)*sin(θ) - 4*π*U^2*b*c3*c4*dr*m*mr*ρ*xbd*sin(α - θ)*sin(θ) - 4*π*U^3*c1*c2*c4*dr*m*mr*ρ*xb*sin(α - θ)*sin(θ) - 4*π*U^3*c2*c3*c4*dr*m*mr*ρ*xb*sin(α - θ)*sin(θ)))/(b*(8*m^2*ra^2 - 8*m^2*xa^2 + 8*dr^2*m^2*mr + 8*m^2*mr*ra^2 + b^4*ρ^2*π^2 + 8*dr^2*m^2*mr^2 - 8*dr^2*m^2*mr^2*sin(θ)^2 - 8*dr^2*m^2*mr^2*sin(α - θ)^2 + π*b^2*m*ρ - 8*m^2*mr*ra^2*sin(θ)^2 - 8*dr^2*m^2*mr*sin(α - θ)^2 - 8*dr^2*m^2*mr^2*cos(α)^2 + 8*π*a^2*b^2*m*ρ + 8*π*b^2*m*ra^2*ρ + 16*dr*m^2*mr*xa*cos(α) + π*b^2*m*mr*ρ - π*b^2*m*mr*ρ*sin(θ)^2 + 16*π*a*b^2*m*ρ*xa + 16*dr*m^2*mr*xa*sin(α - θ)*sin(θ) - 16*dr^2*m^2*mr^2*sin(α - θ)*cos(α)*sin(θ) + 8*π*a^2*b^2*m*mr*ρ + 8*π*b^2*dr^2*m*mr*ρ - 8*π*a^2*b^2*m*mr*ρ*sin(θ)^2 - 8*π*b^2*dr^2*m*mr*ρ*sin(α - θ)^2 - 16*π*a*b^2*dr*m*mr*ρ*cos(α) - 16*π*a*b^2*dr*m*mr*ρ*sin(α - θ)*sin(θ)))
θd
-(- U*αd*b^4*ρ^2*π^2*sin(θ) + 8*αd^2*b*dr^3*m^2*mr*cos(α - θ) + 8*αd^2*b*dr*m^2*ra^2*cos(α - θ) - 8*αd^2*b*dr*m^2*xa^2*cos(α - θ) - 8*b*m^2*ωh^2*ra^2*ξ*sin(θ) + 8*b*ln*m^2*ra^2*rr*θd - 8*b*ln*m^2*rr*θd*xa^2 + αd^2*b^5*dr*ρ^2*π^2*cos(α - θ) + 8*αd^2*b*dr^3*m^2*mr^2*cos(α - θ) + b^5*ln*ρ^2*rr*θd*π^2 - 2*U^2*α*b^3*c0*ρ^2*π^2*sin(θ) + 2*U^2*α*b^3*c1*ρ^2*π^2*sin(θ) + 2*U^2*α*b^3*c3*ρ^2*π^2*sin(θ) + 8*b*dr^2*ln*m^2*mr*rr*θd + 8*b*ln*m^2*mr*ra^2*rr*θd + αd^2*b^3*dr*m*ρ*π*cos(α - θ) + 8*b*dr^2*ln*m^2*mr^2*rr*θd - b^3*m*ωh^2*ρ*ξ*π*sin(θ) - 8*αd^2*b*dr^3*m^2*mr^2*sin(α)*sin(θ) - 4*U*a*αd*b^4*ρ^2*π^2*sin(θ) - U*αd*b^4*c0*ρ^2*π^2*sin(θ) + U*αd*b^4*c1*ρ^2*π^2*sin(θ) + U*αd*b^4*c3*ρ^2*π^2*sin(θ) + b^3*ln*m*ρ*rr*θd*π + 8*b*dr*m^2*ωh^2*xa*ξ*sin(α - θ) - 2*U*b^4*c0*ρ^2*ξd*π^2*sin(θ) + 2*U*b^4*c1*ρ^2*ξd*π^2*sin(θ) + 2*U*b^4*c3*ρ^2*ξd*π^2*sin(θ) - 4*U*αd*b^4*dr*ρ^2*π^2*sin(α - θ) + 8*αd^2*b*dr*m^2*mr*ra^2*cos(α - θ) - 8*b*dr^2*m^2*mr*ωh^2*ξ*sin(θ) + 8*α*b*m^2*ωα^2*ra^2*xa*sin(θ) - 8*αd^2*b*dr^3*m^2*mr^2*cos(α)^2*cos(α - θ) - 8*α*b*dr*m^2*ωα^2*ra^2*sin(α - θ) + 6*U*a*αd*b^4*c0*ρ^2*π^2*sin(θ) - 6*U*a*αd*b^4*c1*ρ^2*π^2*sin(θ) - 6*U*a*αd*b^4*c3*ρ^2*π^2*sin(θ) + 8*b*dr^2*m^2*mr*rr*θd^2*sin(α - θ)*cos(α - θ) + b^3*ln*m*mr*ρ*rr*θd*π + 8*U*a*b^4*c0*ρ^2*ξd*π^2*sin(θ) - 8*U*a*b^4*c1*ρ^2*ξd*π^2*sin(θ) - 8*U*a*b^4*c3*ρ^2*ξd*π^2*sin(θ) - 8*b*dr^2*m^2*mr^2*rr*θd^2*cos(θ)*sin(θ) + 4*U*αd*b^4*c0*dr*ρ^2*π^2*sin(α - θ) - 4*U*αd*b^4*c1*dr*ρ^2*π^2*sin(α - θ) - 4*U*αd*b^4*c3*dr*ρ^2*π^2*sin(α - θ) + 8*a^2*αd^2*b^3*dr*m*ρ*π*cos(α - θ) + 8*αd^2*b^3*dr^3*m*mr*ρ*π*cos(α - θ) + 8*U*b^4*c0*dr*ρ^2*ξd*π^2*sin(α - θ) - 8*U*b^4*c1*dr*ρ^2*ξd*π^2*sin(α - θ) - 8*U*b^4*c3*dr*ρ^2*ξd*π^2*sin(α - θ) + 8*αd^2*b^3*dr*m*ra^2*ρ*π*cos(α - θ) - 8*a^2*b^3*m*ωh^2*ρ*ξ*π*sin(θ) + 8*U^2*a*α*b^3*c0*ρ^2*π^2*sin(θ) - 8*U^2*a*α*b^3*c1*ρ^2*π^2*sin(θ) - 8*U^2*a*α*b^3*c3*ρ^2*π^2*sin(θ) - 8*U*a^2*αd*b^4*c0*ρ^2*π^2*sin(θ) + 8*U*a^2*αd*b^4*c1*ρ^2*π^2*sin(θ) + 8*U*a^2*αd*b^4*c3*ρ^2*π^2*sin(θ) + 8*a^2*b^3*ln*m*ρ*rr*θd*π + 8*b*dr^2*m^2*mr^2*rr*θd^2*sin(α - θ)*cos(α - θ) + 4*U*αd*b^2*m*ρ*xa*π*sin(θ) - 8*α*b*dr*m^2*mr*ωα^2*ra^2*sin(α - θ) - 4*U*αd*b^2*dr*m*ρ*π*sin(α - θ) + 8*b^3*ln*m*ra^2*ρ*rr*θd*π - 2*U^2*b^3*c1*c2*ρ^2*xbd*π^2*sin(θ) - 2*U^2*b^3*c3*c4*ρ^2*xbd*π^2*sin(θ) - 8*αd^2*b*dr^3*m^2*mr^2*sin(α - θ)*cos(α)*sin(α) + 8*U^2*α*b^3*c0*dr*ρ^2*π^2*sin(α - θ) - 8*U^2*α*b^3*c1*dr*ρ^2*π^2*sin(α - θ) - 8*U^2*α*b^3*c3*dr*ρ^2*π^2*sin(α - θ) - 8*b*dr^2*ln*m^2*mr^2*rr*θd*cos(α)^2 + 8*H*α^3*b*m^2*ωα^2*ra^2*xa*sin(θ) - 8*H*α^3*b*dr*m^2*ωα^2*ra^2*sin(α - θ) - 8*αd^2*b*dr*m^2*mr*ra^2*sin(α)*sin(θ) + 16*αd^2*b*dr^2*m^2*mr*xa*cos(α)*cos(α - θ) - 8*U*αd*b^2*m*ra^2*ρ*π*sin(θ) - 8*b*m^2*mr*ra^2*rr*θd^2*cos(θ)*sin(θ) - 8*b*dr^2*m^2*mr*ωh^2*ξ*sin(α - θ)*cos(α) + 8*αd^2*b*dr^2*m^2*mr*xa*sin(α - θ)*sin(α) + αd^2*b^3*dr*m*mr*ρ*π*cos(α - θ) + 8*U^2*b^3*c1*c2*dr*ρ^2*xbd*π^2*sin(α - θ) + 8*U^2*b^3*c3*c4*dr*ρ^2*xbd*π^2*sin(α - θ) - 8*b*dr^2*m^2*mr^2*rr*θd^2*sin(α - θ)*cos(α)*cos(θ) + 8*b*dr^2*m^2*mr^2*rr*θd^2*cos(α)*sin(θ)*cos(α - θ) - 8*α*b*dr*m^2*mr*ωα^2*ra^2*cos(α)*sin(θ) - 8*H*α^3*b*dr*m^2*mr*ωα^2*ra^2*sin(α - θ) - 8*U*αd*b^2*dr^2*m*mr*ρ*π*sin(θ) - 16*U^2*α*b*c0*m*ra^2*ρ*π*sin(θ) + 16*U^2*α*b*c1*m*ra^2*ρ*π*sin(θ) + 16*U^2*α*b*c3*m*ra^2*ρ*π*sin(θ) - 8*U*αd*b^2*c0*m*ra^2*ρ*π*sin(θ) + 8*U*αd*b^2*c1*m*ra^2*ρ*π*sin(θ) + 8*U*αd*b^2*c3*m*ra^2*ρ*π*sin(θ) - 16*U*b^2*c0*m*ra^2*ρ*ξd*π*sin(θ) + 16*U*b^2*c1*m*ra^2*ρ*ξd*π*sin(θ) + 16*U*b^2*c3*m*ra^2*ρ*ξd*π*sin(θ) + 16*a*αd^2*b^3*dr*m*ρ*xa*π*cos(α - θ) - 8*a*b^3*dr*m*ωh^2*ρ*ξ*π*sin(α - θ) + 16*a*b^3*ln*m*ρ*rr*θd*xa*π - 8*U*a*αd*b^4*c0*dr*ρ^2*π^2*sin(α - θ) + 8*U*a*αd*b^4*c1*dr*ρ^2*π^2*sin(α - θ) + 8*U*a*αd*b^4*c3*dr*ρ^2*π^2*sin(α - θ) - 8*a*α*b^3*m*ωα^2*ra^2*ρ*π*sin(θ) + 8*a^2*αd^2*b^3*dr*m*mr*ρ*π*cos(α - θ) - αd^2*b^3*dr*m*mr*ρ*π*sin(α)*sin(θ) - 8*α*b^3*dr*m*ωα^2*ra^2*ρ*π*sin(α - θ) - b^3*m*mr*ρ*rr*θd^2*π*cos(θ)*sin(θ) - 8*U*a*αd*b^2*m*ρ*xa*π*sin(θ) - 8*U^2*α*b*c0*m*ρ*xa*π*sin(θ) + 8*U^2*α*b*c1*m*ρ*xa*π*sin(θ) + 8*U^2*α*b*c3*m*ρ*xa*π*sin(θ) - 4*U*αd*b^2*c0*m*ρ*xa*π*sin(θ) + 4*U*αd*b^2*c1*m*ρ*xa*π*sin(θ) + 4*U*αd*b^2*c3*m*ρ*xa*π*sin(θ) + 8*a^2*b^3*ln*m*mr*ρ*rr*θd*π + 8*U*a*αd*b^2*dr*m*ρ*π*sin(α - θ) + 8*b^3*dr^2*ln*m*mr*ρ*rr*θd*π + 8*U^2*α*b*c0*dr*m*ρ*π*sin(α - θ) - 8*U^2*α*b*c1*dr*m*ρ*π*sin(α - θ) - 8*U^2*α*b*c3*dr*m*ρ*π*sin(α - θ) + 4*U*αd*b^2*c0*dr*m*ρ*π*sin(α - θ) - 4*U*αd*b^2*c1*dr*m*ρ*π*sin(α - θ) - 4*U*αd*b^2*c3*dr*m*ρ*π*sin(α - θ) + 8*U^2*a*b^3*c1*c2*ρ^2*xbd*π^2*sin(θ) + 8*U^2*a*b^3*c3*c4*ρ^2*xbd*π^2*sin(θ) - 2*U^3*b^2*c1*c2*c4*ρ^2*xb*π^2*sin(θ) - 2*U^3*b^2*c2*c3*c4*ρ^2*xb*π^2*sin(θ) - 4*U*αd*b^2*dr*m*mr*ρ*π*sin(α - θ) - 8*U*b^2*c0*m*ρ*xa*ξd*π*sin(θ) + 8*U*b^2*c1*m*ρ*xa*ξd*π*sin(θ) + 8*U*b^2*c3*m*ρ*xa*ξd*π*sin(θ) + 8*U*αd*b^2*dr*m*ρ*xa*π*sin(α - θ) + 8*U*b^2*c0*dr*m*ρ*ξd*π*sin(α - θ) - 8*U*b^2*c1*dr*m*ρ*ξd*π*sin(α - θ) - 8*U*b^2*c3*dr*m*ρ*ξd*π*sin(α - θ) + 16*b*dr*ln*m^2*mr*rr*θd*xa*cos(α) + 8*b*dr*m^2*mr*rr*θd^2*xa*sin(α - θ)*cos(θ) - 8*b*dr*m^2*mr*rr*θd^2*xa*sin(θ)*cos(α - θ) - 8*U*αd*b^2*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) - 16*U^2*a*α*b*c0*m*ρ*xa*π*sin(θ) + 16*U^2*a*α*b*c1*m*ρ*xa*π*sin(θ) + 16*U^2*a*α*b*c3*m*ρ*xa*π*sin(θ) + 16*U^2*a*α*b*c0*dr*m*ρ*π*sin(α - θ) - 16*U^2*a*α*b*c1*dr*m*ρ*π*sin(α - θ) - 16*U^2*a*α*b*c3*dr*m*ρ*π*sin(α - θ) + 8*U^3*a*b^2*c1*c2*c4*ρ^2*xb*π^2*sin(θ) + 8*U^3*a*b^2*c2*c3*c4*ρ^2*xb*π^2*sin(θ) + 8*U*a*αd*b^2*dr*m*mr*ρ*π*sin(α - θ) - 16*U*a*b^2*c0*m*ρ*xa*ξd*π*sin(θ) + 16*U*a*b^2*c1*m*ρ*xa*ξd*π*sin(θ) + 16*U*a*b^2*c3*m*ρ*xa*ξd*π*sin(θ) - 8*H*a*α^3*b^3*m*ωα^2*ra^2*ρ*π*sin(θ) + 8*U^2*α*b*c0*dr*m*mr*ρ*π*sin(α - θ) - 8*U^2*α*b*c1*dr*m*mr*ρ*π*sin(α - θ) - 8*U^2*α*b*c3*dr*m*mr*ρ*π*sin(α - θ) + 4*U*αd*b^2*c0*dr*m*mr*ρ*π*sin(α - θ) - 4*U*αd*b^2*c1*dr*m*mr*ρ*π*sin(α - θ) - 4*U*αd*b^2*c3*dr*m*mr*ρ*π*sin(α - θ) - 8*U^2*b*c1*c2*m*ρ*xa*xbd*π*sin(θ) - 8*U^2*b*c3*c4*m*ρ*xa*xbd*π*sin(θ) - 8*U^3*c1*c2*c4*m*ρ*xa*xb*π*sin(θ) - 8*U^3*c2*c3*c4*m*ρ*xa*xb*π*sin(θ) + 16*U*a*b^2*c0*dr*m*ρ*ξd*π*sin(α - θ) - 16*U*a*b^2*c1*dr*m*ρ*ξd*π*sin(α - θ) - 16*U*a*b^2*c3*dr*m*ρ*ξd*π*sin(α - θ) + 16*U^2*α*b*c0*dr*m*ρ*xa*π*sin(α - θ) - 16*U^2*α*b*c1*dr*m*ρ*xa*π*sin(α - θ) - 16*U^2*α*b*c3*dr*m*ρ*xa*π*sin(α - θ) + 8*U*αd*b^2*c0*dr*m*ρ*xa*π*sin(α - θ) - 8*U*αd*b^2*c1*dr*m*ρ*xa*π*sin(α - θ) - 8*U*αd*b^2*c3*dr*m*ρ*xa*π*sin(α - θ) + 8*U^2*b*c1*c2*dr*m*ρ*xbd*π*sin(α - θ) + 8*U^2*b*c3*c4*dr*m*ρ*xbd*π*sin(α - θ) + 8*U^3*c1*c2*c4*dr*m*ρ*xb*π*sin(α - θ) + 8*U^3*c2*c3*c4*dr*m*ρ*xb*π*sin(α - θ) + 8*U*b^2*c0*dr*m*mr*ρ*ξd*π*sin(α - θ) - 8*U*b^2*c1*dr*m*mr*ρ*ξd*π*sin(α - θ) - 8*U*b^2*c3*dr*m*mr*ρ*ξd*π*sin(α - θ) - 8*a^2*αd^2*b^3*dr*m*mr*ρ*π*sin(α)*sin(θ) - 16*a*αd^2*b^3*dr^2*m*mr*ρ*π*cos(α)*cos(α - θ) + 8*U^3*b^2*c1*c2*c4*dr*ρ^2*xb*π^2*sin(α - θ) + 8*U^3*b^2*c2*c3*c4*dr*ρ^2*xb*π^2*sin(α - θ) + 16*U*b^2*c0*dr*m*ρ*xa*ξd*π*sin(α - θ) - 16*U*b^2*c1*dr*m*ρ*xa*ξd*π*sin(α - θ) - 16*U*b^2*c3*dr*m*ρ*xa*ξd*π*sin(α - θ) - 8*H*α^3*b^3*dr*m*ωα^2*ra^2*ρ*π*sin(α - θ) - 8*a*αd^2*b^3*dr^2*m*mr*ρ*π*sin(α - θ)*sin(α) - 16*U^2*α*b*c0*dr^2*m*mr*ρ*π*sin(θ) + 16*U^2*α*b*c1*dr^2*m*mr*ρ*π*sin(θ) + 16*U^2*α*b*c3*dr^2*m*mr*ρ*π*sin(θ) - 8*U*αd*b^2*c0*dr^2*m*mr*ρ*π*sin(θ) + 8*U*αd*b^2*c1*dr^2*m*mr*ρ*π*sin(θ) + 8*U*αd*b^2*c3*dr^2*m*mr*ρ*π*sin(θ) - 8*a^2*b^3*m*mr*ρ*rr*θd^2*π*cos(θ)*sin(θ) + 16*U*a*αd*b^2*c0*m*ra^2*ρ*π*sin(θ) - 16*U*a*αd*b^2*c1*m*ra^2*ρ*π*sin(θ) - 16*U*a*αd*b^2*c3*m*ra^2*ρ*π*sin(θ) - 4*U*αd*b^2*dr*m*mr*ρ*π*cos(α)*sin(θ) + 16*U*a^2*αd*b^2*c0*m*ρ*xa*π*sin(θ) - 16*U*a^2*αd*b^2*c1*m*ρ*xa*π*sin(θ) - 16*U*a^2*αd*b^2*c3*m*ρ*xa*π*sin(θ) - 16*U*a^2*αd*b^2*c0*dr*m*ρ*π*sin(α - θ) + 16*U*a^2*αd*b^2*c1*dr*m*ρ*π*sin(α - θ) + 16*U*a^2*αd*b^2*c3*dr*m*ρ*π*sin(α - θ) - 16*U*b^2*c0*dr^2*m*mr*ρ*ξd*π*sin(θ) + 16*U*b^2*c1*dr^2*m*mr*ρ*ξd*π*sin(θ) + 16*U*b^2*c3*dr^2*m*mr*ρ*ξd*π*sin(θ) - 16*U^2*b*c1*c2*m*ra^2*ρ*xbd*π*sin(θ) - 16*U^2*b*c3*c4*m*ra^2*ρ*xbd*π*sin(θ) - 16*U^3*c1*c2*c4*m*ra^2*ρ*xb*π*sin(θ) - 16*U^3*c2*c3*c4*m*ra^2*ρ*xb*π*sin(θ) - 8*H*α^3*b*dr*m^2*mr*ωα^2*ra^2*cos(α)*sin(θ) + 8*b^3*dr^2*m*mr*ρ*rr*θd^2*π*sin(α - θ)*cos(α - θ) - 16*U*a^2*αd*b^2*c0*dr*m*mr*ρ*π*sin(α - θ) + 16*U*a^2*αd*b^2*c1*dr*m*mr*ρ*π*sin(α - θ) + 16*U*a^2*αd*b^2*c3*dr*m*mr*ρ*π*sin(α - θ) - 16*U^2*α*b*c0*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) + 16*U^2*α*b*c1*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) + 16*U^2*α*b*c3*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) - 8*U*αd*b^2*c0*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) + 8*U*αd*b^2*c1*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) + 8*U*αd*b^2*c3*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) - 16*U*b^2*c0*dr^2*m*mr*ρ*ξd*π*sin(α - θ)*cos(α) + 16*U*b^2*c1*dr^2*m*mr*ρ*ξd*π*sin(α - θ)*cos(α) + 16*U*b^2*c3*dr^2*m*mr*ρ*ξd*π*sin(α - θ)*cos(α) - 16*a*b^3*dr*ln*m*mr*ρ*rr*θd*π*cos(α) + 16*U^2*a*α*b*c0*dr*m*mr*ρ*π*sin(α - θ) - 16*U^2*a*α*b*c1*dr*m*mr*ρ*π*sin(α - θ) - 16*U^2*a*α*b*c3*dr*m*mr*ρ*π*sin(α - θ) - 16*U^2*a*b*c1*c2*m*ρ*xa*xbd*π*sin(θ) - 16*U^2*a*b*c3*c4*m*ρ*xa*xbd*π*sin(θ) - 16*U^3*a*c1*c2*c4*m*ρ*xa*xb*π*sin(θ) - 16*U^3*a*c2*c3*c4*m*ρ*xa*xb*π*sin(θ) - 8*a*b^3*dr*m*mr*ρ*rr*θd^2*π*sin(α - θ)*cos(θ) + 8*a*b^3*dr*m*mr*ρ*rr*θd^2*π*sin(θ)*cos(α - θ) - 16*U*a*αd*b^2*c0*dr*m*ρ*xa*π*sin(α - θ) + 16*U*a*αd*b^2*c1*dr*m*ρ*xa*π*sin(α - θ) + 16*U*a*αd*b^2*c3*dr*m*ρ*xa*π*sin(α - θ) + 16*U^2*a*b*c1*c2*dr*m*ρ*xbd*π*sin(α - θ) + 16*U^2*a*b*c3*c4*dr*m*ρ*xbd*π*sin(α - θ) + 16*U^3*a*c1*c2*c4*dr*m*ρ*xb*π*sin(α - θ) + 16*U^3*a*c2*c3*c4*dr*m*ρ*xb*π*sin(α - θ) + 16*U*a*b^2*c0*dr*m*mr*ρ*ξd*π*sin(α - θ) - 16*U*a*b^2*c1*dr*m*mr*ρ*ξd*π*sin(α - θ) - 16*U*a*b^2*c3*dr*m*mr*ρ*ξd*π*sin(α - θ) + 8*U^2*b*c1*c2*dr*m*mr*ρ*xbd*π*sin(α - θ) + 8*U^2*b*c3*c4*dr*m*mr*ρ*xbd*π*sin(α - θ) + 8*U^3*c1*c2*c4*dr*m*mr*ρ*xb*π*sin(α - θ) + 8*U^3*c2*c3*c4*dr*m*mr*ρ*xb*π*sin(α - θ) + 16*U^2*b*c1*c2*dr*m*ρ*xa*xbd*π*sin(α - θ) + 16*U^2*b*c3*c4*dr*m*ρ*xa*xbd*π*sin(α - θ) + 16*U^3*c1*c2*c4*dr*m*ρ*xa*xb*π*sin(α - θ) + 16*U^3*c2*c3*c4*dr*m*ρ*xa*xb*π*sin(α - θ) + 16*U*a*αd*b^2*c0*dr^2*m*mr*ρ*π*sin(θ) - 16*U*a*αd*b^2*c1*dr^2*m*mr*ρ*π*sin(θ) - 16*U*a*αd*b^2*c3*dr^2*m*mr*ρ*π*sin(θ) + 8*U*a*αd*b^2*dr*m*mr*ρ*π*cos(α)*sin(θ) + 8*U^2*α*b*c0*dr*m*mr*ρ*π*cos(α)*sin(θ) - 8*U^2*α*b*c1*dr*m*mr*ρ*π*cos(α)*sin(θ) - 8*U^2*α*b*c3*dr*m*mr*ρ*π*cos(α)*sin(θ) + 4*U*αd*b^2*c0*dr*m*mr*ρ*π*cos(α)*sin(θ) - 4*U*αd*b^2*c1*dr*m*mr*ρ*π*cos(α)*sin(θ) - 4*U*αd*b^2*c3*dr*m*mr*ρ*π*cos(α)*sin(θ) - 16*U^2*b*c1*c2*dr^2*m*mr*ρ*xbd*π*sin(θ) - 16*U^2*b*c3*c4*dr^2*m*mr*ρ*xbd*π*sin(θ) - 16*U^3*c1*c2*c4*dr^2*m*mr*ρ*xb*π*sin(θ) - 16*U^3*c2*c3*c4*dr^2*m*mr*ρ*xb*π*sin(θ) + 8*U*b^2*c0*dr*m*mr*ρ*ξd*π*cos(α)*sin(θ) - 8*U*b^2*c1*dr*m*mr*ρ*ξd*π*cos(α)*sin(θ) - 8*U*b^2*c3*dr*m*mr*ρ*ξd*π*cos(α)*sin(θ) + 16*U^2*a*b*c1*c2*dr*m*mr*ρ*xbd*π*sin(α - θ) + 16*U^2*a*b*c3*c4*dr*m*mr*ρ*xbd*π*sin(α - θ) + 16*U^3*a*c1*c2*c4*dr*m*mr*ρ*xb*π*sin(α - θ) + 16*U^3*a*c2*c3*c4*dr*m*mr*ρ*xb*π*sin(α - θ) + 16*U^2*a*α*b*c0*dr*m*mr*ρ*π*cos(α)*sin(θ) - 16*U^2*a*α*b*c1*dr*m*mr*ρ*π*cos(α)*sin(θ) - 16*U^2*a*α*b*c3*dr*m*mr*ρ*π*cos(α)*sin(θ) + 16*U*a*b^2*c0*dr*m*mr*ρ*ξd*π*cos(α)*sin(θ) - 16*U*a*b^2*c1*dr*m*mr*ρ*ξd*π*cos(α)*sin(θ) - 16*U*a*b^2*c3*dr*m*mr*ρ*ξd*π*cos(α)*sin(θ) + 8*U^2*b*c1*c2*dr*m*mr*ρ*xbd*π*cos(α)*sin(θ) + 8*U^2*b*c3*c4*dr*m*mr*ρ*xbd*π*cos(α)*sin(θ) + 8*U^3*c1*c2*c4*dr*m*mr*ρ*xb*π*cos(α)*sin(θ) + 8*U^3*c2*c3*c4*dr*m*mr*ρ*xb*π*cos(α)*sin(θ) - 16*U*a^2*αd*b^2*c0*dr*m*mr*ρ*π*cos(α)*sin(θ) + 16*U*a^2*αd*b^2*c1*dr*m*mr*ρ*π*cos(α)*sin(θ) + 16*U*a^2*αd*b^2*c3*dr*m*mr*ρ*π*cos(α)*sin(θ) + 16*U*a*αd*b^2*c0*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) - 16*U*a*αd*b^2*c1*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) - 16*U*a*αd*b^2*c3*dr^2*m*mr*ρ*π*sin(α - θ)*cos(α) - 16*U^2*b*c1*c2*dr^2*m*mr*ρ*xbd*π*sin(α - θ)*cos(α) - 16*U^2*b*c3*c4*dr^2*m*mr*ρ*xbd*π*sin(α - θ)*cos(α) - 16*U^3*c1*c2*c4*dr^2*m*mr*ρ*xb*π*sin(α - θ)*cos(α) - 16*U^3*c2*c3*c4*dr^2*m*mr*ρ*xb*π*sin(α - θ)*cos(α) + 16*U^2*a*b*c1*c2*dr*m*mr*ρ*xbd*π*cos(α)*sin(θ) + 16*U^2*a*b*c3*c4*dr*m*mr*ρ*xbd*π*cos(α)*sin(θ) + 16*U^3*a*c1*c2*c4*dr*m*mr*ρ*xb*π*cos(α)*sin(θ) + 16*U^3*a*c2*c3*c4*dr*m*mr*ρ*xb*π*cos(α)*sin(θ))/(b*rr*(- 8*m^2*xa^2 + 8*m^2*ra^2 + 8*dr^2*m^2*mr + 8*m^2*mr*ra^2 + b^4*ρ^2*π^2 + 8*dr^2*m^2*mr^2 - 8*dr^2*m^2*mr^2*sin(θ)^2 - 8*dr^2*m^2*mr^2*sin(α - θ)^2 + b^2*m*ρ*π - 8*m^2*mr*ra^2*sin(θ)^2 - 8*dr^2*m^2*mr*sin(α - θ)^2 - 8*dr^2*m^2*mr^2*cos(α)^2 + 8*a^2*b^2*m*ρ*π + 8*b^2*m*ra^2*ρ*π + 16*dr*m^2*mr*xa*cos(α) + b^2*m*mr*ρ*π - b^2*m*mr*ρ*π*sin(θ)^2 + 16*a*b^2*m*ρ*xa*π + 16*dr*m^2*mr*xa*sin(α - θ)*sin(θ) - 16*dr^2*m^2*mr^2*sin(α - θ)*cos(α)*sin(θ) + 8*a^2*b^2*m*mr*ρ*π + 8*b^2*dr^2*m*mr*ρ*π - 8*a^2*b^2*m*mr*ρ*π*sin(θ)^2 - 8*b^2*dr^2*m*mr*ρ*π*sin(α - θ)^2 - 16*a*b^2*dr*m*mr*ρ*π*cos(α) - 16*a*b^2*dr*m*mr*ρ*π*sin(α - θ)*sin(θ)))
xbd
-(2*a*αd*b^2 - 2*b^2*ξd - 2*U*α*b - αd*b^2 + 2*U^2*c2*c4*xb + 2*U*b*c2*xbd + 2*U*b*c4*xbd)/(2*b^2)
]
end
# parameter values
par = (U = 26.0, m = 20.0, a = -0.15, b = 0.5, ra = 0.75, xa = 0.25, ωh = 2.0*π, ωα = 6.0*π, H = 7.5, ρ = 1.0, c0 = 1.0, c1 = 0.165, c2 = 0.0455, c3 = 0.335, c4 = 0.3, mr = 0.1, rr = 0.02, dr = 0.5, ln = 1.5)
# initial condition
Y0 = [1e-2, 1e-2, 1e-2, 1e-2, 1e-2, 1e-2, 1e-2, 1e-2]
# Bifurcation Problem
prob = BifurcationProblem(TS!, Y0, par, (@optic _.U);
record_from_solution = (x, p; k...) -> (ξ = x[1], α = x[3]))
# continuation of equilibria
opts_br = ContinuationPar(p_min = 26.0, p_max = 34.0)
br = continuation(prob, PALC(), opts_br;
# plot = true,
verbosity = 3)
# continuation of periodic orbits
opts_po_cont = ContinuationPar(opts_br, dsmax = 1e-2, ds = 1e-3, dsmin = 1e-4, max_steps = 1000, tol_stability = 1e-10, plot_every_step = 20)
@set opts_po_cont.newton_options.tol = 1e-10
# plot function
function plotPO(x, p; k...)
xtt = get_periodic_orbit(p.prob, x, p.p)
plot!(xtt.t, xtt[1,:]; markersize = 2, k...)
plot!(xtt.t, xtt[2,:]; k...)
plot!(xtt.t, xtt[3,:]; legend = fαse, k...)
end
# record function
function recordPO(x, p; k...)
xtt = get_periodic_orbit(p.prob, x, p.p)
period = getperiod(p.prob, x, p.p)
return (max = maximum(xtt[1,:]), min = minimum(xtt[1,:]), period = period)
end
Mt = 400
br_po = continuation(
br, 1, opts_po_cont,
PeriodicOrbitTrapProblem(M = Mt+100);#PeriodicOrbitTrapProblem(450);#PeriodicOrbitOCollProblem(450);
# plot = true,
verbosity = 3,
record_from_solution = recordPO,
plot_solution = (x, p; k...) -> begin
plotPO(x, p; k...)
## plot previous branch
plot!(br, subplot = 1, putbifptlegend = fαse)
end,
normC = norminf)
plot(br, br_po)
