# The nucleon-core interaction: a nuclear physics simulation suitable for classroom use.

ABSTRACT

The orbits of a nucleon and its respective parent nucleus about their common center of mass are simulated in an effort to provide a pedagogical approach to the understanding of the structure of atomic nuclei. The nuclear exercise is treated by solving the problem with an effective force on the system's reduced mass. The potential governing the mean field is modeled by the Woods-Saxon form-factor with parameters that enable it to describe experimental findings. The Woods-Saxon potential is preferred over the infinite well and harmonic oscillator methods because both require infinite separation energies of the nucleons. The simulation is created using Easy Java Simulations (EJS) which is part of the Open Source Physics project and a MATLAB version (compatible with Octave) is included in the Appendix.

Key Words: Woods-Saxon Potential, Nucleon-Core Binary System, Shell Model, Easy Java Simulations, MATLAB, Octave

INTRODUCTION

Advancements in computing techniques have made it possible to bring about improvements in the classroom pedagogical approaches. This is particularly true in nuclear physics. The idea that nucleons interact through a nuclear force that is many times stronger and more complicated than the electromagnetic force makes it difficult to carry out simulations of their motion. It is possible, however, that with the availability of a potential that has been painstakingly investigated, and which is a good fit to nuclear interactions, such concepts can be employed in the classroom and simulations can be developed for classroom use.

Here we consider the motion of a single nucleon interacting with the core of a nucleus and develop an approach to perform the simulation of the motion while simplifying the process for classroom accessibility purposes. In studying the shell model of an atomic nucleus, the fundamental assumptions are that the motion of a single nucleon is due to the potential caused by all the other nucleons and that the individual nucleons can orbit as if they were transparent to each other, i.e., collisions do not occur (2). If we apply these two properties to the nucleus, a two-body problem emerges, which can be expressed as a one-body problem via the relative coordinate reduced mass concept often studied in classical mechanics (4). Once knowledge of the relative coordinate versus time is obtained, it is possible to have access to information regarding each individual body's position versus time.

The most commonly used shell model potential in modern day nuclear physics is known as the "Woods-Saxon potential." (1) This potential, parameterized by experimental data, provides the most realistic fit to nuclei across the chart of the nuclides (1). Using this potential and the concept of reduced mass, we derive the equations of motion. We then solve the resulting equations numerically with MATLAB's ode45 nonstiff differential equation solver. This is a fourth and fifth order Runge-Kutta method which is the first recommended method to try when solving numerical differential equations. We have thus developed a simulation of the nucleon and its parent nucleus orbiting their common center of mass.

The Binary System and the Concept of Reduced Mass

In this section, for pedagogical reasons, we review relative coordinate, reduced mass approach and show our implementation scheme. We consider a nucleus composed of A nucleons; that is, Z protons and A-Z neutrons. One of the nucleons, with mass [m.sub.2], is removed leaving the parent nucleus with A-1 nucleons and mass [m.sub.2]. We then consider these two particles interacting through a nuclear force. The force itself will later be shown to have origins in a potential of a Woods-Saxon form (2). The nuclear force is an attractive force which is exerted by each mass onto the other. Considering Figure 1, the coordinate of [m.sub.1] is [r.sub.1] while the coordinate of [m.sub.2] is [r.sub.2].

The force on [m.sub.2] due to [m.sub.1] is [f.sub.21], and the force on m due to [m.sub.2] is [f.sub.12]. Therefore the force on each mass is, according to Newton's 2nd law is,

[f.sub.12] = [m.sub.1][d.sup.2][r.sub.1]/d[t.sup.2] (1)

[f.sub.21] = [m.sub.2][d.sup.2][r.sub.2]/d[t.sup.2] (2)

In the binary system, the forces [f.sub.21] and [f.sub.21] point along the vector connecting [m.sub.1] and [m.sub.2]and, because of Newton's 3rd law, we have that [f.sub.21] = -[f.sub.21], where

[f.sub.12] = [f.sub.12] [[^.r].sub.12] = [f.sub.12] [r.sub.12]/[r.sub.12] (3)

[f.sub.21] = [f.sub.21] [[^.r].sub.21] = [f.sub.21] [r.sub.21]/[r.sub.21] (4)

and where we have used the definitions

