Since we know that

we can solve for r'' giving

r'' = (f/μ) r/r

Seperating this into x- and y-components gives

d²x/dt² = (f/μ)*x/√(x²+y²)

d²y/dt² = (f/μ)*y/√(x²+y²)

Next, we need to find the force. Using the WOODS-SAXON spherically symmetric potential which is given by

V(r) = -V0*ff(r,R,a)

where ff(r,R,a) is a Fermi-function form factor given by

ff(r,R,a) = [1+exp(r-R/a)]^-1,

and the parameters are

Vo = 52.06 MeV

a = 0.662 fm

R = R0*A^(1/3), Ro = 1.26 fm.

Since we know that a potential is defined as V(r) = ∫ f(r)∙dr,

we can find f(r) by taking the derivative of V(r) w.r.t. r. We get

f(r) = -V0*exp(r-R/a)/{a*[1+exp(r-R/a)]}.

This is then inserted into the equations of motion giving

d²x/dt² = -V0*exp((r-R)/a)*x/{a*μ*[1+exp((r-R)/a)]*√(x²+y²)}

d²y/dt² = -V0*exp((r-R)/a)*y/{a*μ*[1+exp((r-R)/a)]*√(x²+y²)}

Solving these equations give the x- and y-coordinates of the center of mass which in turn give the x- and y-coodinates of m₁ and m₂, respectively.

To ensure unit compatibility, let

x = X*a0

y = Y*a0

a = A*a0

R = R*a0

t = T*tau

V0 = V0*Eb

m1 = M1*m0

m2 = M2*m0

μ = μ*m0

The equations of motion then become

d²X/dT² = -V0*exp((sqrt(X^2+Y^2)-R)/A)*X/{A*μ*[1+exp((sqrt(X^2+Y^2)-R)/A)]^2*√(X²+Y²)}

d²Y/dT² = -V0*exp((sqrt(X^2+Y^2)-R)/A)*Y/{A*μ*[1+exp((sqrt(X^2+Y^2)-R)/A)]^2*√(X²+Y²)}

after arbitrarily setting (Eb*tau^2 / m0*a0^2) =1.

Solving the above expression for tau, we see that the one second is equivalent to 1.02E-22 seconds of simulation time.

To determine the initial conditions, we set t=0, r_CM=0, and dr_CM/dt=v_CM=0.

This gives

r_1,0 = -(m2/m1)*r_2,0

and

v_1,0 = -(m2/m1)*v_{2,0 ,}

or in component form:

x_1,0 = -(m2/m1)*x_2,0

y_1,0 = -(m2/m1)*y_2,0

vx_1,0 = -(m2/m1)*vx_2,0

vy_1,0 = -(m2/m1)*vy_2,0

Setting either component of x=0 makes the y-component a maximum. For y we choose 1.

Similarly this would make the y-component of velocity zero, the the x-component would be a maximum.

For vx we chose 0.085.