Difference between revisions of "RDYNLAB"

From New IAC Wiki
Jump to navigation Jump to search
 
(5 intermediate revisions by the same user not shown)
Line 1: Line 1:
This is a library for modeling Relativistic DYNamics in the LAB frame. It uses an eigth order Runge-Kutta algorithm to solve the relativistic dynamical equations for particle motion, including radiative power losses for charged particles, as described in [[Simulating Particle Trajectories|this page]].
+
This is a library for modeling Relativistic DYNamics in the LAB frame.  
 +
 
 +
'''The General Problem of Relativistic Particle Dynamics'''
  
 
The general problem of modeling relativistic particle trajectories for a given force boils down to solving the following system of six differential equations:
 
The general problem of modeling relativistic particle trajectories for a given force boils down to solving the following system of six differential equations:
Line 10: Line 12:
  
 
<math>\frac{\vec{v}}{c}=\frac{\vec{p}c}{E}</math>
 
<math>\frac{\vec{v}}{c}=\frac{\vec{p}c}{E}</math>
<math>E^2-c^2\vec{p}\cdot\vec{p}=m_0^2c^4</math>
+
<math>E^2-c^2\vec{p}\cdot\vec{p}=m_0^2c^4</math>,
  
 
and taking the time derivatives of both, we find
 
and taking the time derivatives of both, we find
Line 16: Line 18:
 
<math> \frac{1}{c}\frac{d\vec{v}}{dt}=\frac{c}{E}\frac{d\vec{p}}{dt}-c\vec{p}\frac{1}{E^2}\frac{dE}{dt}</math>
 
<math> \frac{1}{c}\frac{d\vec{v}}{dt}=\frac{c}{E}\frac{d\vec{p}}{dt}-c\vec{p}\frac{1}{E^2}\frac{dE}{dt}</math>
  
<math> 2E\frac{dE}{dt}-2c^2\vec{p}\cdot\frac{d\vec{p}}{dt}=0</math>
+
<math> 2E\frac{dE}{dt}-2c^2\vec{p}\cdot\frac{d\vec{p}}{dt}=0</math>.
  
 
Solving this second equation for <math>dE/dt</math> and plugging into the first, also noting that <math>d\vec{p}/dt=\vec{F}</math>, we obtain
 
Solving this second equation for <math>dE/dt</math> and plugging into the first, also noting that <math>d\vec{p}/dt=\vec{F}</math>, we obtain
 +
 +
<math> \frac{1}{c}\frac{d\vec{v}}{dt}=\frac{c}{E}\vec{F}-\frac{c}{E^3}\vec{p}\left(c^2\vec{p}\cdot\vec{F}\right)</math>.
 +
 +
Noting that <math>c\vec{p}=\frac{E}{c}\vec{v}</math>, we find
 +
 +
<math> \frac{E}{c^2}\frac{d\vec{v}}{dt}=\vec{F}-\frac{\vec{v}}{c}\left(\frac{\vec{v}}{c}\cdot\vec{F}\right)</math>.
 +
 +
We now decompose the force into pieces that are parallel and perpendicular to the velocity of the particle, i.e.
 +
 +
<math> \vec{F}=\vec{F}_{\parallel} + \vec{F}_{\perp} </math>
 +
 +
where
 +
 +
<math> \vec{F}_{\parallel}=\left(\vec{F}\cdot\hat{v}\right)\hat{v} </math>.
 +
 +
This now gives us
 +
 +
<math> \frac{E}{c^2}\frac{d\vec{v}}{dt}=\left(\vec{F}_{\parallel}+\vec{F}_{\perp}\right)-\frac{\vec{v}}{c}\left[\frac{\vec{v}}{c}\cdot\left(\vec{F}_{\parallel}+\vec{F}_{\perp}\right)\right]</math>.
 +
 +
The scalar product of the velocity and the perpendicular component of the force is zero. Furthermore the scalar product involving the velocity and the parallel component of the force, when multiplied again by the velocity, gives us <math>v^2/c^2\vec{F}_{\parallel}</math>, and thus
 +
 +
<math> \frac{E}{c^2}\frac{d\vec{v}}{dt}=\left(\vec{F}_{\parallel}+\vec{F}_{\perp}\right)-\frac{v^2}{c^2}\vec{F}_{\parallel}</math>.
 +
 +
