# Planar Gravitational 3-body problem
In this problem we can separate into four ODE which is
\begin{equation}
\begin{split}
\dfrac{dx}{dt} &= v_x \\
\dfrac{dy}{dt} &= v_y \\
\dfrac{dv_x}{dt} &= 2\omega v_y + \omega^2x
- \dfrac{GM_1(x+a_1)}{[(x+a_1)^2+y^2]^{3/2}}
- \dfrac{GM_2(x-a_2)}{[(x-a_2)^2+y^2]^{3/2}} \\
\dfrac{dv_y}{dt} &= -2\omega v_x + \omega^2y
- \dfrac{GM_1y}{[(x+a_1)^2+y^2]^{3/2}}
- \dfrac{GM_2y}{[(x-a_2)^2+y^2]^{3/2}}
\end{split}
\tag{1}
\end{equation}
Using normalize unit ($[t]=1/\omega$ and $[v]=a_1\omega$) we can reduce above four equation into a set of equation that only contain 3 parameter.
\begin{equation}
\begin{split}
\dfrac{dx}{d\tau} &= v_x \\
\dfrac{dy}{d\tau} &= v_y \\
\dfrac{dv_x}{d\tau} &= 2 v_y + x
- \dfrac{GM_1(x+1)}{[(x+1)^2+y^2]^{3/2}}
- \dfrac{GM_2(x-a_2)}{[(x-a_2)^2+y^2]^{3/2}} \\
\dfrac{dv_y}{d\tau} &= -2 v_x + y
- \dfrac{GM_1y}{[(x+1)^2+y^2]^{3/2}}
- \dfrac{GM_2y}{[(x-a_2)^2+y^2]^{3/2}}
\end{split}
\tag{2}
\end{equation}

## Energy Conservation
Before begin analysis we need to make sure some of the basic physical quantity is conserved, in this case we use energy as a reference for what time step are appropriate in this simulation.


First we randomly place the third mass, and we find out that their asymptotic behavior are very identical.

We then vary the initial energy, again their asymptotic are very similar.

Then we vary the parameter of $a_2$, since we found out that their asymptotic behavior are very similar, we focus on the early dynamics. We find out that with different value of $a_2$, the result can be quite different. But there is still a problem of singularity.


Same in if we vary the ratio of two stationary mass, we fond out that the early stage of the dynamics are very sensitive to parameter as well as initial condition.

But in the end what is the reason that almost every system have very similar asymptotic dynamics. I think the reason can be quite clear if we plot the potential of this restricted 3-body problem.
## Potential
The potential of this system is given by
\begin{equation}
V=-\dfrac{\omega^2}{2}(x^2+y^2)-\dfrac{GM_1}{r_1}
-\dfrac{GM_2}{r_2}.
\tag{3}
\end{equation}
Below I plot different potential for some parameter. If we only consider the effect of the potential, discard the first term in
(2) last two equation, we know that a particle tends to move to the lowest point of potential curve. With that in mind, we found out that the potential in y direction are concave down which means that the particle try to move further away in y direction. Furthermore, potential in x direction also concave down apart from the vicinity of the singularity, so their behavior are very similar to the behavior in y direction.



## Code

