#!/usr/bin/env python
"""Independent off-plane cross-check with AKARI/FIS Bright Source Catalogue
(II/298/fisbsc; 65/90/140/160 um) — colder/FIR, different instrument than IRAS.
Same distance-independent gate: F_bol(65-160um) >= F_bol(1e4 Lsun,3kpc).
Cross-match the IRAS off-plane candidates to see if AKARI independently recovers
them, and whether AKARI flags any off-plane luminous source IRAS missed."""
import numpy as np
from astroquery.vizier import Vizier
from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u

Lsun = 3.828e26; kpc = 3.0856776e19
F_THRESH = 1e4 * Lsun / (4 * np.pi * (3 * kpc) ** 2)
lam = np.array([65., 90., 140., 160.]) * 1e-6
nu = 2.99792458e8 / lam

def fbol(s65, s90, s140, s160):
    Fnu = np.vstack([s65, s90, s140, s160]).T * 1e-26
    return np.trapezoid(Fnu[:, ::-1], nu[::-1], axis=1)

Vizier.ROW_LIMIT = -1
v = Vizier(columns=['objName', 'RAJ2000', 'DEJ2000',
                    'S65', 'S90', 'S140', 'S160', 'q_S65', 'q_S90'], row_limit=-1)
# FIS PSC: pull bright-90um all-sky, filter |b|>5 locally
r = v.query_constraints(catalog='II/298/fis', S90='>30')
t = r[0]
print('AKARI FIS S90>30:', len(t))
if 'GLAT' not in t.colnames:
    sc = SkyCoord(t['RAJ2000'], t['DEJ2000'], unit=u.deg)
    t['GLON'] = sc.galactic.l.deg; t['GLAT'] = sc.galactic.b.deg
t = t[np.abs(t['GLAT']) > 5]
print('AKARI off-plane |b|>5, S90>30:', len(t))

def col(c): return np.nan_to_num(np.array(t[c], float))
s65, s90, s140, s160 = col('S65'), col('S90'), col('S140'), col('S160')
Fb = fbol(s65, s90, s140, s160)
t['F_bol'] = Fb; t['L_at3kpc'] = Fb * 4 * np.pi * (3 * kpc) ** 2 / Lsun
# FIS lacks short-IR colors; use FIR slope S90>S160? keep simple: gate on F_bol only,
# then vet. (FIR-only can't reject galaxies by color; cross-id does that.)
gate = Fb >= F_THRESH
cand = t[gate]
print(f'AKARI off-plane F_bol>=thresh: {len(cand)}')
cand.sort('F_bol'); cand = cand[::-1]
cand.write('offplane_akari_cand.fits', overwrite=True)

# cross-match to IRAS off-plane candidates
iras = Table.read('offplane_iras_v2_cand.fits')
ic = SkyCoord(iras['_ra'], iras['_dec'], unit=u.deg)
ac = SkyCoord(cand['RAJ2000'], cand['DEJ2000'], unit=u.deg)
idx, d2d, _ = ac.match_to_catalog_sky(ic)
cand['iras_sep'] = d2d.arcsec
cand['iras_match'] = [str(iras['IRAS'][i]) if s < 120 else '--'
                      for i, s in zip(idx, d2d.arcsec)]
print(f'\nAKARI off-plane candidates (sorted F_bol); iras_match if <120"):')
for r in cand:
    print(f"  {r['objName']:18s} l={r['GLON']:7.2f} b={r['GLAT']:+6.2f} "
          f"L@3kpc={r['L_at3kpc']:.2e}  IRAS={r['iras_match']}")
# any AKARI candidate with NO IRAS match = potential IRAS-missed source
new = cand[[m == '--' for m in cand['iras_match']]]
print(f'\nAKARI off-plane candidates with NO IRAS off-plane match: {len(new)}')
for r in new:
    print(f"  {r['objName']:18s} l={r['GLON']:7.2f} b={r['GLAT']:+6.2f} L@3kpc={r['L_at3kpc']:.2e}")