#!/usr/bin/env python
"""Analyze ALMA census: observed-vs-unobserved by class; best UC HII targets."""
import numpy as np
from astropy.table import Table
LOG=open('analyze_census.log','w')
def p(*a):
    s=' '.join(str(x) for x in a);print(s);print(s,file=LOG);LOG.flush()

t=Table.read('alma_census.csv',format='csv')
p('Census rows:',len(t))
nobs=np.array(t['n_obs'],float)
obs=nobs>0
rc=np.array([str(x) for x in t['radio_class']])

p('\n=== Observed vs unobserved (all 800 RMS HII regions) ===')
p(f'Observed by ALMA (>=1 archive obs within 30"): {int(obs.sum())} / {len(t)} ({100*obs.mean():.0f}%)')
p(f'Not observed: {int((~obs).sum())}')

p('\n=== By radio_class ===')
for k in ['UCHII','compact','classical','none']:
    sel=rc==k; n=int(sel.sum()); no=int((sel&obs).sum())
    lbl={'UCHII':'UC HII (<=0.1pc)','compact':'compact(0.1-0.5)','classical':'classical(>0.5)','none':'no CORNISH cover'}[k]
    if n: p(f'  {lbl:20s}: {no:3d}/{n:3d} observed ({100*no/n:.0f}%)')

# UC HII specific
uc=t[rc=='UCHII']
ucobs=np.array(uc['n_obs'],float)>0
p(f'\n=== UC HII candidates (strict <=0.1pc): {len(uc)} ===')
p(f'  observed: {int(ucobs.sum())}/{len(uc)}')
# resolution adequacy: need to resolve <=0.1pc at distance d -> ang=0.1/d/0.00485 arcsec
br=np.array(uc['best_res_as'],float)
dist=np.array(uc['Dist'],float)
ang_uc = 0.1/(dist*0.0048481)   # angular size of 0.1pc
resolves = np.isfinite(br) & (br < ang_uc)   # beam finer than the UC HII size
p(f'  with beam finer than the 0.1pc scale (resolves structure): {int((resolves).sum())}')
hires = np.isfinite(br) & (br<0.5)
p(f'  with sub-0.5" data: {int(hires.sum())}')
hasline=np.array(uc["has_line"],float)>0
p(f'  with spectral-line data (<10 km/s): {int((hasline&ucobs).sum())}')

# best targets table
good = uc[(np.array(uc['n_obs'],float)>0)]
good.sort('best_res_as')
cols=['Name','GLON','Dist','Dpc','n_obs','bands','best_res_as','max_mrs_as','has_line','min_sens_10kms','proj_ids']
p('\n=== Top 25 best-observed UC HII (finest resolution) ===')
gg=good[cols][:25]
for r in gg:
    p(f"  {str(r['Name']):22s} d={float(r['Dist']):4.1f} D={float(r['Dpc']):.3f}pc "
      f"res={float(r['best_res_as']):.3f}\" MRS={float(r['max_mrs_as']):5.1f}\" "
      f"b={str(r['bands']):8s} line={int(r['has_line'])} n={int(float(r['n_obs']))} {str(r['proj_ids'])[:30]}")
good.write('uchii_observed_best.fits',overwrite=True)
p(f'\nWrote uchii_observed_best.fits ({len(good)} observed UC HII)')
LOG.close()