apps130704_negativeion.app04_1_initializer

Initial condition maker for batch calculation

''' Initial condition maker for batch calculation
'''
import time

import numpy as np
import matplotlib.pyplot as plt

rm = 1738.
rm2 = rm * rm

def position_creator2(maxsize=1000):
    ''' Initial position at the Moon surface

    :keyword maxsize: Maximum size to create.
        On average, the returned size is
        pi / 4 ~ 78% of the maxsize.
    :returns: Array with (n, 3) shape,
        where n is not more than maxsize.

    This is faster, compared with position_creator,
    particularly when maxsize is bigger.
    The execution time will not be affected much
    by increasing the maxsize.
    '''
    y = np.random.uniform(-rm, rm, maxsize)
    z = np.random.uniform(-rm, rm, maxsize)

    d2 = y * y + z * z
    masked_array = d2 > rm2

    y = np.ma.masked_array(y, mask=masked_array).compressed()
    z = np.ma.masked_array(z, mask=masked_array).compressed()
    x = np.sqrt(rm2 - np.ma.masked_array(d2, mask=masked_array).compressed())

    if x.size != y.size or y.size != z.size:
        raise RuntimeError('Something wrong in the algorithm!')

    return np.array([x, y, z])
        

def position_creator():
    ''' Initial position at the Moon surface
    '''

    y = np.random.uniform(-rm, rm)
    z = np.random.uniform(-rm, rm)

    d2 = y * y + z * z
    if d2 <= rm2:
        x = np.sqrt(rm2 - d2)
        return x, y, z
    else:
        return position_creator()

def direction_creator_rnd2(maxsize=1000):
    ''' Initial velocity vector.

    :keyword maxsize: Maximum size to create.
        On average, the returned size is
        pi / 6 ~ 52% of the maxsize.
    :returns: Array with (n, 3) shape,
        where n is not more than maxsize.

    This is faster, compared with position_creator,
    particularly when maxsize is bigger.
    The execution time will not be affected much
    by increasing the maxsize.
    '''
    x = np.random.uniform(-1, 1, size=maxsize)
    y = np.random.uniform(-1, 1, size=maxsize)
    z = np.random.uniform(-1, 1, size=maxsize)

    d2 = x * x + y * y + z * z
    masked_array = d2 > 1

    x = np.ma.masked_array(x, mask=masked_array).compressed()
    y = np.ma.masked_array(y, mask=masked_array).compressed()
    z = np.ma.masked_array(z, mask=masked_array).compressed()
    d = np.ma.masked_array(np.sqrt(d2), mask=masked_array).compressed()

    if x.size != y.size or y.size != z.size:
        raise RuntimeError('Something wrong in the algorithm!')

    return np.array([x / d, y / d, z / d])



def direction_creator_rnd():
    ''' Create the random direction
    '''
    x = np.random.uniform(-1, 1)
    y = np.random.uniform(-1, 1)
    z = np.random.uniform(-1, 1)

    d2 = x * x + y * y + z * z
    if d2 > 1:
        return direction_creator_rnd()
    else:
        d = np.sqrt(d2)
        return x / d, y / d, z / d


def test_position_creator():
    n = 1000    # 1000 particle examined.

    t0 = time.time()
    
#    xyzlist = []
#    for i in range(n):
#        xyzlist.append(position_creator())
#    xyzlist = np.array(xyzlist).T   # 3x10000

    xyzlist = position_creator2(maxsize=n)

    t1 = time.time()
    print('Creating %d test_position: %f sec' % (xyzlist.shape[1], t1 - t0))

    plt.figure()

    plt.subplot(221)
    plt.scatter(xyzlist[0, :], xyzlist[1, :])
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.gca().set_aspect('equal')

    plt.subplot(222)
    plt.scatter(xyzlist[0, :], xyzlist[2, :])
    plt.xlabel('X')
    plt.ylabel('Z')
    plt.gca().set_aspect('equal')

    plt.subplot(223)
    plt.scatter(xyzlist[1, :], xyzlist[2, :])
    plt.xlabel('Y')
    plt.ylabel('Z')
    plt.gca().set_aspect('equal')

    plt.subplot(224)
#    plt.scatter(xyzlist[0, :], np.sqrt(xyzlist[1, :] ** 2 + xyzlist[2, :] ** 2))
#    plt.xlabel('X')
#    plt.ylabel('sqrt(Y^2+Z^2)')
#    plt.gca().set_aspect('equal')
    lon = np.rad2deg(np.arctan2(xyzlist[1, :], xyzlist[0, :]))
    lat = np.rad2deg(np.arcsin(xyzlist[2, :] / rm))

    plt.scatter(lon, lat)
    plt.xlabel('Long.')
    plt.ylabel('Lat.')
    plt.xlim(-180, 180)
    plt.ylim(-90, 90)
    plt.gca().set_aspect('equal')

def test_direction_creator_rnd():
    n = 2000    # 2000 particle examined.

    t0 = time.time()
    
#    vxyzlist = []
#    for i in range(n):
#        vxyzlist.append(direction_creator_rnd())
#    vxyzlist = np.array(vxyzlist).T   # 3x10000

    vxyzlist = direction_creator_rnd2(maxsize=n)

    t1 = time.time()
    print('Creating %d test_position: %f sec' % (vxyzlist.shape[1], t1 - t0))

    plt.figure()

    plt.subplot(221)
    plt.scatter(vxyzlist[0, :], vxyzlist[1, :])
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.gca().set_aspect('equal')

    plt.subplot(222)
    plt.scatter(vxyzlist[0, :], vxyzlist[2, :])
    plt.xlabel('X')
    plt.ylabel('Z')
    plt.gca().set_aspect('equal')

    plt.subplot(223)
    plt.scatter(vxyzlist[1, :], vxyzlist[2, :])
    plt.xlabel('Y')
    plt.ylabel('Z')
    plt.gca().set_aspect('equal')

    plt.subplot(224)
#    plt.scatter(vxyzlist[0, :], np.sqrt(vxyzlist[1, :] ** 2 + vxyzlist[2, :] ** 2))
#    plt.xlabel('X')
#    plt.ylabel('sqrt(Y^2+Z^2)')
#    plt.gca().set_aspect('equal')
    lon = np.rad2deg(np.arctan2(vxyzlist[1, :], vxyzlist[0, :]))
    lat = np.rad2deg(np.arcsin(vxyzlist[2, :]))

    plt.scatter(lon, lat)
    plt.xlabel('Long.')
    plt.ylabel('Lat.')
    plt.xlim(-180, 180)
    plt.ylim(-90, 90)
    plt.gca().set_aspect('equal')

def main():
    '''Main script'''

    test_position_creator()
    test_direction_creator_rnd()

if __name__ == "__main__":
    main()