%matplotlib inline
from __future__ import division

import numpy as np
import healpy as hp
import astroML

#This removes some nasty deprecation warnings that do not interfere with the execution
import warnings
warnings.filterwarnings('ignore')

# This IPython magic generates a table with version information
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, healpy, astroML, astroPy

for NSIDE in 2.**np.arange(6):
    print 'The number of pixels for NSIDE = %2d is %5d' %(NSIDE, hp.nside2npix(NSIDE))

NSIDE = 1
ipix = np.arange(hp.nside2npix(NSIDE))
print 'The numbering of the pixels for NSIDE = %d is: \n' %NSIDE, ipix

theta, phi = hp.pix2ang(NSIDE,ipix)
print '\n array theta = ', theta/np.pi, 'radians'
print '\n array phi = ', phi/np.pi, 'radians'

NSIDE = 1
hp.ang2pix(NSIDE, np.pi/2, np.pi)

# If arrays of the same length: the matching is element to element
theta = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4])
phi = np.array([0, np.pi/2, np.pi, 3*np.pi/2])
hp.ang2pix(NSIDE, theta, phi)

# Case of arrays of different lengths: we get an exception
theta = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4])
phi2 = np.array([0, np.pi/2, np.pi, 3*np.pi/2, 3.5*np.pi/2])
try:
    hp.ang2pix(NSIDE, theta, phi2)
except Exception as e:
    print "Atttention: ",e

# But if we change the theta array to a column array
# Then the numpy broadcasting is possible,
# and returns a 2D array of pixel numbers
theta2 = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4]).reshape(-1,1)
phi2 = np.array([0, np.pi/2, np.pi, 3*np.pi/2, 3.5*np.pi/2])
hp.ang2pix(NSIDE, theta2, phi2)

npixels = hp.nside2npix(32)
np.random.seed(0)
test = np.random.rand(npixels)
hp.mollview(test, min=0, max=1,
            title = 'Mollweide projection test',
            unit='This is a description of the measurement unit')

NSIDE = 2
npixels = hp.nside2npix(NSIDE)
img = np.linspace(0, 255, num=npixels)
index = np.arange(npixels)
theta, phi= hp.pix2ang(NSIDE,index)
hp.mollview(img, min=0, max = 255, unit='Range of values')
hp.projscatter(theta, phi)
for i in index:
    hp.projtext(theta[i]-0.05, phi[i], i)

lon = [30, 105, 270]
lat = [45, -30, 60]
hp.mollview(title="Mollweide in (lon, lat) coordinates")
hp.graticule()    # add grid
hp.projscatter(lon, lat, lonlat=True)
hp.projtext(30, 45,'(30,45)', lonlat=True)
hp.projtext(105, -30,'(105,-30)', lonlat=True)
hp.projtext(270, 60,'(270,60)', lonlat=True);

from astroML.datasets import fetch_wmap_temperatures
wmap_unmasked = fetch_wmap_temperatures(masked=False)
wmap_unmasked.shape

hp.nside2npix(512)

hp.mollview(wmap_unmasked, min=-1, max=1, 
            title='Unmasked WMAP data \n in galactic coordinates',
            unit=r'$\Delta$T (mK)', xsize = 200)

hp.mollview(wmap_unmasked, min=-1, max=1, 
            title='Unmasked WMAP data \n in equatorial coordinates',
            coord='GC', rot=(180,0),
            unit=r'$\Delta$T (mK)', xsize = 200)
hp.graticule()

hp.projtext(180, 2,'RA=180', lonlat=True)
hp.projtext(90, 2, 'RA=90', lonlat=True)
hp.projtext(270, 2, 'RA=270', lonlat=True);