apps120807_mhddata.mhddata_b_meridian

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

.. image:: ../../../src/scripts/apps120807_mhddata/mhddata_b_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(mf=None):

    if mf == None:
        mf = mhddata.MagneticField1201()

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

    # 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

    # Plot density/velocity vector
    fig = plt.figure()
    ax = fig.add_subplot(111)

#    idx = (np.random.random(4000) * 40000).astype(np.int)

    for xx, zz in zip(x_y0.flatten(), z_y0.flatten()):
        if xx > -3 and xx < 3 and zz > -3 and zz < 3:
            try:
                bx, by, bz = mf.interpolate3d(xx, 0, zz)
                print(bx, by, bz)
            except IndexError:
                continue

            ax.arrow(xx, zz, bx / 1000, bz / 1000)
        
    # 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'))

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

    fig.savefig('mhddata_b_meridian.png')

    return mf

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