apps120807_mhddata.mhddata_meridian

''' A simple sample to plot MHD data in 2D

Plot of the velocity vector above density.

.. image:: ../../../src/scripts/apps120807_mhddata/mhddata_meridian.png


'''

import bisect

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Arrow

from pyana.pep import mhddata

def main(pp=None):

    if pp == None:
        pp = mhddata.PlasmaParameter1205()

    x = pp.xlist()
    y = pp.ylist()
    z = pp.zlist()

    n = pp.nlist()
    lnn = np.log10(n)
    lnn[np.isinf(lnn)] = np.nan

    vx = pp.vxlist()
    vy = pp.vylist()
    vz = pp.vzlist()

    t = pp.tlist()

    # 1-D coordinate value
    xc = x[:, 0, 0]
    yc = y[0, :, 0]
    zc = z[0, 0, :]

    # y=0 index
    y0idx = bisect.bisect(yc, 0) - 1

    # Slice
    x_y0 = x[:, y0idx, :]   # 200, 200
    y_y0 = y[:, y0idx, :]   # 200, 200
    z_y0 = z[:, y0idx, :]   # 200, 200
    n_y0 = n[:, y0idx, :]
    lnn_y0 = lnn[:, y0idx, :]
    vx_y0 = vx[:, y0idx, :]
    vy_y0 = vy[:, y0idx, :]
    vz_y0 = vz[:, y0idx, :]
    t_y0 = t[:, y0idx, :]

    # Plot density/velocity vector
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.contourf(xc, zc, lnn_y0.T, 8)

    tmax= t_y0.max() * 0.2
    print(tmax)

    scale = 1e-2  # 100 km/s corresponds to 0.1 Rg
    nstep = 2     # every nstep data points.
    for xx, yy, uu, vv, tt in zip(x_y0[::nstep, ::nstep].flatten(), z_y0[::nstep, ::nstep].flatten(), vx_y0[::nstep, ::nstep].flatten(), vz_y0[::nstep, ::nstep].flatten(), t_y0[::nstep, ::nstep].flatten()):
        if xx < 1.3 and xx > -1.3 and yy < 1.3 and yy > -1.3:
            arr = Arrow(xx, yy, uu * scale, vv * scale, width=np.clip(tt / tmax, 0.1, 1))
            ax.add_patch(arr)

    ax.set_aspect('equal')
    ax.set_xlim(-1.3, 1.3)
    ax.set_ylim(-1.3, 1.3)

    # Circle
    gany = Circle((0, 0), 1., facecolor='w', edgecolor='k')
    ax.add_patch(gany)

    for rcirc in (1.1, 1.2, 1.3):
        ax.add_patch(Circle((0, 0), rcirc, facecolor='none', edgecolor='gray'))

    fig.savefig('mhddata_meridian.png')

    return pp

if __name__ == "__main__":
    retval = main()