Currently RAPT has a limited selection of field models. However, it is easy to add new field models to it. In this notebook we describe how to add analytic fields (i.e., fields that have an explicit formula in terms of t, x, y, and z).
All new fields are subclasses of the _Field
class. By subclassing, the new field class inherits all the auxiliary field-related functions, so we don't have to write them again. All we need to provide is the new E
and B
. This is called overriding in OOP parlance.
This model is already built into the RAPT package, but it provides a good example for creating new fields. The double-dipole model consists of two identical dipoles, both moments in the negative-z direction, separated by 20$R_e$. This arrangement generates a planar magnetopause between them, and mimics the compression of the dayside magnetosphere by the solar wind.
%matplotlib inline
import pylab as pl
import rapt
import numpy as np
class DoubleDipole(rapt.fields._Field):
def __init__(self, B0=rapt.B0, distance=20*rapt.Re, imagestrength=1):
rapt.fields._Field.__init__(self)
self.gradientstepsize = rapt.Re/1000
self._dd = distance # distance between two dipoles
assert imagestrength >= 1
self._k = imagestrength # >=1. Relative strength of the image dipole
self._coeff = -B0*rapt.Re**3
def B(self, tpos):
t,x,y,z = tpos
B1 = np.array([3*x*z, 3*y*z, (2*z*z -x*x- y*y)]) / pow(x*x+y*y+z*z, 5.0/2.0)
x -= self._dd
B2 = self._k * np.array([3*x*z, 3*y*z, (2*z*z -x*x- y*y)]) / pow(x*x+y*y+z*z, 5.0/2.0)
return self._coeff*(B1+B2)
As an entirely new model, consider a magnetic dipole coincident with an electric charge. To construct a new class, we override both E
and B
.
class ChargedDipole(rapt.fields._Field):
def __init__(self, B0=1, Q=1):
rapt.fields._Field.__init__(self)
self.B0=B0 # field strength (T) at unit distance on the x-y plane
self.Q=Q # charge (C)
self._k = 8.9875517873681764e9 # Coulomb constant
self.static = False # needs to be set to false if E and/or dB/dt is nonzero.
def B(self, tpos):
t,x,y,z = tpos
return self.B0 * np.array([3*x*z, 3*y*z, (2*z*z -x*x- y*y)]) / pow(x*x+y*y+z*z, 5.0/2.0)
def E(self, tpos):
t,x,y,z = tpos
return self._k*self.Q * np.array([x,y,z]) / pow(x*x+y*y+z*z, 3.0/2.0)
Try it out with a small charge.
f = ChargedDipole(Q=1e-6)
Follow a proton's trajectory
p = rapt.Particle([5,0,0], [0,1,0], t0=0, mass=rapt.m_pr, charge=rapt.e, field=f)
p.setke(1) # set kinetic energy to 1 eV
Check the cyclotron period and cyclotron radius to gain an idea about the scale of the motion.
p.cycper()
8.5710888990731629e-06
p.cycrad()
0.098787549981589062
Follow the proton for some time.
p.advance(4e-4)
Display the trajectory on the x-y plane.
pl.plot(p.getx(), p.gety())
pl.xlabel("x")
pl.ylabel("y");
<matplotlib.text.Text at 0x7f5f6f52a630>