Combining the terms with respect to the force vectors, we get
 +
 +
<math> \frac{E}{c^2}\frac{d\vec{v}}{dt}=\vec{F}_{\perp}+\left(1-\frac{v^2}{c^2}\right)\vec{F}_{\parallel}=\vec{F}_{\perp}+\frac{1}{\gamma^2}\vec{F}_{\parallel}</math>.
 +
 +
Last of all, noting that <math>E/c^2=\gamma m_0</math>
 +
 +
where <math>m_0</math> is the rest energy of the particle, and <math>\gamma</math> is the usual relativistic factor, we obtain
 +
 +
<math> \frac{d\vec{v}}{dt}=\frac{1}{\gamma m_0}\vec{F}_{\perp}+\frac{1}{\gamma^3 m_0}\vec{F}_{\parallel}</math>.
 +
 +
This completes (almost) our set of differential equations.
 +
 +
 +
'''Power Loss for Accelerated Charges'''
 +
 +
The code accounts for power radiated away as the charged particles accelerate. Such power loss is given by the Lienard formula,
 +
 +
<math>
 +
\frac{dE}{dt}=-\frac{\mu_0q^2\gamma^6}{6\pi c}\left(\vec{a}\cdot\vec{a}-\left|\frac{\vec{v}\times\vec{a}}{c}\right|^2\right)</math>
 +
 +
However, noting that <math>E=\gamma m_0 c^2</math>, this can also be written as
 +
 +
<math> \frac{d\gamma}{dt}=\frac{1}{m_0 c^2}\frac{dE}{dt}=-\frac{\mu_0q^2\gamma^6}{6\pi m_0c^3}\left(\vec{a}\cdot\vec{a}-\left|\frac{\vec{v}\times\vec{a}}{c}\right|^2\right)</math>.
 +
 +
One can determine the acceleration of the particle using the last formula in the previous section to calculate this quantity.
 +
 +
In order to implement radiative power loss into the dynamical equations, we need to go back to the relationship between force and momentum, i.e.
 +
 +
<math>\vec{F}=\frac{d\vec{p}}{dt}=\frac{d\left(\gamma m_0 \vec{v}\right)}{dt}=m_0\frac{d\gamma}{dt}\vec{v}+\gamma m_0\frac{d\vec{v}}{dt}</math>.
 +
 +
The power losses come through the term involving <math>\frac{d\gamma}{dt}</math>. Including this term into our previous differential equations yields
 +
 +
<math> \frac{d\vec{v}}{dt}=\frac{1}{\gamma m_0}\vec{F}_{\perp}+\frac{1}{\gamma^3 m_0}\vec{F}_{\parallel}-\frac{1}{\gamma}\frac{d\gamma}{dt}\vec{v}</math>.
 +
 +
 +
'''Summary of Relevant Formulae'''
 +
 +
In summary, the problem of relativistic particle dynamics reduces to the solving of the following system of differential equations:
 +
 +
<math>\frac{d\vec{x}}{dt}=\vec{v}</math>
 +
 +
<math> \frac{d\vec{v}}{dt}=\frac{1}{\gamma m_0}\vec{F}_{\perp}+\frac{1}{\gamma^3 m_0}\vec{F}_{\parallel}-\frac{1}{\gamma}\frac{d\gamma}{dt}\vec{v}</math>
 +
 +
<math> \frac{d\gamma}{dt}=\frac{1}{m_0 c^2}\frac{dE}{dt}=-\frac{\mu_0q^2\gamma^6}{6\pi m_0c^3}\left(\vec{a}\cdot\vec{a}-\left|\frac{\vec{v}\times\vec{a}}{c}\right|^2\right)</math>,
 +
 +
where the latter equation is evaluated using the acceleration without the power loss term.
 +
  
  
 +
'''Numerical Method of Solving the Differential Equations'''
  
 +
These equations are solved subject to an appropriate set of boundary conditions, usually the initial position and momentum of the particle.
  
, <math>m_0</math> is the rest energy of the particle, and <math>\gamma</math> is the usual relativistic factor. These equations are solved subject to an appropriate set of boundary conditions, usually the initial position and momentum of the particle.
+
In simple cases, such as when a charged particle travels through uniform electromagnetic fields, the set of dynamical equations can be solved analytically. In more complex cases, the system must be solved numerically.
  
