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

# spherical harmonics tools
import pyshtools.expand 
import pyshtools.spectralanalysis

# 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

# for mppnp runs of N16, there is an associated mass coordinates that define the convection zone
# these are used to use radii within the convection zone only
mtop = 2.991384128949e+04
mbot = 1.149874300454e+04

# create the "hydro" tuple
n16 = hydro(ppm.MomsDataSet(moms_dir,mydump,2,var_list,rprofset=ppm.RprofSet(rprof_dir)),
              ppm.RprofSet(rprof_dir))


    

In [ ]:

    

# get the radius and mass rprof so that we can interpolate to the mtop and mbot of the convection zone
r = n16.rprof.get('R',fname=0,resolution='l')
m = n16.rprof.compute_m(fname=mydump)

# interpolate
rbot = ppm.interpolate(m,r,mbot)
rtop = ppm.interpolate(m,r,mtop)

# spectrogram radii with a spacing of 2*deex (half a moms cell)
step = 2*n16.rprof.get('deex',fname=0)
radii = np.linspace(rbot, rtop, int((rtop - rbot)/step))

# to store the power per l at each and plot with imshow, there must be a constant sized array. For most 
# radii I cannot resolve a large l so those coefficients will be zero
max_l = n16.moms.sphericalHarmonics_lmax(radii.max())[0]

# arrays to hold lmax at a radius and the power at r as a function of l
power_ur_l_r = np.zeros((max_l, len(radii)))
lmax_r = np.zeros(len(radii))


    

In [ ]:

    

# grab the radial velocity on the grid
ur = n16.moms.get_spherical_components('ux','uy','uz',fname=mydump)[0]

# for every radii in the convection zone
for i,radius in enumerate(radii):
    
    # what is the maximum l mode that I can resolve at my radius?
    mylmax, N, npoints = n16.moms.sphericalHarmonics_lmax(radius)
    lmax_r[i] = mylmax
    
    # get the interpolated values of ur on a sphere. These must be sampled a certain way in theta and phi
    # so that pyshtools can efficiently calculate the spherical harmonics (uses a fourier transform)
    ur_interp = n16.moms.sphericalHarmonics_format(ur, radius, lmax=mylmax)
        
    # get coefficients and compute the power, drop l=0, it is just the mean!
    coeffs = pyshtools.expand.SHExpandDH(ur_interp, sampling=2)
    power_ur_l_r[0:mylmax,i] = pyshtools.spectralanalysis.spectrum(coeffs, unit='per_l')[1:]


    

In [ ]:

    

# celln
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=(6,3.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
else:
    ifig += 1; fig = plt.figure(ifig,figsize=(6,3.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)

ax = fig.add_subplot(111)
plt.rc('xtick', labelsize=7)
plt.rc('ytick', labelsize=7)
    
# scale the power by the maximum within that array
scaled_power = power_ur_l_r/ np.max(power_ur_l_r,axis=0)[np.newaxis, :]

# normalize logarithmically
vmax = 1.0
vmin = 1e-4
lognorm = colors.LogNorm(vmin=vmin, vmax=vmax, clip=True)

# imshow of the scaled power
ax.imshow(scaled_power,interpolation='none',origin='lower',norm=lognorm,aspect='auto',
         extent=(radii[0]-0.5*np.diff(radii)[0],radii[-1] + 0.5*np.diff(radii)[0],0.5,0.5+max_l))

# create a colorbar for the axs1 tripcolor, use a new axes for it
cbar1 = fig.colorbar(ax.images[0], ax=ax, orientation='vertical', 
             label=r'S$_{cc}(\ell, r)$ / Max(S$_{cc}(r)$)',aspect=25)

# plot details
ax.set_xlabel('r / Mm')
ax.set_ylabel(r'$\ell$')
fig.tight_layout()


    

In [ ]: