#Set up astropy and astroquery 
from astropy.table import Table,Column
from astropy.coordinates import SkyCoord, Distance, CartesianRepresentation
from astropy import units as u
from astroquery.vizier import Vizier

#Set up Matplotlib
%config InlineBackend.rc = {}
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

v=Vizier(columns=['L']) #Need to ask for L explicitly as it is not part of the standard returned elements
v.ROW_LIMIT = -1
Cats = v.get_catalogs('VII/189')
pulsarCat=Cats[0]
pulsarCat.keep_columns(["RAJ2000","DEJ2000","dist","P","Age","L"])
#Create Column with dist in pc (not kpc) for WWT
distCol=Column(1000*pulsarCat['dist'].filled(0),name='distance',meta={'ucd': 'pos.distance'}, unit='pc')
pulsarCat.add_column(distCol)
del pulsarCat['dist']
pulsarCat.rename_column('RAJ2000', 'RA')
pulsarCat.rename_column('DEJ2000', 'dec')
pulsarCoords=SkyCoord(pulsarCat['RA'],pulsarCat['dec'],unit=(u.degree, u.degree),\
                   distance=Distance(pulsarCat['distance'],u.pc),frame='icrs')
pulsarGalCoords = pulsarCoords.galactic
pulsarGalCoords.representation='cartesian'
#Now we'll add the projection into the Galactic coordinate, minus signs convert from left-handed towards anti-center to right handed towards center
pulsarCat.add_column(Column(-1.0*pulsarCoords.galactic.represent_as(CartesianRepresentation).x,name='x_gal_coord',meta={'ucd': 'pos.cartesian.x'}),0)
pulsarCat.add_column(Column(-1.0*pulsarCoords.galactic.represent_as(CartesianRepresentation).y,name='y_gal_coord',meta={'ucd': 'pos.cartesian.y'}),1)
pulsarCat.add_column(Column(pulsarCoords.galactic.represent_as(CartesianRepresentation).z,name='z_gal_coord',meta={'ucd': 'pos.cartesian.z'}),2)

pulsarCat

import seaborn
import numpy as np
fig = plt.figure (figsize=(13,3))
ax = fig.add_subplot(121,projection="mollweide")
ax.grid(True)
ax.set_title("Pulsar Sky Distribution")
ax.scatter(pulsarCoords.galactic.l.wrap_at(180.*u.degree).radian,pulsarCoords.galactic.b.radian,c='orange')
a2 = fig.add_subplot(122)
a2.hist(pulsarCat['P'],bins=np.logspace(-4.0, 2.0, 100))
a2.set_title("Pulsar period distribution")
a2.set_xlabel("seconds")
a2.set_ylabel("N")
a2.set_xscale("log")


pulsarCat.filled(-9999).write("pulsars.raw",format='ascii.no_header')