#!/usr/bin/env python
"""Publication figure: spectral index at the ionizing source for ALMA-IMF HII
regions (peak vs common-beam), with free-free / dust regimes."""
import numpy as np
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
from astropy.table import Table
LOG=open('fig_alpha.log','w')
def p(*a):
    s=' '.join(str(x) for x in a);print(s);print(s,file=LOG);LOG.flush()

cb=Table.read('almaimf_alpha_commonbeam.fits')
pk=Table.read('almaimf_spectral_index.fits')
pk_by={str(r['Name'])+r['field']:float(r['alpha']) for r in pk}

names=[str(x) for x in cb['Name']]
acb=np.array(cb['alpha_cb'],float)
apk=np.array([pk_by.get(str(cb['Name'][i])+cb['field'][i],np.nan) for i in range(len(cb))])
# nice short labels
def short(n):
    return n.replace('-00','−').replace('+00','+').replace('G0','G').replace('G','G')
labels=[str(n) for n in names]

order=np.argsort(acb)
acb=acb[order]; apk=apk[order]; labels=[labels[i] for i in order]
y=np.arange(len(acb))

fig,ax=plt.subplots(figsize=(9,5.5))
ax.axvspan(-1.5,1.5,color='#3498db',alpha=0.12)
ax.axvspan(2.5,4.2,color='#c0392b',alpha=0.12)
ax.axvspan(1.5,2.5,color='#f39c12',alpha=0.10)
ax.axvline(-0.1,ls=':',color='k',lw=0.8); ax.text(-0.1,len(y)-0.3,'optically-thin\nfree-free',fontsize=7,ha='center',va='top')
ax.axvline(3.7,ls=':',color='k',lw=0.8); ax.text(3.7,len(y)-0.3,'pure\ndust',fontsize=7,ha='center',va='top')
ax.scatter(apk,y,marker='o',facecolor='none',edgecolor='gray',s=70,label='peak α (beam-biased)')
ax.scatter(acb,y,marker='o',color='k',s=70,label='common-beam α (trusted)')
for i in range(len(y)):
    if np.isfinite(apk[i]): ax.plot([apk[i],acb[i]],[y[i],y[i]],color='gray',lw=0.6,zorder=0)
# highlight W51-IRS2 + dust
for i,l in enumerate(labels):
    c='#c0392b' if acb[i]>2.5 else ('#16a085' if 'G049.490' in l or 'G049.488' in l else 'k')
    ax.scatter(acb[i],y[i],color=c,s=80,zorder=5)
ax.set_yticks(y); ax.set_yticklabels(labels,fontsize=8)
ax.set_xlabel('spectral index α  (S ∝ ν^α,  B3 100 GHz → B6 230 GHz)')
ax.set_xlim(-1.5,4.2)
ax.text(0,-1.3,'FREE-FREE\n(no dust at star)',color='#2471a3',ha='center',fontsize=8,fontweight='bold')
ax.text(3.35,-1.3,'DUST\n(accretion reservoir)',color='#922b21',ha='center',fontsize=8,fontweight='bold')
ax.set_title('Dust vs free-free at the ionizing source — ALMA-IMF (U)C HII regions\n'
             '7/8 pure free-free; W51-IRS2 (teal) unambiguous; 1 dust core (red)')
ax.legend(loc='lower right',fontsize=8)
plt.tight_layout(); plt.savefig('fig_alpha.png',dpi=135)
p('wrote fig_alpha.png; common-beam: free-free',int((acb<1.5).sum()),'dust',int((acb>2.5).sum()))
LOG.close()