In [ ]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import matplotlib.colors as colors
import nugridpy.utils as utils
import sys
import logging
import collections
import os
# you can use sys.path.insert to add your own cloned copy of PyPPM which you can modify
# sys.path.insert(0,'/home/user/home/PyPPM/')
from ppmpy import ppm
# plotting utilities
ifig = 0
ptrack = {}
cb = utils.linestylecb # colours
# create smaller labelsizes on axes
plt.rc('xtick', labelsize=7)
plt.rc('ytick', labelsize=7)
# named tuple for using rprofs and momsdata
hydro = collections.namedtuple('hydro', ['moms','rprof'])
# turn off matplotlib messages
logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
In [ ]:
# load in data, use LowZRAWD N16
moms_dir = '/data/ASDR/PPMstar/LowZRAWD/N16-LowZRAWD-1536-10x-burn-moms/myavsbq'
rprof_dir = '/data/ASDR/PPMstar/LowZRAWD/N16-LowZRAWD-1536-10x-burn-moms/prfs/'
# the list of variables that are stored within the momsdata cube.
var_list = ['xc','ux','uy','uz','|ut|','|ur|','|w|','P','rho','fv']
# the dump to plot. read this in immediately
mydump = 650
# create the "hydro" tuple
n16 = hydro(ppm.MomsDataSet(moms_dir,mydump,2,var_list,rprofset=ppm.RprofSet(rprof_dir)),
ppm.RprofSet(rprof_dir))
Some constants for converting $X_{\mathrm{cld}}$ fluid to $X_{\mathrm{H}}$
In [ ]:
atomicnocld = n16.rprof.get('atomicnocld',fname=0)
atomicnoH = 1.
fkcld = n16.rprof.get('fkcld',fname=0)
airmu = n16.rprof.get('airmu',fname=0)
cldmu = n16.rprof.get('cldmu',fname=0)
We want to find a radius in which the H is just beginning to burn which in the LowZRAWD runs is roughly in the center of the convection zone
In [ ]:
# this ensures that I can open an close the figure that was plotted within this cell.
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 0
# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
else:
ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
# set the colour/linestyle counter to 0 and create axes
cbc = 0
ax = fig.add_subplot(111)
# a good radius to use
myradius = 16
# get the Xcld and plot it. We can access the "rprof" intuitively with the hydro tuple
Xcld = n16.rprof.compute_Xcld(mydump)
r = n16.rprof.get('R',fname=mydump,resolution='l')
ax.plot(r, Xcld, color=cb(cbc)[2], ls=cb(cbc)[0])
cbc += 1
# create the line referring to a "suitable" radius
myylim = ax.get_ylim()
ax.vlines(myradius, *myylim, color='k', linestyle=':')
# add limits and labels
ax.set_xlabel('R / Mm')
ax.set_ylabel(r'$X_{\mathrm{cld}}$')
ax.set_yscale('log')
ax.set_ylim([1e-8,1.1])
ax.set_xlim([r.min(),r.max()])
In [ ]:
# number of points for interpolation onto a sphere for plotting and calculations
# I calculate how many points this is at a given radius and what mode l we can resolve
npoints = n16.moms.sphericalHarmonics_lmax(myradius)[-1]
# gridline width for mollweide
gridlines = 0.1
First get all of the appropriate moms data and then interpolate the data onto a sphere for plotting in the mollweide projection. Methods of the MomsDataSet class will handle the appropriate coordinates and computing the interpolations
In [ ]:
# start off with ur
ux = n16.moms.get('ux',fname=mydump)
uy = n16.moms.get('uy',fname=mydump)
uz = n16.moms.get('uz',fname=mydump)
ur, utheta, uphi = n16.moms.get_spherical_components(ux, uy, uz)
ur_r, theta_r, phi_r = n16.moms.get_spherical_interpolation(ur, myradius, npoints=npoints,
plot_mollweide=True)
# convert to km/s
ur_r *= 1e3
In [ ]:
# to get XH we first need Xcld which is derived from fv
fv = n16.moms.get('fv',fname=mydump)
Xcld = fv/((1. - fv)*(airmu/cldmu) + fv)
XH = atomicnoH * (fkcld / atomicnocld) * Xcld
# we use the same number of points so theta_r and phi_r apply to XH_r as well!
XH_r = n16.moms.get_spherical_interpolation(XH, myradius, npoints=npoints)
We will showcase two different types of interpolation methods for $\rho$, the linear and moments interpolation types
In [ ]:
# get rho on the sphere with linear interpolation
rho = n16.moms.get('rho',fname=mydump)
rho_trilinear_r = n16.moms.get_spherical_interpolation(rho, myradius, npoints=npoints, method='trilinear')
rho_moments_r = n16.moms.get_spherical_interpolation(rho, myradius, npoints=npoints, method='moments')
# estimate spherical average of rho on myradius
avg_rho_trilinear = rho_trilinear_r.mean()
avg_rho_moments = rho_moments_r.mean()
# compute perturbation on the sphere
rho_pert_trilinear_r = (rho_trilinear_r - avg_rho_trilinear) / avg_rho_trilinear
rho_pert_moments_r = (rho_moments_r - avg_rho_moments) / avg_rho_moments
We have the $\theta$ and $\phi$ coordinates which are defined from $\pi/2, -\pi/2$ for $\theta$ and $-\pi, \pi$ for $\phi$ instead of spherical coordinates in order to plot the data with matplotlib. We can simply plot with a scatter plot keeping in mind that we want "pixel" sized markers so that it will show most/all of the data without any overlap. With the 384$^{3}$ data there are enough points at the radii within the convection zone so that there are no blank spots in the plots
In [ ]:
# this ensures that I can open an close the figure that was plotted within this cell.
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 1
# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
else:
ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
# create mollweide projection axes. Leave room for a colourbar
ax = plt.axes([0.01, 0.01, 0.98, 0.88], projection='mollweide')
# ensure even bounds
ur_bounds = np.max(np.abs(ur_r))
# plot as a scatter plot. Make the markersize, s, a pixel (look at the definition of s)
ax.scatter(phi_r, theta_r, s=(72./fig.dpi)**2, marker=',', c=ur_r,
cmap='RdBu_r',vmin=-ur_bounds,vmax=ur_bounds)
# create a colorbar for the scatter plot data, use a new axes for it
cax = fig.add_axes([0.01, 0.03, 0.98, 0.03])
cbar1 = fig.colorbar(ax.collections[0], cax=cax, orientation='horizontal',format='%.1f',
label=r'v$_{\mathrm{r}}$ / km s$^{-1}$')
# remove ticklabels and thin out grid
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(color='k', linewidth=gridlines)
This is an older technique that uses the $\theta$ and $\phi$ points and then creates roughly equal area trianges that are filled in with colour. Because it needs to create these triangles, it is a much slower plotting technique and it can only perform reasonably with $\approx 10,000$ points
In [ ]:
# get less points
npoints_plot = 10000
ur_r_triang, theta_r_triang, phi_r_triang = n16.moms.get_spherical_interpolation(ur, myradius,
npoints=npoints_plot, plot_mollweide=True)
# convert to km/s
ur_r_triang *= 1e3
In [ ]:
# this ensures that I can open an close the figure that was plotted within this cell.
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 2
# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
else:
ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
ax = plt.axes([0.01, 0.01, 0.98, 0.88], projection='mollweide')
# create mollweide projection axes. Leave room for a colourbar
triang = tri.Triangulation(phi_r_triang,theta_r_triang)
# ensure even bounds
ur_bounds = np.max(np.abs(ur_r))
# tripcolor plots colour values of unordered data points of the triangles. This is why it takes so long
ax.tripcolor(triang,ur_r_triang,cmap='RdBu_r',vmin=-ur_bounds,vmax=ur_bounds)
# create a colorbar for the ax tripcolor, use a new axes for it
cax = fig.add_axes([0.01, 0.03, 0.98, 0.03])
cbar1 = fig.colorbar(ax.collections[0], cax=cax, orientation='horizontal',format='%.1f',
label=r'v$_{\mathrm{r}}$ / km s$^{-1}$')
# remove ticklabels and thin out grid
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(color='k', linewidth=gridlines)
In [ ]:
# this ensures that I can open an close the figure that was plotted within this cell.
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 3
# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
else:
ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
# create mollweide projection axes. Leave room for a colourbar
ax = plt.axes([0.01, 0.01, 0.98, 0.88], projection='mollweide')
# plot with particular bounds. Anything below 1e-7 is 0
XH_log = [1e-7, 1e-4]
XH_r[np.where(XH_r < 1e-7)] = 0.
# plot as a scatter plot. Make the markersize, s, a pixel (look at the definition of s)
ax.scatter(phi_r, theta_r, s=(72./fig.dpi)**2, marker=',', c=XH_r,
cmap='Spectral_r',norm=colors.LogNorm(*XH_log))
# create a colorbar for the scatter plot data, use a new axes for it
cax = fig.add_axes([0.01, 0.03, 0.98, 0.03])
cbar1 = fig.colorbar(ax.collections[0], cax=cax, orientation='horizontal',format='%.1e',
label=r'X$_{\mathrm{H}}$')
# remove ticklabels and thin out grid
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(color='k', linewidth=gridlines)
# remove ticklabels and thin out grid
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(color='k', linewidth=gridlines)
In [ ]:
# this ensures that I can open an close the figure that was plotted within this cell.
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 4
# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
else:
ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
# create mollweide projection axes. Leave room for a colourbar
ax = plt.axes([0.01, 0.01, 0.98, 0.88], projection='mollweide')
# ensure even bounds
rho_bounds = np.abs(rho_pert_trilinear_r).max()
# plot as a scatter plot. Make the markersize, s, a pixel (look at the definition of s)
ax.scatter(phi_r, theta_r, s=(72./fig.dpi)**2, marker=',', c=rho_pert_trilinear_r,
cmap='BrBG',vmin=-rho_bounds,vmax=rho_bounds)
# create a colorbar for the scatter plot data, use a new axes for it
cax = fig.add_axes([0.01, 0.03, 0.98, 0.03])
cbar1 = fig.colorbar(ax.collections[0], cax=cax, orientation='horizontal',format='%.0e',
label=r'$\left( \rho - \langle \rho \rangle \right) / \langle \rho \rangle$')
# remove ticklabels and thin out grid
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(color='k', linewidth=gridlines)
Because linear interpolation is used on a quantity that is very close to spherically symmetric there are systematic ringing effects seen. Instead, for a quantity like $\rho$ it is more appropriate to use the moments type interpolation as it can be much "smoother" IF the underlying data is nearly continuous. This is not true for quantities like $u_{\mathrm{r}}$ or X$_{\mathrm{H}}$
In [ ]:
# this ensures that I can open an close the figure that was plotted within this cell.
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 5
# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
else:
ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
ppm.add_plot(celln,ifig,ptrack)
# create mollweide projection axes. Leave room for a colourbar
ax = plt.axes([0.01, 0.01, 0.98, 0.88], projection='mollweide')
# ensure even bounds
rho_bounds = np.abs(rho_pert_moments_r).max()
# plot as a scatter plot. Make the markersize, s, a pixel (look at the definition of s)
ax.scatter(phi_r, theta_r, s=(72./fig.dpi)**2, marker=',', c=rho_pert_moments_r,
cmap='BrBG',vmin=-rho_bounds,vmax=rho_bounds)
# create a colorbar for the scatter plot data, use a new axes for it
cax = fig.add_axes([0.01, 0.03, 0.98, 0.03])
cbar1 = fig.colorbar(ax.collections[0], cax=cax, orientation='horizontal',format='%.0e',
label=r'$\left( \rho - \langle \rho \rangle \right) / \langle \rho \rangle$')
# remove ticklabels and thin out grid
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.grid(color='k', linewidth=gridlines)
In [ ]: