apps120807_mhddata.mhddata_line_meridian

''' A line plot along the specific meridian.

.. image:: ../../../src/scripts/apps120807_mhddata/mhddata_line_meridian_1.2_180.png
'''

import matplotlib.pyplot as plt
import numpy as np

import pyana.pep.mhddata


def main(longitude=0., height=1.2, pp=None):

    r = height
#    longitude = 180.    # Longitude in degrees.

    longitude_r = np.deg2rad(longitude)

    fig = plt.figure()
    ax1 = fig.add_subplot(311)
    ax2 = fig.add_subplot(312)
    ax3 = fig.add_subplot(313)

    latlist = np.arange(-90, 91, 1)   # Every 1 deg.
    latlist_r = np.deg2rad(latlist)

    x = r * np.cos(latlist_r) * np.cos(longitude_r)
    y = r * np.cos(latlist_r) * np.sin(longitude_r)
    z = r * np.sin(latlist_r)

    n = np.zeros_like(latlist_r)
    vx = np.zeros_like(latlist_r)
    vy = np.zeros_like(latlist_r)
    vz = np.zeros_like(latlist_r)
    t = np.zeros_like(latlist_r)
    pth = np.zeros_like(latlist_r)

    if pp == None:
        pp = pyana.pep.mhddata.PlasmaParameter1205()
    
    for i, lat_r in enumerate(latlist_r):
        n[i], vx[i], vy[i], vz[i], t[i], pth[i] = pp.interpolate3d(x[i], y[i], z[i])
        
    vr = (x * vx[i] + y * vy[i] * z * vz[i]) / r

    ax1.plot(latlist, n, label='N')
    ax1.set_ylim(0, n.max() * 1.1)
    ax1.legend(loc='best')
    ax1.set_xlim(-90, 90)
    ax1.set_title('Longitude=%g' % longitude)
    ax1.set_ylabel('N [/cm3]')

    v = np.sqrt(vx * vx + vy * vy + vz * vz)
    vth = pyana.pep.mhddata.t2vth(t, mass=16.0) / 1e3   # in km/s
    ax2.plot(latlist, v, label='|V|')
    ax2.plot(latlist, vth, label='Vth')
    ax2.set_ylim(0, np.max([v.max(), vth.max()]) * 1.1)
    ax2.legend(loc='best')
    ax2.set_xlim(-90, 90)
    ax2.set_ylabel('Vel [km/s]')

    ax3.plot(latlist, vx, label='Vx')
    ax3.plot(latlist, vy, label='Vy')
    ax3.plot(latlist, vz, label='Vz')
    ax3.plot(latlist, vr, label='Vr')
    ax3.set_ylim(np.min([vx.min(), vy.min(), vz.min()]), np.max([vx.max(), vy.max(), vz.max()]))
    ax3.legend(loc='best')
    ax3.set_xlabel('Latitude')
    ax3.set_xlim(-90, 90)
    ax3.set_ylabel('Vel [km/s]')

    fig.savefig('mhddata_line_meridian_%g_%g.png' % (height, longitude))

    return pp

if __name__ == "__main__":
    pp = main(longitude=180)