program sysrk4
!========================================================
!   purpose: solve system of ode
!
!   methods: rk4
!
!   input: M1, M2, a2
!
!   output: x, y, vx, vy
!
!       date        programer       description of change
!       ====        =========       =====================
!     11/17/20       morrish        original code
!========================================================
    implicit none
    double precision, external :: f1, f2, f3, f4, potential, kinetic
    double precision :: w(4, 4), h, para(3)
    double precision, dimension(:), allocatable :: x, y, vx, vy, t, U, K
    integer :: i, n
    character(len=30) :: filename

    ! read input
    write(*, *) 'M1, M2, a2:'
    read(*, *) para(1), para(2), para(3)
    write(*, *) 'filename'
    read(*, *) filename

    ! time step
    h = 1.d-5
    n = nint(10/h)
    allocate(t(n))
    allocate(x(n))
    allocate(y(n))
    allocate(vx(n))
    allocate(vy(n))
    allocate(U(n))
    allocate(K(n))
    ! parameter
    ! para(1) = 50.d0
    ! para(2) = 30.d0
    ! para(3) = 1.d0
    ! inital value
    t(1) = 0.d0
    x(1) = 0.d0
    y(1) = 5.d0
    vx(1) = 5.d0
    vy(1) = 0.d0
    U(1) = potential(x(1), y(1), para)
    K(1) = kinetic(vx(1), vy(2))

    do i = 1, n-1
        w(1, 1) = h * f1(vx(i))
        w(1, 2) = h * f2(vy(i))
        w(1, 3) = h * f3(x(i), y(i), vy(i), para)
        w(1, 4) = h * f4(x(i), y(i), vx(i), para)

        w(2, 1) = h * f1(vx(i)+.5d0*w(1, 3))
        w(2, 2) = h * f2(vy(i)+.5d0*w(1, 4))
        w(2, 3) = h * f3(x(i)+.5d0*w(1, 1), &
                         y(i)+.5d0*w(1, 2), &
                         vy(i)+.5d0*w(1, 4), &
                         para)
        w(2, 4) = h * f4(x(i)+.5d0*w(1, 1), &
                         y(i)+.5d0*w(1, 2), &
                         vx(i)+.5d0*w(1, 3), &
                         para)

        w(3, 1) = h * f1(vx(i)+.5d0*w(2, 3))
        w(3, 2) = h * f2(vy(i)+.5d0*w(2, 4))
        w(3, 3) = h * f3(x(i)+.5d0*w(2, 1), &
                         y(i)+.5d0*w(2, 2), &
                         vy(i)+.5d0*w(2, 4), &
                         para)
        w(3, 4) = h * f4(x(i)+.5d0*w(2, 1), &
                         y(i)+.5d0*w(2, 2), &
                         vx(i)+.5d0*w(2, 3), &
                         para)

        w(4, 1) = h * f1(vx(i)+w(3, 3))
        w(4, 2) = h * f2(vy(i)+w(3, 4))
        w(4, 3) = h * f3(x(i)+w(3, 1), &
                         y(i)+w(3, 2), &
                         vy(i)+w(3, 4), &
                         para)
        w(4, 4) = h * f4(x(i)+w(3, 1), &
                         y(i)+w(3, 2), &
                         vx(i)+w(3, 3), &
                         para)

        x(i+1) = x(i) + (w(1, 1) + 2.d0 * w(2, 1) + &
                         2.d0 * w(3, 1) + w(4, 1)) / 6.d0
        y(i+1) = y(i) + (w(1, 2) + 2.d0 * w(2, 2) + &
                         2.d0 * w(3, 2) + w(4, 2)) / 6.d0
        vx(i+1) = vx(i) + (w(1, 3) + 2.d0 * w(2, 3) + &
                           2.d0 * w(3, 3) + w(4, 3)) / 6.d0
        vy(i+1) = vy(i) + (w(1, 4) + 2.d0 * w(2, 4) + &
                           2.d0 * w(3, 4) + w(4, 4)) / 6.d0
        U(i+1) = potential(x(i+1), y(i+1), para)
        K(i+1) = kinetic(vx(i+1), vy(i+1))
        t(i+1) = t(i) + h
        ! write(*, '(1a5, 1f5.1)') 'left:', dble(i)/dble(n) * 100
    end do

    open(66, file=trim(filename)//'.dat', status='unknown')
    do i = 1, n
        write(66, *) t(i), x(i), y(i), vx(i), vy(i), U(i), K(i), U(i)+K(i)
    end do
    close(66)
end program sysrk4

function f1(vx)
    implicit none
    double precision :: f1, vx
    f1 = vx
end function

function f2(vy)
    implicit none
    double precision :: f2, vy
    f2 = vy
end function

function f3(x, y, vy, para)
    implicit none
    double precision :: f3, para(3), m1, m2, a2
    double precision :: x, y, vy, G=1.d0
    m1 = para(1)
    m2 = para(2)
    a2 = para(3)
    f3 = 2.d0 * vy + x - G * m1 * (x + 1.d0) / dsqrt(((x + 1.d0)**2 + y**2)**3) - &
         - G * m2 * (x - a2) / dsqrt(((x - a2)**2 + y**2)**3)
end function

function f4(x, y, vx, para)
    implicit none
    double precision :: f4, para(3), m1, m2, a2
    double precision :: x, y, vx, G=1.d0
    m1 = para(1)
    m2 = para(2)
    a2 = para(3)
    f4 = - 2.d0 * vx + y - G * m1 * y / dsqrt(((x + 1.d0)**2 + y**2)**3) - &
         - G * m2 * y / dsqrt(((x - a2)**2 + y**2)**3)
 end function

function potential(x, y, para)
    implicit none
    double precision :: potential, x, y, r1, r2, para(3), m1, m2, a2, G=1.d0
    m1 = para(1)
    m2 = para(2)
    a2 = para(3)
    r1 = dsqrt((x + 1.d0)**2 + y**2)
    r2 = dsqrt((x - a2)**2 + y**2)
    potential = - (x**2 + y**2) / 2.d0 - G * m1 / r1 - G * m2 / r2
end function

function kinetic(vx, vy)
    implicit none
    double precision :: kinetic, vx, vy
    kinetic = .5d0 * (vx**2 + vy**2)
end function