In simple cases, such as when the fields are uniform, this set of equations can be solved analytically. However, for the "silver" permanent magnet we plan on using for the pair spectrometer, the fields are anything but uniform. In such a case, the system must be solved numerically.
+
The code solves the system of equations using an embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimates with adaptive stepsize control, as implemented in the GNU Scientific Library [http://www.gnu.org/software/gsl].
  
I am unable to post code to the wiki, so I have copied and pasted it all here. You will need all of these files in order to compile the code. GSL must be installed before you can compile this library.
+
I am unable to post code to the wiki, so I have copied and pasted it all [[RDYNLAB Code|here]]. You will need all of the files in order to compile the code. GSL must be installed before you can compile this library.
  
 
The code can be compiled into a shared object library using the autoconf, automake, and libtool tools. Optionally, you could just compile these source files into other programs.
 
The code can be compiled into a shared object library using the autoconf, automake, and libtool tools. Optionally, you could just compile these source files into other programs.
  
 
Return to [[Simulating Particle Trajectories]]
 
Return to [[Simulating Particle Trajectories]]

Latest revision as of 17:54, 18 June 2009

This is a library for modeling Relativistic DYNamics in the LAB frame.

The General Problem of Relativistic Particle Dynamics

The general problem of modeling relativistic particle trajectories for a given force boils down to solving the following system of six differential equations:

[math]\frac{d\vec{x}}{dt}=\vec{v}[/math]

[math]\frac{d\vec{v}}{dt}=\vec{a}[/math]

where [math]\vec{x}[/math] is the position of the particle, [math]\vec{v}[/math] is the momentum of the particle, and [math]\vec{a}[/math] is the acceleration of the particle which can be derived from the force. Classically, we would just use [math]\vec{a}=\vec{F}/m[/math]. Such a description, however, is not consistent with special relativity. Indeed, starting from the relativistic relations

[math]\frac{\vec{v}}{c}=\frac{\vec{p}c}{E}[/math] [math]E^2-c^2\vec{p}\cdot\vec{p}=m_0^2c^4[/math],

and taking the time derivatives of both, we find

[math] \frac{1}{c}\frac{d\vec{v}}{dt}=\frac{c}{E}\frac{d\vec{p}}{dt}-c\vec{p}\frac{1}{E^2}\frac{dE}{dt}[/math]

[math] 2E\frac{dE}{dt}-2c^2\vec{p}\cdot\frac{d\vec{p}}{dt}=0[/math].

Solving this second equation for [math]dE/dt[/math] and plugging into the first, also noting that [math]d\vec{p}/dt=\vec{F}[/math], we obtain

[math] \frac{1}{c}\frac{d\vec{v}}{dt}=\frac{c}{E}\vec{F}-\frac{c}{E^3}\vec{p}\left(c^2\vec{p}\cdot\vec{F}\right)[/math].

Noting that [math]c\vec{p}=\frac{E}{c}\vec{v}[/math], we find

[math] \frac{E}{c^2}\frac{d\vec{v}}{dt}=\vec{F}-\frac{\vec{v}}{c}\left(\frac{\vec{v}}{c}\cdot\vec{F}\right)[/math].

We now decompose the force into pieces that are parallel and perpendicular to the velocity of the particle, i.e.

[math] \vec{F}=\vec{F}_{\parallel} + \vec{F}_{\perp} [/math]

where

[math] \vec{F}_{\parallel}=\left(\vec{F}\cdot\hat{v}\right)\hat{v} [/math].

This now gives us

[math] \frac{E}{c^2}\frac{d\vec{v}}{dt}=\left(\vec{F}_{\parallel}+\vec{F}_{\perp}\right)-\frac{\vec{v}}{c}\left[\frac{\vec{v}}{c}\cdot\left(\vec{F}_{\parallel}+\vec{F}_{\perp}\right)\right][/math].

The scalar product of the velocity and the perpendicular component of the force is zero. Furthermore the scalar product involving the velocity and the parallel component of the force, when multiplied again by the velocity, gives us [math]v^2/c^2\vec{F}_{\parallel}[/math], and thus

[math] \frac{E}{c^2}\frac{d\vec{v}}{dt}=\left(\vec{F}_{\parallel}+\vec{F}_{\perp}\right)-\frac{v^2}{c^2}\vec{F}_{\parallel}[/math].

Combining the terms with respect to the force vectors, we get

[math] \frac{E}{c^2}\frac{d\vec{v}}{dt}=\vec{F}_{\perp}+\left(1-\frac{v^2}{c^2}\right)\vec{F}_{\parallel}=\vec{F}_{\perp}+\frac{1}{\gamma^2}\vec{F}_{\parallel}[/math].

Last of all, noting that [math]E/c^2=\gamma m_0[/math]

where [math]m_0[/math] is the rest energy of the particle, and [math]\gamma[/math] is the usual relativistic factor, we obtain

[math] \frac{d\vec{v}}{dt}=\frac{1}{\gamma m_0}\vec{F}_{\perp}+\frac{1}{\gamma^3 m_0}\vec{F}_{\parallel}[/math].

This completes (almost) our set of differential equations.


Power Loss for Accelerated Charges

The code accounts for power radiated away as the charged particles accelerate. Such power loss is given by the Lienard formula,

[math] \frac{dE}{dt}=-\frac{\mu_0q^2\gamma^6}{6\pi c}\left(\vec{a}\cdot\vec{a}-\left|\frac{\vec{v}\times\vec{a}}{c}\right|^2\right)[/math]

However, noting that [math]E=\gamma m_0 c^2[/math], this can also be written as

[math] \frac{d\gamma}{dt}=\frac{1}{m_0 c^2}\frac{dE}{dt}=-\frac{\mu_0q^2\gamma^6}{6\pi m_0c^3}\left(\vec{a}\cdot\vec{a}-\left|\frac{\vec{v}\times\vec{a}}{c}\right|^2\right)[/math].

One can determine the acceleration of the particle using the last formula in the previous section to calculate this quantity.

In order to implement radiative power loss into the dynamical equations, we need to go back to the relationship between force and momentum, i.e.

[math]\vec{F}=\frac{d\vec{p}}{dt}=\frac{d\left(\gamma m_0 \vec{v}\right)}{dt}=m_0\frac{d\gamma}{dt}\vec{v}+\gamma m_0\frac{d\vec{v}}{dt}[/math].

The power losses come through the term involving [math]\frac{d\gamma}{dt}[/math]. Including this term into our previous differential equations yields

[math] \frac{d\vec{v}}{dt}=\frac{1}{\gamma m_0}\vec{F}_{\perp}+\frac{1}{\gamma^3 m_0}\vec{F}_{\parallel}-\frac{1}{\gamma}\frac{d\gamma}{dt}\vec{v}[/math].


Summary of Relevant Formulae

In summary, the problem of relativistic particle dynamics reduces to the solving of the following system of differential equations:

[math]\frac{d\vec{x}}{dt}=\vec{v}[/math]

[math] \frac{d\vec{v}}{dt}=\frac{1}{\gamma m_0}\vec{F}_{\perp}+\frac{1}{\gamma^3 m_0}\vec{F}_{\parallel}-\frac{1}{\gamma}\frac{d\gamma}{dt}\vec{v}[/math]

[math] \frac{d\gamma}{dt}=\frac{1}{m_0 c^2}\frac{dE}{dt}=-\frac{\mu_0q^2\gamma^6}{6\pi m_0c^3}\left(\vec{a}\cdot\vec{a}-\left|\frac{\vec{v}\times\vec{a}}{c}\right|^2\right)[/math],

where the latter equation is evaluated using the acceleration without the power loss term.


Numerical Method of Solving the Differential Equations

These equations are solved subject to an appropriate set of boundary conditions, usually the initial position and momentum of the particle.

In simple cases, such as when a charged particle travels through uniform electromagnetic fields, the set of dynamical equations can be solved analytically. In more complex cases, the system must be solved numerically.

The code solves the system of equations using an embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimates with adaptive stepsize control, as implemented in the GNU Scientific Library [1].

I am unable to post code to the wiki, so I have copied and pasted it all here. You will need all of the files in order to compile the code. GSL must be installed before you can compile this library.

The code can be compiled into a shared object library using the autoconf, automake, and libtool tools. Optionally, you could just compile these source files into other programs.

Return to Simulating Particle Trajectories