import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from optparse import OptionParser
from irfpy.util.leapflog import LeapFlog
def traceing_setting():
bsw = np.array([0, 10e-9, 0])
vsw = np.array([-400e3, 0, 0])
m = 1.67e-27
q = 1.60e-19
dt = 0.001 # seconds.
return bsw, vsw, m, q, dt
class Tracer:
def __init__(self, bsw, vsw, m, q, dt):
self.bsw = bsw
self.vsw = vsw
self.esw = np.cross(bsw, vsw) # E = -vxB
self.m = m
self.q = q
self.dt = dt
self.force = lambda pp, vv: (self.q / self.m *
(self.esw + np.cross(vv, self.bsw)))
def run(self, p0, v0):
''' Run the tracing from p0 and v0
'''
bsw = self.bsw
vsw = self.vsw
esw = self.esw
force = self.force
dt = self.dt
p = p0
v = v0
t = 0
### Half step forward for velocity
vh = LeapFlog.progress_velocity(v0, force(p0, v0), dt / 2.)
tv = dt / 2.
plist = [p, ] # List for position at t=0, dt, 2dt, ...
vhlist = [v, vh] # List for velocity at t=0, 0.5dt, 1.5dt, ...
tlist = [0, ]
thlist = [0, tv]
inner2 = 1738e3 ** 2 # At the lunar surface
outer2 = 1838e3 ** 2 # At 100 km
while abs(t) <= 5: # 10 sec ~ 4000 km
p = LeapFlog.progress_position(p, vh, dt)
t += dt
plist.append(p)
tlist.append(t)
r2 = (p * p).sum()
if r2 > outer2 or r2 < inner2:
break
vh = LeapFlog.progress_velocity(vh, force(p, vh), dt)
tv += dt
vhlist.append(vh)
thlist.append(tv)
plist = np.array(plist) # (N, 3) with unit of km
vhlist = np.array(vhlist) # (N+1, 3) with unit of km/s
return tlist, plist, thlist, vhlist
def tracing_test5():
'''Tracing test 5 is under the following condition.
**Fields**
B-field: Bsw=(0, 10, 0) nT
V-sw: Vsw=(-400, 0, 0) km/s
E-field: Esw=(0, 0, 4) mV/m
**Particle to trace**
Species: proton
V(0) = (300, 0, 0)
X(0) = (0, 0, 0)
'''
bsw, vsw, m, q, dt = traceing_setting()
# dt = -dt
# q = -q
# bsw = np.array((5, 5, 0)) * 1e-9
# bsw = np.array((0, 8, 8)) * 1e-9
tracer_p = Tracer(bsw, vsw, m, q, dt)
tracer_n = Tracer(bsw, vsw, m, -q, dt)
p0 = np.array([1738.1e3, 0, 0])
# v = v0 = np.array([0., 0, 0])
v0 = np.array([300., 0, 0]) * 1e3
tlist, plist, thlist, vhlist = tracer_p.run(p0, v0)
tlist2, plist2, thlist2, vhlist2 = tracer_n.run(p0, v0)
from .app03_trace_moon import plot_xy, plot_xyz
fig = plt.figure(figsize=(6, 10))
for i in range(3):
ax = fig.add_subplot(6, 1, i + 1)
plot_xy(tlist, plist[:, i], ax=ax, xlabel='t', ylabel='pos [km]')
plot_xy(tlist2, plist2[:, i], ax=ax, xlabel='t', ylabel='pos [km]')
for i in range(3):
ax = fig.add_subplot(6, 1, i + 4)
plot_xy(thlist, vhlist[:, i], ax=ax, xlabel='t', ylabel='vel [km/s]')
plot_xy(thlist2, vhlist2[:, i], ax=ax, xlabel='t', ylabel='vel [km/s]')
ax = plot_xyz(plist[:, 0], plist[:, 1], plist[:, 2])
plot_xyz(plist2[:, 0], plist2[:, 1], plist2[:, 2], ax=ax)
ax.plot([plist[0, 0]], [plist[0, 1]], [plist[0, 2]], marker='o')
# ax.set_xlim(-100, 100)
# ax.set_ylim(-100, 100)
# ax.set_zlim(-100, 100)
ax = plot_xyz(vhlist[:, 0], vhlist[:, 1], vhlist[:, 2], xlabel='Vx', ylabel='Vy', zlabel='Vz')
plot_xyz(vhlist2[:, 0], vhlist2[:, 1], vhlist2[:, 2], xlabel='Vx', ylabel='Vy', zlabel='Vz', ax=ax)
ax.plot([vhlist[0, 0]], [vhlist[0, 1]], [vhlist[0, 2]], marker='o')
def main():
usage = "usage: %prog [options] arg"
parser = OptionParser(usage)
parser.add_option("-f", "--file", dest="filename",
help="read data from FILENAME")
parser.add_option("-v", "--verbose",
action="store_true", dest="verbose")
parser.add_option("-q", "--quiet",
action="store_false", dest="verbose")
(options, args) = parser.parse_args()
# if len(args) != 1:
# parser.error("incorrect number of arguments")
# if options.verbose:
# print "reading %s..." % options.filename
tracing_test5()
if __name__ == "__main__":
main()