<jupyter_start><jupyter_text>Electron equation of motion due to Lorentz force and [Ford&O'Connell1991](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.524.6982&rep=rep1&type=pdf) form of the Abraham-Lorentz force\begin{align*}m_e \frac{d \vec{v}}{d t} &= \vec{F}_{Lorentz} + \frac{q^2}{6 \pi \epsilon_0 c^3 m_e}\left(\frac{d \vec{F}_{Lorentz}}{dt}\right)\\&= \vec{F}_{Lorentz} + \tau\left(\frac{d \vec{F}_{Lorentz}}{dt}\right)\end{align*}Note: Above equation is in SI. Ford and O'Connell is in CGS.<jupyter_code>%run config.py
from ford1991 import Ford1991Solver
from utils import error_plot<jupyter_output><empty_output><jupyter_text>Rewrite RHS of Ford 1991 from:\begin{align*}m \frac{dv_x}{dt} &= q v_y B_z + q B_z \tau \frac{d}{dt}v_y \\m \frac{dv_y}{dt} &= - q v_x B_z - q B_z \tau \frac{d}{dt}v_x \\\end{align*} To:\begin{align*}\frac{d v_x}{dt} &= \left(\omega v_y - \omega^2 \tau v_x\right)\left(\frac{1}{1 + \omega^2 \tau^2}\right)\\\frac{d v_y}{dt} &= \left(-\omega v_x - \omega^2 \tau v_y\right)\left(\frac{1}{1 + \omega^2 \tau^2}\right)\end{align*}Where $\omega = \frac{qB_z}{m}$ is the ***signed*** angular frequency.  Writing $\mu = \tau \omega^2$, the velocity solution is:\begin{align*}v_x &= \frac{e^{-\mu t} \sin(\omega t)}{1 + \tau^2 \omega^2} \\v_y &= \frac{e^{-\mu t} \cos(\omega t)}{1 + \tau^2 \omega^2}\end{align*} Integrate to obtain position solutions:\begin{align*}x &= \frac{-e^{-\mu t} \left(\mu\sin(\omega t) + \omega \cos(\omega t)\right)}{\left(\omega^2 + \mu^2\right)\left(1 + \tau^2 \omega^2\right)} \\y &= \frac{e^{-\mu t} \left(\omega\sin(\omega t) - \mu \cos(\omega t)\right)}{\left(\omega^2 + \mu^2\right)\left(1 + \tau^2 \omega^2\right)}\end{align*}<jupyter_code># Problem parameters. Not real values
B = 1.5
m = 3
q = -1

# Note, this is a signed value
omega = B * q / m

tau = 0.025 # Not the real value. Chosen to test scheme

solver = Ford1991Solver(charge=q, mass=m, b_field=B, tau=tau)
# Solve equation numerically
res = solver.solve_1d(5)

# Get analytic solution
x_soln, y_soln, vx_soln, vy_soln = solver.analytic_solution_1d(res.t, res.y[0:2,0], res.y[2:4,0])

# Plot result
error_plot(res.y[0], res.y[1], x_exact=x_soln, y_exact=y_soln,
          title='Electron Trajectory', xlabel='x', ylabel='y')
# Kinetic energy of electron as function of time
ke = 0.5 * m * (res.y[2]**2 + res.y[3]**2)
# Analytic solution
mu = omega**2 * tau
taue = - 2.0 * mu

error_plot(res.t, ke, y_exact=ke[0]*np.exp(taue * res.t), title='Energy decay',
          xlabel='t', ylabel='Electron energy')
et = ke + res.y[4]
error_plot(res.t, np.abs(et - et[0]) / et[0], xlabel='t', ylabel='Fractional Error',
          title='Energy Conservation Error')
error_plot(res.t, res.y[2], y_exact=vx_soln, xlabel='t', ylabel=r'$v_x$')
error_plot(res.t, res.y[3], y_exact=vy_soln, xlabel='t', ylabel=r'$v_y$')
error_plot(res.t, res.y[0], y_exact=x_soln, xlabel='t', ylabel=r'$x$')
error_plot(res.t, res.y[1], y_exact=y_soln, xlabel='t', ylabel=r'$y$')
rad = np.sqrt(res.y[1]**2 + res.y[0]**2)
rad_exact = np.sqrt(y_soln**2 + x_soln**2)
error_plot(res.t, rad, y_exact=rad_exact, xlabel='t', ylabel=r'radius')
# Error in x,y as a function of time
plt.plot(res.t, res.y[0] - x_soln)
plt.plot(res.t, res.y[1] - y_soln)
# Error in v_x, v_y as a function of time
plt.plot(res.t, res.y[2] - vx_soln)
plt.plot(res.t, res.y[3] - vy_soln)
# Change v0
v0 = 3.15
res = solver.solve_1d(5, v0=v0)

# Get analytic solution
x_soln, y_soln, vx_soln, vy_soln = solver.analytic_solution_1d(res.t, res.y[0:2,0], res.y[2:4,0])

# Plot result
error_plot(res.y[0], res.y[1], x_exact=x_soln, y_exact=y_soln,
          title='Electron Trajectory', xlabel='x', ylabel='y')
# Error in v_x, v_y as a function of time
plt.plot(res.t, res.y[2] - vx_soln)
plt.plot(res.t, res.y[3] - vy_soln)
# Error in x,y as a function of time
plt.plot(res.t, res.y[0] - x_soln)
plt.plot(res.t, res.y[1] - y_soln)
# Change v0 direction
vel0 = [1, 1]
res = solver.solve_1d(5, v0=vel0)

# Get analytic solution
x_soln, y_soln, vx_soln, vy_soln = solver.analytic_solution_1d(res.t, res.y[0:2,0], res.y[2:4,0])

# Plot result
error_plot(res.y[0], res.y[1], x_exact=x_soln, y_exact=y_soln,
          title='Electron Trajectory', xlabel='x', ylabel='y')
# Error in x,y as a function of time
plt.plot(res.t, res.y[0] - x_soln)
plt.plot(res.t, res.y[1] - y_soln)
# Change v0 direction
vel0 = [1, 0]
res = solver.solve_1d(5, v0=vel0)

# Get analytic solution
x_soln, y_soln, vx_soln, vy_soln = solver.analytic_solution_1d(res.t, res.y[0:2,0], res.y[2:4,0])

# Plot result
error_plot(res.y[0], res.y[1], x_exact=x_soln, y_exact=y_soln,
          title='Electron Trajectory', xlabel='x', ylabel='y')
# Error in x,y as a function of time
plt.plot(res.t, res.y[0] - x_soln)
plt.plot(res.t, res.y[1] - y_soln)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20428987921699293, 'fraction_numerical': 0.03519366930445648, 'mean_word_length': 4.489142857142857, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '160091', 'license': 'MIT', 'n_tokens_mistral': 2246, 'n_tokens_neox': 2132, 'n_words': 448, 'path': 'python/Ford1991.ipynb', 'repo_name': 'QTNM/Electron-Tracking', 'size': '4802'}