<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Defining some utilitary functions and importing some modules
Step2:  Query Gaia tables in VizieR
Query table I/337/gaia
Let us query table I/337/gaia 20 arcmin around Messier 44
Step3:  By default, output is limited to 50 rows. Let's change this and redo the query.
Step4:  We add now a constraint to keep only sources with a parallax > 5 mas
Step5:  Send result to VO tools through SAMP
Step6:  Query RRLyrae table and retrieve a light curve
Retrieve list of RRLyrae and sort it by period
Step7:  Retrieve light curve points for Source = 4660711867541110784
Step8:  Plot the light curve
Step9:  Fold the light curve, folded by the period (P1 = 0.21271385 d)
Step11:  Recreate HR diagram from Gaia DR1 paper
Get the data
Step12:  Plot the HR diagram
Step13:  Find Gaia couterparts for a list of AllWISEsources and draw a color-color diagram
Load sample of AllWISE sources in Small Magellanic Cloud
Step14:  Cross-match with Gaia DR1
Step15:  Create color-color plot
Step16:  Retrieve Gaia sources within HST and SDSS coverage
Retrieve SDSS DR9 and HST-R MOCs (coverage)
Step17:  Compute intersection
Step18:  Retrieve Gaia DR1 sources in the coverage of the computed MOC
Step20:  Query SIMBAD to retrieve stars member of a cluster and cross-match with TGAS
Step21:  We now query the xmatch service to find counterparts in TGAS (max search radius
<ASSISTANT_TASK:>
Python Code:
def query_TAP(tap_endpoint, adql_query, table_to_upload=None):
    
    Query a TAP service (designated by its tap_endpoint)
    with a given ADQL query
    
    Query is performed synchronously
    
    Return an AstroPy Table object
    
    import requests
    from astropy.table import Table
    from astropy.io.votable import parse_single_table
    import os
    import tempfile
    import warnings
    
    r = requests.post(tap_endpoint + '/sync', data={'query': adql_query, 'request': 'doQuery', 'lang': 'adql', 'format': 'votable', 'phase': 'run'})
    
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        
    tmp_vot = tempfile.NamedTemporaryFile(delete = False)
    with open(tmp_vot.name, 'w') as h:
        for line in r.iter_lines():
            if line:
                h.write(line.decode(r.encoding)+'\n')
    table = parse_single_table(tmp_vot.name).to_table()
    # finally delete temp files
    os.unlink(tmp_vot.name)
    return table
from astroquery.xmatch import XMatch
import types
import six
from astropy.io import ascii
# monkey patching XMatch
def patched_is_table_available(self, table_id):
        if isinstance(table_id, six.string_types) and (table_id[:7] == 'vizier:'):
            table_id = table_id[7:]
        if not isinstance(table_id, six.string_types):
            return False
        
        return table_id in self.get_available_tables()
    
    
XMatch.is_table_available = types.MethodType(patched_is_table_available, XMatch)
    
from astropy import units as u
    
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib
def plot_scatter_density(xdata, ydata, xlabel, ylabel, title, xlim=None, ylim=None, cmap='viridis', invert_yaxis = True, s=2, grid=False):
    from scipy.stats import gaussian_kde
    import numpy as np
    x = np.reshape(np.array(xdata, copy=False).astype('float'), (len(xdata)))
    y = np.reshape(np.array(ydata, copy=False).astype('float'), len(ydata))
    xy = np.vstack([x,y])
    z = gaussian_kde(xy)(xy)
    
    # Sort the points by density, so that the densest points are plotted last
    idx = z.argsort()
    x, y, z = x[idx], y[idx], z[idx]
    
    
    w = h = 8
    fig, ax = plt.subplots(figsize = (w, h))
    
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    if xlim:
        ax.set_xlim(*xlim)
    if ylim:
        ax.set_ylim(*ylim)
    ax.set_title(title)
    ax.scatter(x, y, c=z, s=s, edgecolor='', cmap=cmap)
    if invert_yaxis:
        ax.invert_yaxis()
    ax.grid(grid)
# disable warnings as not to clutter notebook
import warnings
warnings.simplefilter("ignore")
Explanation: Defining some utilitary functions and importing some modules
End of explanation
from astroquery.vizier import Vizier
result = Vizier.query_region('Messier 44', radius='0d20m0s', catalog='I/337/gaia')
print(result[0])
Explanation: Query Gaia tables in VizieR
Query table I/337/gaia
Let us query table I/337/gaia 20 arcmin around Messier 44
End of explanation
Vizier.ROW_LIMIT = 100000
result = Vizier.query_region('Messier 44', radius='0d20m0s', catalog='I/337/gaia')
print(result[0]['RA_ICRS', 'DE_ICRS', 'Source', 'Plx'])
Explanation: By default, output is limited to 50 rows. Let's change this and redo the query.
End of explanation
result = Vizier(column_filters={"Plx":">5"}).query_region('Messier 44', radius='0d20m0s', catalog='I/337/gaia', )
print(result[0]['RA_ICRS', 'DE_ICRS', 'Source', 'Plx'])
Explanation: We add now a constraint to keep only sources with a parallax > 5 mas
End of explanation
# TODO
Explanation: Send result to VO tools through SAMP
End of explanation
rrlyrae = Vizier.get_catalogs(catalog='I/337/rrlyrae')[0]
rrlyrae.sort('P1')
print(rrlyrae['Source', 'P1', 'RA_ICRS', 'DE_ICRS'])
Explanation: Query RRLyrae table and retrieve a light curve
Retrieve list of RRLyrae and sort it by period
End of explanation
lc_points = Vizier(column_filters={"Source":"4660711867541110784"}).get_catalogs(catalog='I/337/fov')[0]
print(lc_points['ObsTime', 'FG'])
Explanation: Retrieve light curve points for Source = 4660711867541110784
End of explanation
x = lc_points['ObsTime']
y = lc_points['FG']
y_err = lc_points['e_FG']
plt.scatter(x, y, color='r')
plt.errorbar(x, y, yerr=y_err, color='r', linestyle='None')
plt.grid(True)
Explanation: Plot the light curve
End of explanation
period = 0.21271385
#(time+8)%16 - 8
x = (lc_points['ObsTime'] % period) / period
y = lc_points['FG']
y_err = lc_points['e_FG'] 
plt.scatter(x, y, color='r')
plt.errorbar(x, y, yerr=y_err, color='r', linestyle='None')
plt.grid(True)
Explanation: Fold the light curve, folded by the period (P1 = 0.21271385 d)
End of explanation
adql_query = SELECT gaia.ra, gaia.dec, gaia.source_id, gaia.hip, gaia.phot_g_mean_mag+5*log10(gaia.parallax)-10 as g_mag_abs,
hip."B-V"
                FROM "I/337/tgas" as gaia
                inner join "I/311/hip2" as hip
                      on gaia.hip= hip.HIP
                where gaia.parallax/gaia.parallax_error >= 5 and hip."e_B-V" > 0.0 and hip."e_B-V" <= 0.05 and
2.5/log(10)*gaia.phot_g_mean_flux_error/gaia.phot_g_mean_flux <= 0.05
tap_vizier_endpoint = 'http://tapvizier.u-strasbg.fr/TAPVizieR/tap/'
hr_data = query_TAP(tap_vizier_endpoint, adql_query)
print(hr_data)
Explanation: Recreate HR diagram from Gaia DR1 paper
Get the data: query TAP VizieR
End of explanation
plot_scatter_density(hr_data['B-V'], hr_data['g_mag_abs'],
                     title='Gaia DR1 HR diagram', xlabel = 'B-V', ylabel=r'$M_G$', 
                     cmap='viridis', invert_yaxis = True, s=6, grid=True)
Explanation: Plot the HR diagram
End of explanation
from astropy.io import votable
allwise_sources = votable.parse_single_table('allwise-SMC-sample.vot').to_table()
print(len(allwise_sources))
Explanation: Find Gaia couterparts for a list of AllWISEsources and draw a color-color diagram
Load sample of AllWISE sources in Small Magellanic Cloud
End of explanation
allwise_X_gaia = XMatch.query(cat1 = allwise_sources,
                             cat2 = 'vizier:I/337/gaia',
                             max_distance = 1 * u.arcsec,                             
                             colRA1 = 'RAJ2000',
                             colDec1= 'DEJ2000', cache=False)
print(len(allwise_X_gaia))
Explanation: Cross-match with Gaia DR1
End of explanation
plot_scatter_density(allwise_X_gaia['phot_g_mean_mag'] - allwise_X_gaia['W1mag'],
                     allwise_X_gaia['W1mag'] - allwise_X_gaia['W4mag'], title='Color-color plot',
                     xlabel = 'G-W1', ylabel='W1-W4', 
                     cmap='viridis', invert_yaxis = True, s=6, grid=True)
Explanation: Create color-color plot
End of explanation
from mocpy import MOC
moc_sdss = MOC.from_ivorn('ivo://CDS/P/SDSS9/u', nside=512)
moc_hst = MOC.from_ivorn('ivo://CDS/P/HST/R', nside=512)
moc_sdss.plot('SDSS coverage')
moc_hst.degrade_to_order(6).plot('HST coverage')
Explanation: Retrieve Gaia sources within HST and SDSS coverage
Retrieve SDSS DR9 and HST-R MOCs (coverage)
End of explanation
common_coverage = moc_sdss.intersection(moc_hst)
print('Common coverage represents %.5f%% of the sky' % (common_coverage.sky_fraction))
common_coverage.degrade_to_order(6).plot('Common coverage between SDSS and HST')
Explanation: Compute intersection
End of explanation
gaia_sources = common_coverage.query_vizier_table('GAIA-DR1', max_rows=200000)
print(gaia_sources['ra', 'dec'])
Explanation: Retrieve Gaia DR1 sources in the coverage of the computed MOC
End of explanation
stars_in_cluster = query_TAP('http://simbad.u-strasbg.fr/simbad/sim-tap', 
                             SELECT  ra, dec, pmra, pmdec, plx_value, plx_err FROM basic
                                WHERE otype_txt = '*iC'
                                and plx_value>1 and plx_value/plx_err > 5)
print(stars_in_cluster)
Explanation: Query SIMBAD to retrieve stars member of a cluster and cross-match with TGAS
End of explanation
xmatch_result = XMatch.query(cat1 = stars_in_cluster,
                             cat2 = 'vizier:I/337/tgasptyc',
                             max_distance = 1 * u.arcsec,                             
                             colRA1 = 'ra',
                             colDec1= 'dec', cache=False)
print(xmatch_result['ra', 'dec', 'plx_value', 'Plx', 'plx_err', 'e_Plx'])
histvals, binvals, patches = plt.hist(xmatch_result['Plx'], bins=100, facecolor='g', alpha=0.75)
plt.xlabel('Plx')
plt.ylabel('count')
plt.title('Gaia parallax distribution for cluster members')
plt.grid(True)
histvals, binvals, patches = plt.hist(100 * abs(xmatch_result['Plx'] - xmatch_result['plx_value']) / xmatch_result['Plx'], bins=100, facecolor='g', alpha=0.75)
plt.xlabel('Percentage')
plt.ylabel('count')
plt.title('Relative difference between parallax in Simbad and parallax measures by Gaia')
plt.grid(True)
plot_scatter_density(xmatch_result['ra'] - xmatch_result['_RAJ2000'],
                     xmatch_result['dec']  - xmatch_result['_DEJ2000'],
                     r'$\Delta$RA', r'$\Delta$dec', 'Position differences between Simbad and Gaia', xlim=[-5e-5, 5e-5], ylim=[-5e-5, 5e-5],
                     cmap='viridis', invert_yaxis = False, s=6, grid=True)
Explanation: We now query the xmatch service to find counterparts in TGAS (max search radius: 1 arcsec)
End of explanation
<END_TASK>