Add quadrupole surface vibrations system (#77)

* Add quadrupole surface vibrations system

* Improve docstring

- Add more references
- Add initial conditions

Author: SebastianM-C

src/continuous_famous_systems.jl

@@ -219,6 +219,51 @@ function hhjacob!(J, u, p, t)
 end
 
 
+"""
+    qbh(u0=[]; A=1., B=0.55, D=0.4)
+
+```math
+\\begin{aligned}
+\\dot{q}_0 &= A p_0 \\\\
+\\dot{q}_2 &= A p_2 \\\\
+\\dot{p}_0 &= -A q_0 -3 \\frac{B}{\\sqrt{2}} (q_2^2 - q_1^2) - D q_1 (q_1^2 + q_2^2) \\\\
+\\dot{p}_2 &= -q_2 (A + 3\\sqrt{2} B q_1 + D (q_1^2 + q_2^2)) (x^2 - y^2)
+\\end{aligned}
+```
+
+These equations of motion correspond to a Hamiltonian used [1,2] in nuclear
+physics to study the quadrupole vibrations of the nuclear surface.
+
+```math
+H(p_0, p_2, q_0, q_2) = \\frac{A}{2}\\left(p_0^2+p_2^2\\right)+\\frac{A}{2}\\left(q_0^2+q_2^2\\right)
+			 +\\frac{B}{\\sqrt{2}}q_0\\left(3q_2^2-q_0^2\\right) +\\frac{D}{4}\\left(q_0^2+q_2^2\\right)^2
+```
+
+The Hamiltonian has a similar structure with the Henon-Heiles one, but it has an added fourth order term
+and presents a nontrivial dependence of chaoticity with the increase of energy [3].
+The default initial condition is chaotic.
+
+[1]: Eisenberg, J.M., & Greiner, W., Nuclear theory 2 rev ed. Netherlands: North-Holland pp 80 (1975)
+[2]: Baran V. and Raduta A. A., International Journal of Modern Physics E, **7**, pp 527--551 (1998)
+[3]: Micluta-Campeanu S., Raportaru M.C., Nicolin A.I., Baran V., Rom. Rep. Phys. **70**, pp 105 (2018)
+"""
+function qbh(u0=[0., -2.5830294658973876, 1.3873470962626937, -4.743416490252585];  A=1., B=0.55, D=0.4)
+    return CDS(qeom, u0, [A, B, D])
+end
+function qeom(z, p, t)
+    @inbounds begin
+        A, B, D = p
+        p₀, p₂ = z[1], z[2]
+        q₀, q₂ = z[3], z[4]
+
+        return SVector{4}(
+            -A * q₀ - 3 * B / √2 * (q₂^2 - q₀^2) - D * q₀ * (q₀^2 + q₂^2),
+            -q₂ * (A + 3 * √2 * B * q₀ + D * (q₀^2 + q₂^2)),
+            A * p₀,
+            A * p₂
+        )
+    end
+end
 
 """
     lorenz96(N::Int, u0 = rand(M); F=0.01)

test/continuous_systems.jl

@@ -64,3 +64,9 @@ end
         @test xyz[i] == data[:, i]
     end
 end
+
+@testset "Quadrupole boson Hamiltonian" begin
+    ds = Systems.qbh()
+    tr = trajectory(ds, 20.0)
+    @test size(tr) == (2001, 4)
+  end