nuclear_binary_system_EJS.xml  by Javier E. Hasbun and Benjamin E. Hogan (2012)
For a binary system, such as the one shown below, we can use the following equations to determine the motion of m_{1} and m_{2}:
r_{1} = r_{cm}  (m_{2}/M) r r_{2} = r_{cm} + (m_{1}/M) r r_{cm} = (m_{1}r_{1}+m_{2}r_{2})/M μr'' = f r/r μ = m_{1}m_{2}/m_{1}+m_{2 } 
Since we know that 
μr'' = f r/r 
we can solve for r'' giving
r'' = (f/μ) r/r
Seperating this into x and ycomponents gives
d^{2}x/dt^{2} = (f/μ)*x/√(x^{2}+y^{2}) d^{2}y/dt^{2 }= (f/μ)*y/√(x^{2}+y^{2})
Next, we need to find the force. Using the WOODSSAXON spherically symmetric potential which is given by
V(r) = V0*ff(r,R,a)
where ff(r,R,a) is a Fermifunction form factor given by
ff(r,R,a) = [1+exp(rR/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(rR/a)/{a*[1+exp(rR/a)]}.
This is then inserted into the equations of motion giving

d^{2}x/dt^{2} = V0*exp((rR)/a)*x/{a*μ*[1+exp((rR)/a)]*√(x^{2}+y^{2})} d^{2}y/dt^{2 }= V0*exp((rR)/a)*y/{a*μ*[1+exp((rR)/a)]*√(x^{2}+y^{2})} 
Solving these equations give the x and ycoordinates of the center of mass which in turn give the x and ycoodinates of m_{1} and m_{2}, 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^{2}X/dT^{2} = V0*exp((sqrt(X^2+Y^2)R)/A)*X/{A*μ*[1+exp((sqrt(X^2+Y^2)R)/A)]^2*√(X^{2}+Y^{2})} d^{2}Y/dT^{2 }= V0*exp((sqrt(X^2+Y^2)R)/A)*Y/{A*μ*[1+exp((sqrt(X^2+Y^2)R)/A)]^2*√(X^{2}+Y^{2})}
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.02E22 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 ycomponent a maximum. For y we choose 1. Similarly this would make the ycomponent of velocity zero, the the xcomponent would be a maximum. For vx we chose 0.085. 