#!/usr/bin/env python
"""Crossmatch UC HII (and all RMS HII) to the 15 ALMA-IMF protocluster fields."""
import numpy as np
from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u
LOG=open('almaimf_match.log','w')
def p(*a):
    s=' '.join(str(x) for x in a);print(s);print(s,file=LOG);LOG.flush()

# ALMA-IMF 15 fields (Motte+2022): name, GLON, GLAT (deg), approx field radius
IMF=[('W43-MM1',30.817,-0.057),('W43-MM2',29.961,-0.022),('W43-MM3',29.963,-0.020),
     ('W51-E',49.490,-0.388),('W51-IRS2',49.490,-0.369),('G008.67',8.671,-0.358),
     ('G010.62',10.624,-0.384),('G012.80',12.805,-0.200),('G327.29',327.293,-0.580),
     ('G328.25',328.251,-0.532),('G333.60',333.604,-0.213),('G337.92',337.922,-0.456),
     ('G338.93',338.926,0.394),('G351.77',351.776,-0.538),('G353.41',353.409,-0.360)]
inames=[x[0] for x in IMF]
ic=SkyCoord([x[1] for x in IMF]*u.deg,[x[2] for x in IMF]*u.deg,frame='galactic')

for fn,lab,hf in [('uchii_candidates.fits','UC HII',True),('master_phase1_classified.fits','all HII',False)]:
    t=Table.read(fn)
    if 'simple_type' in t.colnames and not hf:
        t=t[np.array([str(x) for x in t['simple_type']])=='HII region']
    gl=np.array(t['GLON'],float); gb=np.array(t['GLAT'],float)
    c=SkyCoord(gl*u.deg,gb*u.deg,frame='galactic')
    idx,d2d,_=c.match_to_catalog_sky(ic)
    # ALMA-IMF 12m field of view ~ 1' per pointing, mosaics up to ~2-3'. use 90"
    m=d2d<90*u.arcsec
    p(f'\n{lab}: in an ALMA-IMF field (<90"): {int(m.sum())}')
    for i in np.where(m)[0]:
        p(f"  {str(t['Name'][i]):22s} l={gl[i]:7.3f} b={gb[i]:+.3f} -> {inames[idx[i]]:10s} ({d2d.arcsec[i]:.0f}\")")
LOG.close()