compute_jacobi_constant Function

    pure function compute_jacobi_constant(mu,x) result(c)

    implicit none

    real(wp),intent(in)              :: mu   !! CRTBP parameter (See [[compute_crtpb_parameter]])
    real(wp),dimension(6),intent(in) :: x    !! normalized state vector
    real(wp)                         :: c    !! Jacobi constant

    !local variables:
    real,dimension(3) :: r,v,rb1,rb2
    real(wp) :: omm,r1,r2

    !extract variables from x vector:
    r = x(1:3)    ! position
    v = x(4:6)    ! velocity

    !other parameters:
    omm = one - mu
    rb1 = [-mu,zero,zero]   ! location of body 1
    rb2 = [omm,zero,zero]   ! location of body 2
    r1  = norm2(r - rb1)    ! body1 -> sc distance
    r2  = norm2(r - rb2)    ! body2 -> sc distance

    !compute Jacobi integral:
    ! [ See: http://cosweb1.fau.edu/~jmirelesjames/hw4Notes.pdf ]
    if (r1==zero .or. r2==zero) then
        c = huge(one)   ! a large value
    else
        c = r(1)**2 + r(2)**2 + &
            two*omm/r1 + two*mu/r2 - &
            (v(1)**2 + v(2)**2 + v(3)**2)
    end if

    end function compute_jacobi_constant