Hover over the menu bar to pick a physics animation.

Energy Distribution of 3D Charged Rigid Arrays In this animation, we will examine the energy distributions of 3D rigid arrays that interact with each other via potentials between the their elements and elements of other arrays. As such, these arrays will be subject to both repelling forces on their centers of mass (CM) as well as torques that tend to rotate the array about its CM. Therefore, for the 3 dimension case that we sill study here, we will have 6 degrees of freedom. These are position in the x, y and z directions as well as rotation about the (x,y,z) axes with repect to the CM. In order to compute the rotations we first compute the moments of inertia, Ix, Iy, and Iz about the three principle axes. Then we will have the equations for changes of the spin rates about these axes

`I_xdomega_x=(bbrxxbbF)bbx, I_y*domega_y=(bbxxbbF)bby, I_z*domega_z=(bbrxbbF)bbz`

where the vector cross product `(bbrxbbF)`is the torque due to force `bbF` at particle vector position `bbr` of the array. To keep the equations simple we will choose the particle locations so that the array is symmetric about its center of mass (CM). Then the expressions for moments of inertia are`I_x=sum_im_ix_i`

, etc. Having computed `domega_x`, etc., we then need to increment `omega_x`, etc., by that amount and then increment the `x` angle of the `jth` array, `Phix_j`, of the array about the `x` axis by the angle `dPhi_x=omega_x`. Then we need to compute the positions of all masses with respect to the center of mass by these equations. This is done using the (x,y,z) rotation matrices for the angle increments `(omega_x,omega_y,omega_z)` at each increment of time where the time increments are video frames. It is important to preserve the order of rotation of the `(x,y,z)` angles. Since there are 6 degrees of freedom, we should therefore expect that the energy distributions would be`D(E)=E^2*exp(-E/E_c)`

where `E_c` is some constant that is related to the average energy, E_(avg), and determined by the value of the ratio of the following integrals`E_(avg)=(intED(E)dE)/(intD(E)dE)`.

It turns out the above ratio is `E_(avg)=3E_c` which means that `E_c=E_(avg)/3` and then `D(E)` becomes`D_p(E)=CE^2exp(-(3E)/E_(avg))`

where C can be chosen so that `D(E)` represents either a probability (and integrates to 1.0) or, more conveniently for plotting, `C` can cause the peak value of `D(E)` to equal 1. To make `D(E)` equal 1 at its peak value, we differentiate the original expression with respect to E and find that the peak value of `D(E)` is at `Ep=E_c/2`. Then the value of `D(E)` at its peak is`D_(max)(Ep)=E_c^2exp(-2E_c/E_c)`

Then we choose`C=exp(2)/(4Ec*Ec)`

wich causes the peak value of Dp(E) to equal 1. So`D_p(E)=exp(2)(E^2/(4Ec^2))exp(-E/Ec)=exp(2)(9E^2/(4E_(avg)^2))exp(-3E/E_(avg))`

The plots show the distributions of the total kinetic energy (black), the translational kinetic energy (green), the rotational kinetic energy (blue), and the expression above (black) which is a smooth curve whereas the others are stepped since they are histograms. Because it takes a long time to compute any single plot, I have chosen to present some Final Distribution Pictures (Dist i..) which the learner is free to click on and view. Just close the picture window after viewing to return to the application.In this case, I have had great difficulty with the usual `1/r^2` (where r is the distance between the elements) forces between elements. I have instead used the expression `F=q^2exp(-r/r_C)` where `r_C` is a constant chosen by the learner and q represents the charge on each element.

** Over limited range of the variable parameters this force value is stable but cannot be expected to be stable for all parameters.
If an unstable set of parameters is chosen, at the left middle of the page a message 'Maximum Energy Exceeded...' will appear and you should change at least parameter and press Start again.
**