#!/usr/bin/env python
"""Phase 6: crossmatch UC HII (and all HII) to Wells+2024 accretion-flow sample."""
import numpy as np
from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u
LOG=open('phase6b.log','w')
def p(*a):
    s=' '.join(str(x) for x in a);print(s);print(s,file=LOG);LOG.flush()

w=Table.read('wells2024_tabled1.fits')   # 87 sources w/ accretion-flow modelling
p('Wells cols',w.colnames)
from collections import Counter
p('Wells Class:',dict(Counter([str(x) for x in w['Class']]).most_common()))
cw=SkyCoord(np.array(w['_RA.icrs'],float)*u.deg,np.array(w['_DE.icrs'],float)*u.deg)

def sky(t,ra,de,unit):
    r=t[ra];d=t[de]
    if r.dtype.kind in ('S','U'):
        return SkyCoord([x.decode() if isinstance(x,bytes) else str(x) for x in r],
                        [x.decode() if isinstance(x,bytes) else str(x) for x in d],unit=unit)
    return SkyCoord(np.array(r,float)*u.deg,np.array(d,float)*u.deg)

for fn,lab in [('uchii_candidates.fits','UC HII (102)'),('master_phase1_classified.fits','all HII (800)')]:
    t=Table.read(fn)
    if 'simple_type' in t.colnames:
        t=t[np.array([str(x) for x in t['simple_type']])=='HII region']
    c=sky(t,'RAJ2000','DEJ2000',(u.hourangle,u.deg))
    idx,d2d,_=c.match_to_catalog_sky(cw)
    m=d2d<10*u.arcsec
    p(f'\n{lab}: matched to Wells accretion-flow sample (<10"): {int(m.sum())}')
    if m.sum():
        for i in np.where(m)[0]:
            j=idx[i]
            p(f"  {str(t['Name'][i]):22s} -> Wells {str(w['ALMAGAL'][j])} Class={str(w['Class'][j])} "
              f"Mass={float(w['Mass'][j]):.0f} Lum={float(w['Lum'][j]):.0f}")
LOG.close()