In [2]:

    

import sys, os
from os.path import join
from mayavi import mlab
import healpy as hp
import matplotlib
#matplotlib.use('Qt4Agg')
matplotlib.interactive(True)
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


    

    




/global/homes/f/fantaye/.conda/envs/work/lib/python2.7/site-packages/IPython/kernel/__init__.py:13: ShimWarning: The `IPython.kernel` package has been deprecated. You should import from ipykernel or jupyter_client instead.
  "You should import from ipykernel or jupyter_client instead.", ShimWarning)

In [ ]:

    

get_ipython().magic('matplotlib inline')

#sys.path
ld_lib=os.getenv('LD_LIBRARY_PATH')
scratch='/global/cscratch1/sd/fantaye' #os.getenv('SCRATCH')
home=os.getenv('HOME')
print('home:',home)
print('scratch',scratch)
romepy=home+"/romepy/"
sys.path.append(romepy)


#---------------------
#my routines
import dotunf as unf

def threed_sphere(m,vmin=None,vmax=None):
    nside = hp.npix2nside(len(m))
    xsize = ysize = 1000

    theta = np.linspace(np.pi, 0, ysize)
    phi   = np.linspace(-np.pi, np.pi, xsize)
    longitude = np.radians(np.linspace(-180, 180, xsize))
    latitude = np.radians(np.linspace(-90, 90, ysize))

    # project the map to a rectangular matrix xsize x ysize
    PHI, THETA = np.meshgrid(phi, theta)
    grid_pix = hp.ang2pix(nside, THETA, PHI)
    grid_map = m[grid_pix]

    # Create a sphere
    r = 0.3
    x = r*np.sin(THETA)*np.cos(PHI)
    y = r*np.sin(THETA)*np.sin(PHI)
    z = r*np.cos(THETA)

    mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
    mlab.clf()

    mlab.mesh(x, y, z, scalars=grid_map, colormap="jet", vmin=vmin, vmax=vmax) 



srcdir=join(scratch,'seedlets/fits4096/SurfDensMap/')
rootdir='../'
fdir=rootdir+'figures/'
if not os.path.exists(fdir):
    os.makedirs(fdir)
    

font = {'family' : 'fantasy',
        #'weight' : 'bold',
        'size'   : 16}

matplotlib.rc('font', **font)

#print os.path.abspath(pview.__file__)

nj=4
nmax=128
lmax=2000
nshell=2*nmax
nside=1024
npix=hp.nside2npix(nside)

#---------------------

#os.listdir(home+'/Radial3Dneedlet/output/seedlets')
fname=join(srcdir,'a2beta/gln2.unf')
cnt=(lmax+1)*nshell*nj
gln=unf.runf(fname,count=cnt, shape=(lmax+1,nshell,nj))

fig = plt.figure()
ax3d = fig.gca(projection='3d')
x = np.linspace(0, lmax, lmax+1)
y = np.linspace(0, nmax, nmax+1)
x,y=np.meshgrid(y,x)
z = gln[:,0:nmax+1,1]
ax3d.plot_surface(x, y, z,  rstride=4, cstride=4,cmap='jet',
    linewidth=0, antialiased=False)



mapj=[]
for j in range(1,2):
	m=np.ndarray((npix,10))
	for ir in range(10):
		fname=join(srcdir,'a2beta/maps/map.unf_gln2_j'+str(j)+'_r'+str(ir)+'.unf')
		#mr=np.fromfile(fname, count=npix,dtype=np.float64)
		mr=unf.runf(fname,count=npix)
		if ir<10: print('j, ir, map.shape,max,min:',j, ir,mr.shape,mr.max(),mr.min())
		m[:,ir]=mr 

	mapj.append(m)
    
mapj[0].shape

#---------------------

#fig, (ax1, ax2)=plt.subplots(ncols=2,figsize=(15,6))
ir=5
hp.mollview(mapj[0][:,ir],title='shell '+str(ir))
ir=8
#hp.mollview(mapj[0][:,ir],title='shell '+str(ir))

threed_sphere(mapj[0][:,ir])


    

In [ ]:

    

%qtconsole