apps120807_mhddata.mhddata_isoheight

''' A simple sample to plot radial flux at 1.2 Rm.

.. image:: ../../../src/scripts/apps120807_mhddata/mhddata_isoheight.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()

    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, :]

    r = 1.2  # At which altitude?

    lonlist = np.arange(0, 361, 2) * np.pi / 180.
    latlist = np.arange(-90, 90, 2) * np.pi / 180.

    lonmesh, latmesh = np.meshgrid(lonlist, latlist)
    datamesh = np.zeros_like(lonmesh)
    
    
    for ilat, lat in enumerate(latlist):
        for ilon, lon in enumerate(lonlist):
            x0 = r * np.cos(lat) * np.cos(lon)
            y0 = r * np.cos(lat) * np.sin(lon)
            z0 = r * np.sin(lat)

            ix, iy, iz = pp.nearest_neighbor(x0, y0, z0)

            xr = xc[ix]
            yr = yc[iy]
            zr = zc[iz]
#            print lon * 180. / np.pi, lat * 180. / np.pi, '/', x0, xr, '/', y0, yr, '/', z0, zr
            distancer = np.sqrt(xr * xr + yr * yr + zr * zr)

            nr = n[ix, iy, iz]
            vxr = vx[ix, iy, iz]
            vyr = vy[ix, iy, iz]
            vzr = vz[ix, iy, iz]

            vradial = (xr * vxr + yr * vyr + zr * vzr ) / distancer
#            vradial = np.sqrt(vxr * vxr + vyr * vyr + vzr * vzr )
            vradial *= 1e5   # km/s to be cm/s

#            datamesh[ilat, ilon] = nr  # Sorry, first density.
#            datamesh[ilat, ilon] = vradial  # Sorry, second radial velocity.
            datamesh[ilat, ilon] = -nr * vradial  # Then, multiple

#            datamesh[ilat, ilon] = xr
#            datamesh[ilat, ilon] = yr
#            datamesh[ilat, ilon] = zr

    datamesh = np.ma.masked_less(datamesh, 0)
    datamesh = np.ma.log10(datamesh)
    vmin = 7

    plt.figure()
    pc = plt.pcolor(np.rad2deg(lonmesh), np.rad2deg(latmesh), datamesh, vmin=vmin)
    plt.colorbar(pc)

    plt.savefig('mhddata_isoheight.png')

    return pp

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