[r.sub.21] = [r.sub.2] - [r.sub.1] (5)

[r.sub.12] = [r.sub.1] - [r.sub.2] (6)

We see that [r.sub.12] = -[r.sub.21] and [r.sub.12] = [r.sub.21], therefore

[f.sub.12] = -[f.sub.12] [r.sub.21]/[r.sub.21] (7)

Equation (7) permits us to rearrange Equations (3) and (4) as

[r.sub.21]/[r.sub.21] = -[f.sub.21]/[f.sub.12] = [f.sub.21]/[f.sub.12] (8)

which is consistent with Newton's 3rd law and showing that for the magnitudes [f.sub.12] = [f.sub.21] Since the magnitudes of the forces are equal, we will refer to them as a single central force f later.

Combining Equations (1) and (2) with Equations (3) and (4) gives

[m.sub.1][d.sup.2][r.sub.1]/d[t.sub.2] = [f.sub.12]/[r.sub.12] [r.sub.12] (9)

and

[m.sub.2][d.sup.2][r.sub.2]/d[t.sub.2] = [f.sub.21]/[r.sub.21] [r.sub.21]. (10)

The center of mass for a system of two particles is

[r.sub.CM] = ([m.sub.1][r.sub.1] + [m.sub.2][r.sub.2])/M, (13)

where M = [m.sub.1] + [m.sub.2]; i.e., the total mass. We also define the relative coordinate variable as

r [equivalent to] [r.sub.21] = [r.sub.2] - [r.sub.1]. (12)

With the above definitions we can write Equations (11) and (12) in form as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (13)

which with the help of the inverse of the 2x2 matrix,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

allows us to express [r.sub.1] and [r.sub.2] in terms of [r.sub.CM] and r as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (14)

or

[r.sub.1] = [r.sub.CM] - [m.sub.2]/M r (15)

and

[r.sub.2] = [r.sub.CM] - [m.sub.1]/M r. (16)

We next multiply Equation (9) by [m.sub.2] and Equation (10) by [m.sub.1] to get

[m.sub.1][m.sub.2][d.sup.2][m.sub.1]/d[t.sub.2] = [m.sub.2][f.sub.12]/[r.sub.12][r.sub.12], (17)

[m.sub.1][m.sub.2][d.sup.2][m.sub.1]/d[t.sub.2] = [m.sub.1][f.sub.21]/[r.sub.21][r.sub.21]. (18)

We then subtract Equation (17) from Equation (18) and obtain

[m.sub.1][m.sub.2][d.sup.2][m.sub.1]/d[t.sub.2][r.sub.2] - [r.sub.1] = [m.sub.1][f.sub.21]/[r.sub.21][r.sub.21] - [m.sub.2][f.sub.12]/[r.sub.12][r.sub.12]. (19)

Using the fact that r = [r.sub.2] - [r.sub.1] = [r.sub.12] = -[r.sub.21] so that [r.sub.21]= [r.sub.12] = and that the magnitudes [f.sub.21] = [f.sub.12] = f, we can rewrite Equation (19) as

[m.sub.1][m.sub.2][d.sup.2][m.sub.1]/d[t.sub.2] = f/r r([m.sub.1] + [m.sub.2]). (20)

or

[mu] [d.sup.2]r/d[t.sup.2] = f [^.r] (21)

where we have used r = r[^.r] and the definition for the reduced mass

[m.sub.1][m.sub.2]/[m.sub.1] + [m.sub.2] = [mu]. (22)

Solving the single Equation (21) for the motion of the reduced mass [??] with relative coordinate r is equivalent to solving Equations (9) and (10) separately for the motion of [m.sub.1] and [m.sub.2] with respective coordinates [r.sub.1] and [r.sub.2]. This is so because, given that the center of mass of the system, once r is known we can use Equations (15) and (16) to obtain the separate coordinates.

In order to develop a program code for classroom simulation purposes, we will work with the following two-dimensional Cartesian coordinate system:

r = x[^.[iota]] + y[^.J], (23)

r = [square root of ([x.sup.2] + [y.sup.2])], (24)

[^.r] = r/r = x/[square root of ([x.sup.2] + [y.sup.2])][^.[iota]] + y/[square root of ([x.sup.2] + [y.sup.2])][^.J], (25)

