Healpy Coordinate Transformation Residuals

coord_transform_test.py

import healpy as hp
import numpy as np
import matplotlib
matplotlib.rc('text', usetex=True)
import matplotlib.pyplot as plt


euler_transf_index = {
    'CG': 1,
    'GC': 2,
    'CE': 3,
    'EC': 4,
    'EG': 5,
    'GE': 6
}

def transform_coords_euler(lon, lat, coord_in, coord_out):
    '''
    Convert between astronomical coordinate systems.
    
    All inputs and outputs are in degrees.
    
    Valid values for coord_in and coord_out are
      'G' for 'Galactic'                (l, b)
      'C' for 'Celestial'/'Equatorial'  (ra, dec)
      'E' for 'Ecliptic'                (lon, lat)
    '''
    
    return hp.rotator.euler(lon, lat, euler_transf_index[coord_in + coord_out])


def transform_coords_rotator(lon, lat, coord_in, coord_out):
    '''
    Convert between astronomical coordinate systems.
    
    All inputs and outputs are in degrees.
    
    Valid values for coord_in and coord_out are
      'G' for 'Galactic'                (l, b)
      'C' for 'Celestial'/'Equatorial'  (ra, dec)
      'E' for 'Ecliptic'                (lon, lat)
    '''
    
    rot = hp.rotator.Rotator(coord=[coord_in, coord_out])
    
    t_in = np.radians(90. - lat)
    p_in = np.radians(lon)
    
    t_out, p_out = rot(t_in, p_in)
    lon_out = np.degrees(p_out)
    lat_out = 90. - np.degrees(t_out)
    
    return lon_out, lat_out


def rand_lon_lat(n):
    '''
    Generate n random longitudes and latitudes, drawn uniformly from the
    surface of a sphere.
    '''
    u = np.random.random(n)
    v = np.random.random(n)
    
    lon = 360. * u
    lat = 90. - np.degrees(np.arccos(2. * v - 1.))
    
    return lon, lat


def dtheta_hist(n):
    l,b = rand_lon_lat(n)
    
    ra0, dec0 = transform_coords_euler(l, b, 'G', 'C')
    ra1, dec1 = transform_coords_rotator(l, b, 'G', 'C')
    
    dir0 = np.vstack([ra0, dec0])
    dir1 = np.vstack([ra1, dec1])
    
    d = hp.rotator.angdist(dir0, dir1, lonlat=True)
    
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.hist(d*3600000., bins=25, normed=True)
    ax.set_xlabel(r'$\Delta \theta \ \left( \mathrm{mas} \right)$', fontsize=18)
    ax.set_ylabel(r'$\mathrm{frequency}$', fontsize=18)
    ax.set_title(r'$\mathrm{Error \ Histogram}$', fontsize=22)
    ax.set_ylim(0., ax.get_ylim()[1]*1.1)


def dtheta_map(nside):
    npix = hp.pixelfunc.nside2npix(nside)
    t,p = hp.pixelfunc.pix2ang(nside, np.arange(npix))
    
    # Galactic -> Celestial
    l = np.degrees(p)
    b = 90. - np.degrees(t)
    
    ra0, dec0 = transform_coords_euler(l, b, 'G', 'C')
    ra1, dec1 = transform_coords_rotator(l, b, 'G', 'C')
    
    dir0 = np.vstack([ra0, dec0])
    dir1 = np.vstack([ra1, dec1])
    
    d = hp.rotator.angdist(dir0, dir1, lonlat=True)
    
    hp.visufunc.mollview(3600000.*d, coord='G', unit='mas', title=r'$\Delta \theta$', format='%.2g')
    
    # Celestial -> Galactic
    ra = np.degrees(p)
    dec = 90. - np.degrees(t)
    
    l0, b0 = transform_coords_euler(ra, dec, 'C', 'G')
    l1, b1 = transform_coords_rotator(ra, dec, 'C', 'G')
    
    dir0 = np.vstack([l0, b0])
    dir1 = np.vstack([l1, b1])
    
    d = hp.rotator.angdist(dir0, dir1, lonlat=True)
    
    hp.visufunc.mollview(3600000.*d, coord='C', unit='mas', title=r'$\Delta \theta$', format='%.2g')


def main():
    dtheta_hist(500000)
    dtheta_map(512)
    plt.show()
    
    return 0

if __name__ == '__main__':
    main()