apps130704_negativeion.app03_trace_moon

Then, approaching to real image, circular moon is assumed.

The solar wind velocity with -400 km/s in x direction. The moon is in y-z plane, a 1738 km sphere. 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.

''' Then, approaching to real image, circular moon is assumed.

The solar wind velocity with -400 km/s in x direction.
The moon is in y-z plane, a 1738 km sphere.
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 = 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.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]

    inner2 = 1738e3 ** 2   # At the lunar surface
    outer2 = (1738e3 * 3) ** 2   # At the 3Rm

    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) * 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_test4():
    '''Tracing test 4 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([1738.1e3, 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)
    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_test4()


if __name__ == "__main__":
    main()