#!/usr/bin/env python
"""Build the H-alpha surface-brightness cache for the OB sample (resumable).
SkyView 'H-Alpha Comp' (Finkbeiner) per star -> _xm_halpha_ob.ecsv. Slow
(~6800 queries); run in background. circ_mass_ob.py folds it in when present.
"""
import os
import numpy as np
from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u
from astroquery.skyview import SkyView
from circ_mass_ob_common import load_ob

ra, dec, gl, gb, d, logL, isO = load_ob()
N = len(d)
print(f'OB sample: N={N}')
CACHE = '_xm_halpha_ob.ecsv'
I_R = {}
if os.path.exists(CACHE):
    for r in Table.read(CACHE):
        I_R[int(r['id'])] = float(r['I_R'])
for i in range(N):
    if i in I_R:
        continue
    try:
        imgs = SkyView.get_images(position=SkyCoord(ra[i], dec[i], unit=u.deg),
                                  survey=['H-Alpha Comp'], pixels=20,
                                  radius=0.06 * u.deg, cache=False)
        I_R[i] = float(np.nanmedian(imgs[0][0].data))
    except (ConnectionError, TimeoutError, ValueError, IndexError, OSError):
        I_R[i] = np.nan
    if (i + 1) % 100 == 0:
        Table(dict(id=list(I_R), I_R=[I_R[k] for k in I_R])).write(CACHE, overwrite=True)
        print(f'  {i+1}/{N}')
Table(dict(id=list(I_R), I_R=[I_R[k] for k in I_R])).write(CACHE, overwrite=True)
print(f'done: {sum(np.isfinite(v) for v in I_R.values())}/{N} measured')