%matplotlib inline
%config InlineBackend.figure_formats=['svg']

import matplotlib
matplotlib.rcParams['figure.figsize'] = (9, 6)

import numpy as np
from matplotlib import pyplot as plt

from astropy import units as u
from astropy import coordinates as coord
from astropy.table import Column
from astroquery.vizier import Vizier

number_names = np.array(['zero', 'one', 'two', 'three', 'four', 'five'])
print(number_names[[3, 0, 0, 1, 2, 5, 5, 4]])

Vizier.ROW_LIMIT = -1  # If we don't do that, Vizier only send 50 rows.
region_center = coord.SkyCoord(34.8*u.deg, -4.4*u.deg)
sdss9, tycho2 = Vizier.query_region(region_center, width="40m", catalog=['sdss9', 'tycho2'])

plt.plot(sdss9['RAJ2000'], sdss9['DEJ2000'], "bo", 
         markeredgecolor='None', alpha=.5, label="SDSS9 galaxies")
plt.plot(tycho2['RAmdeg'], tycho2['DEmdeg'], "y*", 
         label="Tycho2 stars", markersize=14)
plt.xlabel("RA")
plt.ylabel("Dec")
plt.legend()

print("The SDSS9 data was observed in {}.".format(np.int(np.mean(sdss9['ObsDate']))))

# Correcting the star position with their proper motions.
# Note: we convert the columns to arrays because of a bug
star_ra = (np.array(tycho2['RAmdeg']) + np.array(tycho2['pmRA']) /
           1000. / 3600. / np.cos(np.deg2rad(np.array(tycho2['DEmdeg'])))
           * (2008-2000))
star_dec = (np.array(tycho2['DEmdeg']) + np.array(tycho2['pmDE']) / 
            1000. / 3600. * (2008-2000))

# Star positions.
star_pos = coord.SkyCoord(star_ra * u.degree, star_dec * u.degree)

# Star r magnitudes
star_rmag = tycho2['VTmag'] - 0.01 - 0.57545 * (tycho2['BTmag'] - tycho2['VTmag'])

# Mask radius associated to each star
star_rmag.unit = None
star_mask_rad = np.exp(-0.43 * star_rmag + 8.4) * u.arcsec

# The SDSS9 catalogue positions
galaxy_pos = coord.SkyCoord(sdss9['RAJ2000'], sdss9['DEJ2000'])

coord.SkyCoord.search_around_sky?

idx_galaxy, idx_star, d2d, _ = star_pos.search_around_sky(galaxy_pos, np.max(star_mask_rad))

to_flag_idx_galaxy = 

sdss9.add_column(Column(data=np.zeros(len(sdss9)), name="mask_tycho2", dtype=np.int))
sdss9['mask_tycho2'][to_flag_idx_galaxy] = 1

plt.plot(sdss9['RAJ2000'], sdss9['DEJ2000'], "bo", 
         markeredgecolor='None', alpha=.5, label="SDSS9 galaxies")
plt.plot(sdss9['RAJ2000'][sdss9['mask_tycho2'] == 1], 
         sdss9['DEJ2000'][sdss9['mask_tycho2'] == 1], 
         "go", markeredgecolor='None', alpha=.5,
         label="Flagged galaxies")
plt.plot(tycho2['RAmdeg'], tycho2['DEmdeg'], "y*", 
         label="Tycho2 stars", markersize=14)
plt.xlabel("RA")
plt.ylabel("Dec")
plt.legend()