Hover over the menu bar to pick a physics animation.
In order to have an ion trap, we need to compute the vectors that are normal to the boundary and which pass through any given (x,y) pair in the interior. One method is to define a potential and inward force distribution at the boundaries of the ellipse. Another method is to find the location where a particle incident on the boundary will collide with it and to reverse the normal component of its velocity. Here we chose to use the potential and force method.
The particles (discs, since this is a 2D animation) need to be kept inside the boundaries of the containing structure. The equation for the ellipse's boundaries is
`x_E^2/a^2+y_E^2/b^2=1`
or `y_E=+-bsqrt(1-(x_E/a)^2)`
In order to conserve the disc energies, when the discs hit the boundaries they must be reflected so that the normal velocity component is reversed and the tangential velocity component is not changed. An expression for the slope of the boundaries is
`(dy_E)/(dx_E)=-(+-((bx_E)/a^2)/sqrt(1-(x_E/a)^2))=-((y_Ex_E)/a^2)/(1-(x_E/a)^2)`
The `x` component of the unit vector from this is`n_x=-((dy_E)/(dx_E))/sqrt(1+((dy_E)/(dx_E))^2)`
The normal to our ellipse function described by `g(x,y)=x_E^2/a^2+y_E^2/b^2=1` is given by its gradient:`bbvecn(x,y)=grad(g(x,y))=2(x_E)/a^2bbhatx+2y_E/b^2bbhaty`
For collisions of particles inside the boundary we need the inward unit normal:`bbhatn(x,y)=(-x_E/a^2bbhatx-y_E/b^2bbhaty)/sqrt((x_E/a^2)^2+(y_E/b^2)^2)`
If we have determined that the boundary collision will occur at position (x,y) then the new velocity components will be calculated from:`b=vecn*vecv=(n_x*v_x+n_y*v_y)`
`(v_x'bbhatx,v_y'bbhaty)=(v_x(1-2n_x),v_ybbhaty(1-2bn_y)`
where the `2` is required to conserve momentum magnitude (or energy) since we've removed that incoming normal momentum component. To see this, consider the case where we have normal incidence along `bbhatn_x`. Then the result is`bbvecv'=bbvecv(1-2)=-bbvecv`
so the magnitude squared `bbvecv*bbvecv` is not changed. In general the values of `v_x` and `v_y` will be altered by different fractional amounts. It turns out that an accurate projection of a collision point on the boudary is a difficult problem. So this animation uses exponential boundary potentials and forces to contain the particles.d
`bbvecn=grad(f(x_E,y_E))=-2x_E/a^2bbhatx-(2y_E(x_E))/b^2bbhaty`
Converting this to a unit vector we have`bbhatn(x_E,y_E)=(-x_E/a^2bbhatx-(y_E(x_E))/b^2bbhaty)/sqrt((x_E/a^2)^2+((y_E(x_E))/b^2)^2)`
So that using the `bbhatx` component of the above gradient we can write:`(-x_E/a^2)/sqrt((x_E/a^2)^2+((y_E(x_E))/b^2)^2)=-((x-x_E))/sqrt((x-x_E)^2+(y-y_E(x_E)^2)^2)`
This gives us an equation in `(x,y,x_E)` that we can solve for `x_E` in terms of `(x,y)`. Since we are blessed with computers, we might use Newton's Method to solve for `x_E`. However this method is very time intensive and sometimes fails to give reasonable results. See Appendix 1 for details.What is actually done here is to first define a large array of boundary normal unit vectors as well as their roots on the boundary. Then, for each interior (x,y) pair required, we find the distance along the boundary array unit vector from its root to that (x,y) pair. The array element that we choose is the one of shortest distance, `deltas`. That represents a reasonable choice for the normal vector direction and its boundary position. Then we define to potential for this (x,y) pair as:
`V(x,y)=V_(max)exp(-(deltas)/lambda)`
where `lambda` is the `1/e` scaling length of the potential.The particle force distribution is define in an exactly analogous way.
`bbvecsF(x,y)=F_(max)exp(-(deltas)/lambda)bbhatn`
where, since the force is `bbgrad.V`, `F_(max)=-V_(max)/lambda`.The viewer may control the exponential width of the boundary potential or, equivalently, the force distribution width. The viewer may also control the number of force points along the x direction using its slider. Controls were also needed for the number of discs, the maximum potential (which regulates the maximum force) and the Initial speed of all of the discs. Radio buttons allow selection of either the display of the force or the potential distributions.
The intensity of the potential and the force distributions can be seen by comparing the color of the entity to the color bar on the left.
The animation shows that the force distribution is intense enough to keep 1000 particles trapped. These particles are repelled by the force at the boundary. The also collide with each other when sufficiently close. The kinetic energy of the particles is slightly variable because part of the kinetic energy is converted to potential energy at the boundaries. With the choice that average particle energy at the start is 1/2 you may note that the time average of the sum of the particle kinetic energies is slightly less than 500 as shown at the top left of the canvas.
notes
The basic form of 'Newton's Method' is the following:
`x'_g=x_g-f(x_g)/((df(x_g))/dx)`
where `x_g` is an initial guess and `x'_g` is the updated value of `x_g`. This process is itererated by setting `x_g=x'_g` until `|x'_g-x_g|` is small enough for the user's requirements. For our problem we have `f(x_E)`:`f(x_E)= (x_E/a^2)/sqrt((x_E/a^2)^2+((y_E(x_E))/b^2)^2)-((y-y_E(x_E))/(x-x_E))/sqrt((1+((y-y_E(x_E))/(x-x_E))^2))=0`
For the initial guess for `x_E` we might just choose `x_Eg=x`.The huge advantage of using a computer here is that we can compute the value of `(df)/dx` by making just two calculations of `f(x_g+dx/2)` and `f(x_g-dx/2)` where `dx` is very small compare to `x_g`. Then we just use the definition of a derivative:
`(df(x))/dx=(f(x+dx/2)-f(x-dx/2))/dx`
where `deltax` is about `sqrt(10^(-16)` of the value that we expect `f(x)` to be.