''' Plot MHD data in 2D in the noon-midnight plane
.. image:: ../../../src/scripts/apps120807_mhddata/mhddata_meridian_500km_1.png
.. image:: ../../../src/scripts/apps120807_mhddata/mhddata_meridian_500km_2.png
.. image:: ../../../src/scripts/apps120807_mhddata/mhddata_meridian_500km_3.png
'''
import bisect
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Arrow
from matplotlib.ticker import LogLocator
from pyana.pep import mhddata
def main(pp=None, mf=None, r=1.2):
if pp == None:
pp = mhddata.PlasmaParameter1205()
if mf == None:
mf = mhddata.MagneticField1201()
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)
fig3 = plt.figure()
ax3 = fig3.add_subplot(111)
cntr = ax.contourf(xc, zc, n_y0.T, 8)
clb = fig.colorbar(cntr)
clb.set_label('Density [/cm3]')
vth_y0 = mhddata.t2vth(t_y0, mass=16.) / 1e3
cntr = ax3.contourf(xc, zc, np.log10(vth_y0.T), np.linspace(0.5, 3.5, 16))
clb = fig3.colorbar(cntr)
clb.set_label('LOG Vth [km/s]')
ax.set_aspect('equal')
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-1.5, 1.5)
ax3.set_aspect('equal')
from .mhddata_bsimple_trace import main as trace
trace(mf=mf, ax=ax3)
# Circle
lowbnd = Circle((0, 0), 1.1, facecolor='gray', edgecolor='gray')
gany = Circle((0, 0), 1., facecolor='w', edgecolor='k')
ax.add_patch(lowbnd)
ax.add_patch(gany)
lowbnd = Circle((0, 0), 1.1, facecolor='gray', edgecolor='gray')
gany = Circle((0, 0), 1., facecolor='w', edgecolor='k')
ax3.add_patch(lowbnd)
ax3.add_patch(gany)
# Line plot for context.
fig2 = plt.figure(figsize=(8, 14))
ax21 = fig2.add_subplot(511)
ax22 = fig2.add_subplot(512)
ax23 = fig2.add_subplot(513)
ax24 = fig2.add_subplot(514)
ax25 = fig2.add_subplot(515)
latsdeg = np.arange(0, 360, 5)
ns = []
vxs = []
vys = []
vzs = []
ts = []
pths = []
vrs = []
for latdeg in latsdeg:
latrad = np.deg2rad(latdeg)
x = np.cos(latrad) * r
y = 0
z = np.sin(latrad) * r
n, vx, vy, vz, t, pth = pp.interpolate3d(x, y, z)
ns.append(n)
vxs.append(vx)
vys.append(vy)
vzs.append(vz)
ts.append(t)
pths.append(pth)
vrs.append((vx * x + vz * z) / r)
ax21.plot(latsdeg, ns)
ax21.set_ylabel('N')
ax22.plot(latsdeg, vxs, label='Vx')
ax22.plot(latsdeg, vys, label='Vy')
ax22.plot(latsdeg, vzs, label='Vx')
ax22.plot(latsdeg, vrs, lw=3, label='Vr')
ax22.grid()
ax22.set_ylabel('Vx,Vy,Vz,Vr')
ax22.legend(loc='best')
ax23.plot(latsdeg, np.array(vrs) * np.array(ns) * 1e5, label='F=n*Vr') # 1e5 comes from the unit to be /cm2 s
ax23.grid()
ax23.axhline(0, c='k')
ax23.set_ylabel('F=N x Vr')
vths = mhddata.t2vth(np.array(ts), mass=16.) / 1000.
ax24.plot(latsdeg, vths, label='Vth')
vs = np.sqrt(np.array(vxs) ** 2 + np.array(vys) ** 2 + np.array(vzs) ** 2)
ax24.plot(latsdeg, vs, label='|V|')
ax24.set_ylabel('Vth,|V|')
ax24.legend(loc='best')
ax25.plot(latsdeg, pths)
ax25.set_ylabel('Pth')
scale = 3e-2 # 100 km/s corresponds to 0.1 Rg
for latdeg in latsdeg:
latrad = np.deg2rad(latdeg)
x = np.cos(latrad) * r
y = 0
z = np.sin(latrad) * r
n, vx, vy, vz, t, pth = pp.interpolate3d(x, y, z)
print(x, z, vx, vz)
# arr = Arrow(x, z, vx * scale, vz * scale)
# ax.add_patch(arr)
ax.arrow(x, z, vx * scale, vz * scale, shape='full', length_includes_head=True, head_width=0.05, facecolor='k')
ax3.arrow(x, z, vx * scale, vz * scale, shape='full', length_includes_head=True, head_width=0.05, facecolor='k')
### I tried to plot thermal velocity as a circle from the edge of the bulk velocity, but the presentation is not intuitive enough.
#vth = mhddata.t2vth(t, mass=16) / 1e3 # in km/s
#ax.add_patch(Circle((x + vx * scale, z + vz * scale), vth * scale, facecolor='none', edgecolor='gray'))
ax.arrow(1.0, -1.43, 10 * scale, 0 * scale, shape='full', length_includes_head=True, head_width=0.05, facecolor='k')
ax.text(1.0, -1.43, '10 km/s', horizontalalignment='right')
ax3.arrow(1.0, -1.43, 10 * scale, 0 * scale, shape='full', length_includes_head=True, head_width=0.05, facecolor='k')
ax3.text(1.0, -1.43, '10 km/s', horizontalalignment='right')
ax3.set_xlim(-1.5, 1.5)
ax3.set_ylim(-1.5, 1.5)
fig.savefig('mhddata_meridian_500km_1.png')
fig2.savefig('mhddata_meridian_500km_2.png')
fig3.savefig('mhddata_meridian_500km_3.png')
return pp
if __name__ == "__main__":
retval = main()