Simulation of the Rutherford Experiment

Make changes in the parameters qc, xc, yc, x0, y0, time, px0, xmax and ymax as necessary.

Off[General::spell1];

This notebook uses the Runge-Kutta (RK) algorithm for the numerical integration of "Hamilton's Equations", a set of first order partial differential equations describing the motion of a massive particle in the presence of a potential field.  It then produces an animation of that motion.  The user has control over the potential, initial conditions, number of frames and delta time values (both between frames and between RK iterations).

The notebook requires the animation package.

Needs["Graphics`Animation`"];

The independent variables are defined in the list v.  They are the standard Hamiltonian phase space variables, consisting of the coordinates and momenta in each dimension.  h is the Hamiltonian (sum of kinetic and potential energies).  rk is a Runge-Kutta function which iteratively integrates Hamilton's equations of motion.  Note that we attempt to minimize the number of expressions which need to be evaluated at each step.  The spaces surrounding the replacement operator are significant.  Also, the last expression is not followed by a semicolon so that its value will be passed back to the caller (otherwise NULL is passed back).

r is defined to be the distance from the origin.  It is provided as a convenience for the user who needs to use it in the definition of the potential.

v:={x,y,px,py};

h:=(px^2+py^2)/(2*N[m]) + u;

rk[v0_]:=(a=dt*N[s /. Thread[v->v0]    ];

		b=dt*N[s /. Thread[v->v0+a/2]];
			    			
		c=dt*N[s /. Thread[v->v0+b/2]];
			    			
		d=dt*N[s /. Thread[v->v0+c]  ];
			    			
		v0+(a+2*b+2*c+d)/6
	);
									
r:=Sqrt[x^2+y^2];

qc = {79} * 1.6 * 10^-19;

xc = {0} * 10^-15;

yc = {0} * 10^-15;

x0 = -100 * 10^-15;        	y0 = 2.5 * 10^-15;

time = 1. * 10^-20;

px0 = 1.92 * 10^7 * m;

xmax = 100 * 10^-15;      ymax = 100 * 10^-15;

After changing the above values as appropriate, evaluate the entire notebook.

Red particle trajectories indicate acceleration, while blue indicate deceleration.

DO NOT CHANGE ANYTHING IN THE CELLS BELOW !

eps0 = 8.854*10^-12;
V = Sum[qc[[b]]/Sqrt[(x-xc[[b]])^2+(y-yc[[b]])^2], {b,Length[qc]}] / (4 Pi * eps0);
q = 2 * 1.6 * 10^-19;
u := q * V;
m = 6.69 * 10^-27;
t0 = 0.0;  	
py0 = 0 * 10^-10 * m;
f = 2;
inc = time/2;
ns = 200;

s={D[h,px], D[h,py], -D[h,x], -D[h,y]};
path={};
dt=N[inc/ns];
Animate[({xa,ya,pxa,pya}=Nest[rk,N[{x0,y0,px0,py0}],ns];
	AppendTo[path,If[(pxa^2+pya^2)>(px0^2+py0^2),
				RGBColor[1,0,0],RGBColor[0,0,1]]];
	AppendTo[path,Line[{{x0,y0},{xa,ya}}]];
	{ax,ay}={(pxa-px0)/m*inc,(pya-py0)/m*inc};
	{x0,y0,px0,py0}={xa,ya,pxa,pya};
	Graphics[{path},
		AspectRatio->Automatic,
		 PlotRange->{{-xmax,xmax},{-ymax,ymax}},
		 Axes->True,Ticks->None]),
	{t,t0,t0+(f-1)*inc,inc}];