# apps130704_negativeion.app02_trace_flatmoonΒΆ

First, to get the image, flat moon is assumed.

The solar wind velocity with -400 km/s in x direction. The moon is in y-z plane, a flat plane. The b field is changeable in the simulation. The default is +y direction, resulting the electric field +z.

The backscattering may occur +300 km/s initial velocity, but also changeable.

Still single particle.

''' First, to get the image, flat moon is assumed.

The solar wind velocity with -400 km/s in x direction.
The moon is in y-z plane, a flat plane.
The b field is changeable in the simulation.
The default is +y direction, resulting the electric field +z.

The backscattering may occur +300 km/s initial velocity, but
also changeable.

Still single particle.
'''

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.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

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.001  # seconds.

return bsw, vsw, m, q, dt

def tracing_main(bsw, vsw, m, q, dt, p0, v0):
esw = np.cross(bsw, vsw)   # E = -VxB

t = 0.
p = p0
v = v0

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) <= 5:  # 10 sec ~ 4000 km
p = LeapFlog.progress_position(p, vh, dt)
t += dt
plist.append(p)
tlist.append(t)

if p[0] > 200e3 or abs(p[1]) > 2000e3 or abs(p[2]) > 2000e3:
break

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

return tlist, plist, thlist, vhlist

def tracing_test3():
'''Tracing test 3 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

p0 = np.array([0., 0, 0])
#    v = v0 = np.array([0., 0, 0])
v0 = np.array([300., 0, 0]) * 1e3

tlist, plist, thlist, vhlist = tracing_main(bsw, vsw, m, q, dt, p0, v0)
tlist2, plist2, thlist2, vhlist2 = tracing_main(bsw, vsw, m, -q, dt, p0, v0)

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)
action="store_true", dest="verbose")
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_test3()

if __name__ == "__main__":
main()