#!/usr/bin/env python
"""Phase 1: build master catalog from RMS survey (Lumsden+2013) reproducing
the Padoan & Gieles 2026 embedded-source sample, then classify."""
import numpy as np
from astropy.table import Table
from astroquery.vizier import Vizier

def fetch_rms():
    v = Vizier(row_limit=-1)
    t = v.get_catalogs('J/ApJS/208/11')[0]
    return t

def build():
    try:
        t = Table.read('rms_table1_raw.fits')
    except FileNotFoundError:
        t = fetch_rms(); t.write('rms_table1_raw.fits', overwrite=True)

    typ = np.array([str(x) for x in t['Type']])
    young = np.isin(typ, ['HII region', 'YSO', 'HII/YSO'])
    lbol = np.ma.filled(t['Lbol'].astype(float), np.nan)
    dist = np.ma.filled(t['Dist'].astype(float), np.nan)
    logL = np.full(len(t), np.nan)
    pos = np.isfinite(lbol) & (lbol > 0)
    logL[pos] = np.log10(lbol[pos])

    # Padoan & Gieles cuts: positive Lbol, d<=18 kpc, 2.9 <= logL <= 6.6
    sel = (young & pos & np.isfinite(dist) & (dist <= 18)
           & (logL >= 2.9) & (logL <= 6.6))
    m = t[sel].copy()
    m['logLbol'] = logL[sel]
    m['simple_type'] = typ[sel]
    m.write('master_phase1.fits', overwrite=True)
    print(f'Master sample: {len(m)} sources')
    for ty in ['HII region', 'YSO', 'HII/YSO']:
        print(f'  {ty}: {(typ[sel]==ty).sum()}')
    return m

if __name__ == '__main__':
    build()