so that Equation (21) becomes

[d.sup.2]/d[t.sup.2](x[^.[iota]] + y[^.j] = f/[mu](x/[square root of ([x.sup.2] + [y.sup.2])][^.[iota]] + y/[square root of ([x.sup.2] + [y.sup.2])][^.J]), (26)

or for each coordinate

[d.sup.2]x/d[t.sup.2] = f/[mu](x/[square root of [x.sup.2] + [y.sup.2])], (27)

[d.sup.2]y/d[t.sup.2] = f/[mu] y/[square root of ([x.sup.2] + [y.sup.2])]. (28)

We will be seeking the numerical solutions equations. From these, with [r.sub.CM] = [x.sub.CM][^.[iota]] + [y.sub.CM]^j, we will obtain the individual coordinates of Equations (15) and (16) as follows.

[x.sub.1] = [x.sub.CM] - [m.sub.2]/([m.sub.2] + [m.sub.2])x (29)

[y.sub.1] = [y.sub.CM] - [m.sub.2]/([m.sub.2] + [m.sub.2])y (30)

[x.sub.2] = [x.sub.CM] - [m.sub.2]/([m.sub.2] + [m.sub.2])x (31)

[y.sub.2] = [y.sub.CM] - [m.sub.2]/([m.sub.2] + [m.sub.2])y (32)

To summarize, our approach is to solve Equations (27) and (28) numerically as a function of time. From these solutions, we will apply Equations (29-32) to obtain the motion for the binary system's individual masses [m.sub.1] and [m.sub.2]. Before we can effect this process, however, we need to develop the force f that's responsible for the nucleon-nucleus core interaction. This brings us directly to the Woods-Saxon potential.

The Woods-Saxon Form Factor

The much stronger nuclear force overcomes the repulsive Coulomb force between the constituent protons and thus the Coulomb interaction is ignored. The protons are accompanied by a number of neutrally charged particles--the neutrons. Unlike the gravitational and Coulomb forces, the strong force does not have a definitely well-known simple form. The Woods-Saxon potential models the nuclear mean field of the atomic nucleus, which is governed by the Fermi-function form factor (1); that is,

f(r, R, a) = [[1 + exp(r-R/a)].sup.-1]. (33)

where the size R and diffuseness of the surface a are fixed parameters. The total nuclear potential is given by

V(r) = -[V.sub.0]f(r, R, a), (34)

where [V.sub.0] is the total potential strength and the minus sign indicates its attractive nature. Since

V(r) = -[integral] f(r) * dr. (35)

we know that

f(r) = -dv(r)/dr, (36)

which, with the use of Equations (33) and (34) into Equation (36), yields the Woods-Saxon nuclear force

f(r) = -[v.sub.0][[1 + exp(r-R/a)].sup.-2] exp(r-R/a). (37)

Substituting this force, into Equations (27) and (28), we find the Cartesian form for each component

[d.sup.2]x/d[t.sup.2] = [v.sub.0]/[[mu]a[1 + exp([square root of ([x.sup.2] + [y.sup.2]] - R/a)].sup.-2]exp([square root of ([x.sup.2] + [y.sup.2]] - R/a) x/[square root of ([x.sup.2] + [y.sup.2])], (38)

[d.sup.2]y/d[t.sup.2] = [v.sub.0]/[[mu]a[1 + exp([square root of ([x.sup.2] + [y.sup.2]] - R/a)].sup.-2]exp([square root of ([x.sup.2] + [y.sup.2]] - R/a) x/[square root of ([x.sup.2] + [y.sup.2])], (39)

The universal values of the parameters R. a, and Vo provided by Schwierz et. al. (1) are

[v.sub.0] = 52.06 MeV, (40)

a = 0.662 fm, (41)

R = [R.sub.0] [A.sup.1/3], (42)

where

[R.sub.0] = 1.26 fm. (43)

In calculations of this type and to ensure proper functionality of the simulating program, all quantities should be expressed in appropriate dimensions. We define our distances in units of [a.sub.0] = 1fm; that is,

x = [bar.x][a.sub.0] y = [bar.y][a.sub.0] a = [bar.a][a.sub.0] and R = [bar.R][a.sub.0], (44)

with the barred quantities being dimensionless throughout this paper. Time will be expressed in units of -r and energy in units of [E.sub.b] = 1MeV, to write

t = [bar.t][tau] and [V.sub.0] = [bar.[V.sub.0]][E.sub.b]. (45)

The units of mass will be given in terms of the proton mass [m.sub.0] = 938MeV, to have

[m.sub.1] = [bar.[m.sub.1]][m.sub.0] [m.sub.2] = [bar.[m.sub.2]][m.sub.0], (46)

and

[mu] = [m.sub.1][m.sub.2]/[m.sub.1] + [m.sub.2] = [bar.[m.sub.1]][bar.[m.sub.2]]/[bar.[m.sub.1]] + [bar.[m.sub.2]][m.sub.0] = [bar.[mu]][m.sub.0]. (47)

Since we are dealing with a nucleon-core interaction where the nucleon mass is [m.sub.nucleon] the values of [bar.[m.sub.1]] and [bar.[m.sub.2]] are

[bar.[m.sub.1]] = [m.sub.nucleon]/[m.sub.0], and [bar.[m.sub.2]] = [m.sub.2]/[m.sub.0]. (48)

By making these replacements into Equations (38) and (39) and grouping units we have

[d.sup.2][bar.x]/d[t.sub.2] = -([E.sup.b][[tau].sub.2]/[m.sub.0][a.sub.0.sup.2][bar.[v.sub.0]]/[bar.[mu][bar.a]][1 + exp([square root of ([[bar.x].sup.2] + [[bar.y].sup.2] - [bar.R])]/[bar.a])].sup.-2]exp([square root of ([[bar.x].sup.2] + [[bar.y].sup.2] - [bar.R])]/[bar.a])[.bar.x]/[square root of ([[bar.x].sup.2] + [[bar.y].sup.2] - [bar.R])], (49)

Arbitrarily we next let

[d.sup.2][bar.y]/d[t.sub.2] = -([E.sup.b][[tau].sub.2]/[m.sub.0][a.sub.0.sup.2][bar.[v.sub.0]]/[bar.[mu][bar.a]][1 + exp([square root of ([[bar.x].sup.2] + [[bar.y].sup.2] - [bar.R])]/[bar.a])].sup.-2]exp([square root of ([[bar.x].sup.2] + [[bar.y].sup.2] - [bar.R])]/[bar.a])[.bar.y]/[square root of ([[bar.x].sup.2] + [[bar.y].sup.2] - [bar.R])]. (50)

([E.sup.b][[tau].sub.2]/[m.sub.0][a.sub.0.sup.2]) = 1, (51)

so that Equations (49) and (50) yield the final equations that will be used in the simulation: that is,

[d.sup.2][bar.x]/d[t.sup.2] = -[V.sub.0]/[bar.u][bar.a][[1+exp(([square root of ([[bar.x].sup.2] + [[bar.y].sup.2])-[bar.R])/[bar.a])].sup.-2]exp(([square root of ([[bar.x].sup.2] + [[bar.y].sup.2])-[bar.R])/[bar.a])][bar.x]/[square root of ([[bar.x].sup.2] + [[bar.y].sup.2], (52)

[d.sup.2][bar.y]/d[t.sup.2] = -[V.sub.0]/[bar.u][bar.a][[1+exp(([square root of ([[bar.x].sup.2] + [[bar.y].sup.2])-[bar.R])/[bar.a])].sup.-2]exp(([square root of ([[bar.x].sup.2] + [[bar.y].sup.2])-[bar.R])/[bar.a])][bar.x]/[square root of ([[bar.x].sup.2] + [[bar.y].sup.2], (53)

It is useful to know the time scale of the simulation. We can solve (51) for [tau] to find

[tau] = [square root of ([[m.sub.0][a.sub.0.sup.2])/[E.sub.b] = [square root of((938MeV/[c.sup.2])(1f[m.sup.2])/1MeV]) = [square root of (938[(fm/c).sup.2])] = 30.6 fm/c = 1.02 x [10.sup.-22]s (54)

From this, we see that each second of simulation time is approximately equivalent to about [1O.sup.-22] seconds in real time. This is indeed the reason why

dimensionless units are used.

Initial Conditions

To get the simulation process started, we need the initial relative coordinate r and the relative velocity [*.r] = v at t=0. The relative and the center of mass coordinates are given by Equations (11) and (12); that is,

[r.sub.CM] = ([m.sub.2][r.sub.2] + [m.sub.2][r.sub.2])/M, r [equivalent to] [r.sub.21] = [r.sub.2]-[r.sub.1]. (55)

We also need the velocities

(d/dt)r = d/dt([r.sub.2] - [r.sub.1]) = v = [v.sub.2] - [v.sub.1], (56)

Similarly, for the center of mass velocity (d [r.sub.CM]/dt),

[v.sub.CM] = ([m.sub.1][v.sub.1]+ [m.sub.2][v.sub.2])/M. (57)

With these, if we pick the initial coordinates [r.sub.1,0] and [r.sub.2,0], then the relative coordinate and center of mass positions at t=0 are

[r.sub.0] = [r.sub.2,0] - [r.sub.1,0] (58)

[r.sub.CM,0] = [m.sub.2][r.sub.2,0] + [m.sub.2][r.sub.2,0]/M (59)

Similarly, letting the initial particles velocities be [v.sub.1,0] and [v.sub.2,0], we have for the relative and center of mass velocities

[v.sub.0] = [v.sub.2,0] - [v.sub.1,0], (60)

[v.sub.CM,0] = ([m.sub.1][v.sub.1,0] + [m.sub.2][v.sub.2,0])/M. (61)

We further let the center of mass be located at the origin with zero velocity; that is,

[r.sub.CM,0] = 0, [v.sub.CM,0] = 0. (62)

With this understanding, Equation (59) gives

[m.sub.1][r.sub.1,0] = -[m.sub.2][r.sub.2,0] (63)

and Equation (61) gives

[m.sub.1][r.sub.1,0] = -[m.sub.2][r.sub.2,0] (64)

We see that Equation (63) yields the initial position of [m.sub.1] as

[r.sub.1,0] = -[m.sub.2]/[m.sub.1] [r.sub.2,0] (65)

or in Cartesian coordinates

[x.sub.1,0] = -[m.sub.2]/[m.sub.1] [x.sub.2,0], [y.sub.1,0] = -([m.sub.2]/[m.sub.2] [y.sub.2,0] (66)

Likewise, from Equation (64) we have

[v.sub.x1,0] = -([m.sub.2]/[m.sub.1] [v.sub.x2,0]), [v.sub.y1,0] = -([m.sub.2]/[m.sub.1] [v.sub.y2,0]) (67)

Equations (66) and (67) consist of our initial conditions for the masses. We essentially choose [x.sub.2,0], [y.sub.2,0], [v.sub.x2,0] and [v.sub.y2,0] then use the formulas above to solve for the unknowns. These are then inserted into the relative coordinate formulas of Equations (58) and (60) giving an initial position and velocity for the reduced mass orbit of Equations (52) and (53); that is,

[v.sub.x,0] = [v.sub.x2,0] - [v.sub.x1,0], [v.sub.y,0] = [v.sub.y2,0] - [v.sub.y1,0], (68)

[x.sub.0] = [x.sub.2,0] - [x.sub.1,0], [y.sub.0] = [y.sub.2,0] - [y.sub.1,0]. (69)

Simulations and Results

The actual simulation involves the numerical solutions of the Equations (52) and (53). The sample nucleus treated in our code is helium, but other nuclei can be considered equally well by changing the value of the nucleus' mass as well as the revolving particle's mass (neutron or proton). As mentioned before, the nucleus is considered as a system of two particles. One of the particles is a nucleon ([m.sub.1]) that's considered to be interacting with the second particle; i.e., the remaining core ([m.sub.2]). In our simulation the revolving particle is a neutron, thus the remaining core consists of two protons and a neutron. We have carried out the solutions in Java (using EJS) and MATLAB or Octave. The implementation is such that one writes

d[bar.x]/dt = [[bar.v].sub.x], d[[bar.v].sub.x]/dt = [[bar.a].sub.x] = [bar.f]/[mu][bar.r] [bar.x], d[bar.y]/dt = [[bar.v].sub.y], d[bar.y]/dt = [[bar.a].sub.y] = [bar.f]/[mu][bar.r] [bar.y] (70)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and [bar.r] = [square root of ([[bar.x].sup.2] + [[bar.y].sup.2])]. So that we actually solve four first order differential equations with the initial conditions outlined in the previous sections. A plot of [bar.f] versus [bar.r] is shown in Figure 2.

The solution of these equations yields the sought reduced mass coordinates [bar.x]([bar.t]) and [bar.y]([bar.t]), which we plot below in Figure 3. Included in the plot are the coordinate positions [r.sub.1] and [r.sub.2] for the masses [m.sub.1] and [m.sub.2], respectively, using the dimensionless units discussed previously. The figure's legend shows the coordinates plotted. The dotted black curve represents the coordinates of the reduced mass versus time. The middle sized orbit represents the motion of [m.sub.1] (the neutron) while the smaller sized orbit refers to the motion of [m.sub.2] (the remaining core). The center of mass is located at the origin. For colors refer to the code provided in the appendix or the website of reference (5).

To produce Figure 3, in addition to the conditions mentioned before, the initial values of [x.sub.2] = 0, [y.sub.2] = 1.0, [v.sub.y2] = 0, and [v.sub.x2] = 0.085 were also used. The two masses are seen to orbit about their common center of mass, from which the larger mass deviates less, producing a smaller orbit. The number of revolutions of either mass before returning to their respective starting positions is determined by the initial velocities. A larger tangential velocity will result in fewer revolutions.

The MATLAB/Octave code used to carry out the simulation is reproduced in the Appendix and the results are shown on the left panel of the figure. The EJS (3), (5) code was used to create the figure's right panel. While the right panel of Figure 3 only shows the resulting plot, the EJS simulation contains a panel with which the user can change the atomic mass number and mass in order to analyze the orbitals of other atoms. The MATLAB simulation can be altered in the same way by simply changing the same values in the code.

CONCLUSION

In this paper, the orbits of a nucleon and its respective parent nucleus about their common center of mass has been simulated in an effort to provide a pedagogical approach to the understanding of the structure of atomic nuclei. We have treated the problem by modeling the interaction between the particles using an effective force derived from a Woods-Saxon potential. The parameters used in the simulation are those proposed in a recent study (1). We have presented full details of the simulation to serve teachers and students alike. The approach is suitable for classroom use for those interested in learning more about nuclear physics. We have developed Java code that will be made accessible in the future (5) and which runs using the EJS open source project (3). A suitable version of the simulation has also been developed using MATLAB (6) code, which is compatible with Octave (7). This has been included in the Appendix.

ACKNOWLEDGEMENTS

We thank Francisco Esquembre (3) and the Open Source Physics project for the use of the EJS software. The OSP project can be found at http://www.opensourcephysics.org/.

APPENDIX

This appendix contains the MATLAB/Octave code "nuclear_binary_system.m" that reproduces the text's Woods-Saxon nucleon-nucleus-core force shown in Figure 2 as well as the simulation results shown in Figure 3.
```
%nuclear_binary_system_matlab_octave.m--by Javier E. Hasbun and
%Benjamin E. Hogan (2012)
%Program to simulate the nucleon-core motion as a binary system
%The code is setup to run compatibly with Matlab as well as Octave.
%Matlab (commercial) is found here http://www.mathworks.com
%Octave (open source) is found here http://www.gnu.org/software/octave/

function nuclear_binary_system_matlab_octave
clear;
global a VO R mu
mp=1.00727647; %proton mass in atomic mass units
mn=1.00866492; %neutron mass in atomic mass units
m1 = mn/mp; %ml in units of one proton (the neutron here)
A= 4; %atomic mass number (A=4 => helium)
m2 = 4.002603-mn; %core: mass of He nucleus minus one neutron
m2 = m2/mp; %m2 in "proton units"
M = ml+m2; %total mass
VO= 52.06; %potential of the system (units of Eb=1 MeV)
a= 0.662; %skin thickness of the nucleus (units of fm)
R= 1.26 * A^(1/3); %radius of the nucleus (units of fm)
mu= ml * m2/M; %reduced mass (units of mp)
%center of mass set at zero
C[M.sub.x] = 0;
C[M.sub.y] = 0;
%initial positions for ml and m2
x2 = 0;
xl = (-m2/m1)*x2;
y2=1.0;
yl = (-m2/m1)*y2;
vy2 = 0;
vyl = (-m2/m1)*vy2;
vx2 = 0.085;
vxl = (-m2/m1)*vx2;
%initial conditions for the reduce mass orbit
vx0 = vx2-vx1;
vy0 = vy2-vy1;
x0 = x2-x1;
y0 = y2-y1;
ic=[x0;vx0;y0;vy0]; %initial conditions
tmax=50; %max simulation time in units of tau (tau_0~1e-20 s)
%Use Runge-Kutta (4,5) formula with the following settings
opt=odeset('AbsTor ,1.e-7,'RelTol',1.e-4,'InitialStep',3,'MaxStep',5);
[t,w]=ode45(@nuclear_binary_der,[0.0,tmax],ic,opt);
%the calculated reduced mass coordinates and velocities
x=w(:,1);
vx=w(:,2);
Y=4,3);
vy=w(:,4);
[nr,nc]=size(w); %get rows and columns of r
%Uncomenting the lines in the following loop allow real time view
%but works best only in MATLAB.
% figure
% for i=1:nr
% clf;
% plot(CMx,CMy,'.')% center of mass
% hold on
% plot(x,y,'k:')
% plot(x(i),y(i),'k.','MarkerSize',20)
% pause(0.5)
% end
figure
plot(CMx,CMy,'k.')% center of mass
hold on
plot(x,y,'k:')
xi = CMx - m2*x/M;
yl = CMy - m2*y/M;
x2 = CMx + ml*x/M;
y2 = CMy + m1y*/M;
plot(xl,y1,'r-')
plot(x2,y2,'b--')
xlabel('x (fm)','FontSize',14)
ylabel('y (fm) ,'FontSize',14)
str=cat(2,'Nuclear Binary System m_1 = num2str(m1,3),' u, m_2 = ', ...
num2str(m2,3),' u');
title(str,'FontSize',12)

%Plot of the Woods-Saxon nuclear force
figure
r=0:01:5;
plot(r,woods_saxon(r,a,VO,R),'k )
xlabel('r (fm)','FontSize',12)
ylabel(f(r) (MeV/fm)','FontSize',12)
title('Woods-Saxon nuclear force versus r','FontSize',12)

%onuclear_binary_der: returns the derivatives for the binary mass
system
function derivs = nuclear_binary_der(t,w)
%This calculates the derivatives:
%vx=dx/dt. ax=dvx/dt, vy=dy/dt and ay=dvy/dt
%Entries in the vector of dependent variables are:
%w(1,2,3,4)=x,vx,y,vy
global a VO R mu
x=w(1);
vx=w(2);
y=w(3);
vy=w(4);
r=sqrt(x. ^2+y.^2);
fw=woods_saxon(r,a,VO,R)/mu;
derivs = [vx;fw.*x./r;vy;fw.*y./r]; %contains the derivatives

function fws=woods_saxon(r,a,VO,R)
%The Woods-Saxon nuclear force
ep=exp((r-R)/a);
eps=a*(1+ep).^2;
fws=-V0*ep./eps;
```

REFERENCES

(1.) Schwierz N, Wiedenhover I, and Volya A: "Parameterization of the Woods-Saxon Potential for Shell-Model Calculations." Cornell University Library, p5, 2007. The manuscript is found here: http://arxiv.org/pdf/0709.3525v1.

(2.) Krane KS: Introductory Nuclear Physics. John Wiley & Sons, NY, p118-121, 1988.

(4.) Hasbun, Javier E: Classical Mechanics with MATLAB Applications. Jones and Bartlett, Sudbury, MA, 2009.

(5.) The EJS open source code to run the Java version of the simulation will be available in the future through this website: http://www.westga.edu/~jhasbun/osp/osp.htm (item 25); or directly here: http://www.westga.edu/~jhasbun/osp/nuclear_binary_system_EJS.html.

(6.) MATLAB. Natick, MA: The Mathworks Inc., 2012. Available at http://www.mathworks.com.

(7.) GNU Octave. John W Eaton, 1998-2012. Available at http://www.gnu.org/software/octave/.

Benjamin E. Hogan and Javier E. Hasbun *

Department of Physics

University of West Georgia

Carrollton, GA 30118

* Corresponding Author

E-mail: jhasbun@westga.edu