apps130704_negativeion.app01_trace_protonΒΆ
As a start, I will trace a single particle under the solar wind.
This is to confirm the way of using tracing package and its validity.
''' As a start, I will trace a single particle under the solar wind.
This is to confirm the way of using tracing package and its validity.
'''
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 plot_xy(xpos, ypos, ax=None, xlabel='X [km]', ylabel='Y [km]', **kwds):
''' X_Y slice (or X-Z, Y-Z)
'''
if ax == None:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xpos, ypos, **kwds)
# ax.set_aspect('equal')
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
return ax
def plot_xyz(xpos, ypos, zpos, ax=None, xlabel='X', ylabel='Y', zlabel='Z', **kwds):
''' XYZ representation
ax should be instanced as
``ax = fig.gca(projection='3d')``
or equivalent. If None given to ax, new figure
and axes is generated.
'''
if ax == None:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(xpos, ypos, zpos, **kwds)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_zlabel(zlabel)
return ax
def add_sphere(ax, r=1738, **kwds):
''' Add sphere.
rstride=4 and cstride=4 is recommended.
'''
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = r * np.outer(np.cos(u), np.sin(v))
y = r * np.outer(np.sin(u), np.sin(v))
z = r * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, **kwds)
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.01 # seconds.
return bsw, vsw, m, q, dt
def tracing_test2():
'''Tracing test 1 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) = (0, 0, 0)
X(0) = (0, 0, 0)
**Expected result**
Proton moves along the so-called pick up ions in x-z plane.
'''
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
esw = np.cross(bsw, vsw) # E = -VxB
p = p0 = np.array([0., 0, 0])
# v = v0 = np.array([0., 0, 0])
v = v0 = np.array([300., 0, 0]) * 1e3
t = .0
force = lambda pp, vv: (q / m * (esw + np.cross(vv, bsw)))
### 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]
while abs(t) <= 20: # 10 sec ~ 4000 km
p = LeapFlog.progress_position(p, vh, dt)
t += dt
plist.append(p)
tlist.append(t)
vh = LeapFlog.progress_velocity(vh, force(p, vh), dt)
tv += dt
vhlist.append(vh)
thlist.append(tv)
plist = np.array(plist) * 1e-3 # (N, 3) with unit of km
vhlist = np.array(vhlist) * 1e-3 # (N+1, 3) with unit of km/s
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]')
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]')
ax = plot_xyz(plist[:, 0], plist[:, 1], plist[:, 2])
ax.plot([plist[0, 0]], [plist[0, 1]], [plist[0, 2]], marker='o')
ax = plot_xyz(vhlist[:, 0], vhlist[:, 1], vhlist[:, 2], xlabel='Vx', ylabel='Vy', zlabel='Vz')
ax.plot([vhlist[0, 0]], [vhlist[0, 1]], [vhlist[0, 2]], marker='o')
def tracing_test1():
'''Tracing test 1 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) = (-400, 0, 0)
X(0) = (0, 0, 0)
**Expected result**
Proton moves with velocity -400 km/s in the x direction without any change.
'''
bsw, vsw, m, q, dt = traceing_setting()
p = p0 = np.array([0, 0, 0])
v = v0 = np.array([-400e3, 0, 0])
t = 0
force = lambda pp, vv: (q / m * (esw + np.cross(vv, bsw)))
### 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]
while t <= 10: # 10 sec ~ 4000 km
p = LeapFlog.progress_position(p, vh, dt)
t += dt
plist.append(p)
tlist.append(t)
vh = LeapFlog.progress_velocity(vh, force(p, vh), dt)
tv += dt
vhlist.append(vh)
thlist.append(tv)
plist = np.array(plist) * 1e-3 # (N, 3) with unit of km
vhlist = np.array(vhlist) * 1e-3 # (N+1, 3) with unit of km/s
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]')
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]')
ax = plot_xyz(plist[:, 0], plist[:, 1], plist[:, 2])
ax.plot([plist[0, 0]], [plist[0, 1]], [plist[0, 2]], marker='o')
ax = plot_xyz(vhlist[:, 0], vhlist[:, 1], vhlist[:, 2], xlabel='Vx', ylabel='Vy', zlabel='Vz')
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_test2()
if __name__ == "__main__":
main()