def get_seppa(source_id1, source_id2):
    """Use Gaia DR2 source ids to return the position angle between two sources in degrees and 
    separation in arcsec"""
    import warnings
    warnings.filterwarnings('ignore')
    import numpy as np
    from astroquery.gaia import Gaia
    from myastrotools.gaia_tools import to_polar
    deg_to_mas = 3600000.0
    mas_to_deg = 1.0 / 3600000.0
    job = Gaia.launch_job(
        'SELECT * FROM gaiadr2.gaia_source WHERE source_id = ' + str(
        source_id1))
    j = job.get_results()
    job = Gaia.launch_job(
        'SELECT * FROM gaiadr2.gaia_source WHERE source_id = ' + str(
        source_id2))
    k = job.get_results()
    RAa, RAaerr = j[0]['ra'], j[0]['ra_error'] * mas_to_deg
    RAb, RAberr = k[0]['ra'], k[0]['ra_error'] * mas_to_deg
    Deca, Decaerr = j[0]['dec'], j[0]['dec_error'] * mas_to_deg
    Decb, Decberr = k[0]['dec'], k[0]['dec_error'] * mas_to_deg
    raa_array = np.random.normal(RAa, RAaerr, 10000)
    rab_array = np.random.normal(RAb, RAberr, 10000)
    deca_array = np.random.normal(Deca, Decaerr, 10000)
    decb_array = np.random.normal(Decb, Decberr, 10000)
    rho_array, pa_array = to_polar(raa_array, rab_array, deca_array, decb_array
        )
    pa, paerr = np.mean(pa_array).value, np.std(pa_array).value
    rho, rhoerr = np.mean(rho_array).value, np.std(rho_array).value
    return rho / 1000, rhoerr / 1000, pa, paerr
