#!/usr/bin/env python """Final P&G26 reassessment: volume-limited (<3 kpc) luminosity functions of (a) our embedded UC/compact-HII+MYSO sample (N=92), (b) WISE-catalog known HII regions (all sizes, mostly extended) with L from our own WISE22 aperture photometry, calibrated L22->Lbol on our sample, (c) ALS II OB stars (SpType -> L). Deficit test: phi_embedded/phi_OB vs L (M11/P&G lifetime ratio); does adding extended HII regions fill the high-L end?""" import csv import pickle import re import numpy as np import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt # ---------- (a) our embedded sample: replicate write_sample_lists loading ---------- def fnum(x): try: return float(x) except (TypeError, ValueError): return np.nan ours = [] for o in csv.DictReader(open('master_distances.csv')): ours.append(dict(Name=o['Name'], d=float(o['d_best']), L=float(o['Lbol_best']), typ=o['Type'])) ext = {o['Name']: o for o in csv.DictReader(open('distances_ext.csv'))} cand = {o['Name']: o for o in csv.DictReader(open('candidate_distances_ext.csv'))} for nm, e in ext.items(): if nm not in cand: continue c = cand[nm] db, dr, Lr = fnum(e['d_best']), fnum(c['Dist_RMS']), fnum(c['Lbol_RMS']) if not (np.isfinite(db) and np.isfinite(dr) and dr > 0): continue ours.append(dict(Name=nm, d=db, L=Lr * (db / dr) ** 2, typ=c['Type'])) for o in csv.DictReader(open('curated_nearby_sample.csv')): if 'REVEALED-EXCLUDE' in str(o.get('note', '')) or 'marginal' in str(o.get('note', '')).lower(): continue ours.append(dict(Name=o['CommonName'], d=float(o['Dist_kpc']), L=fnum(o['Lbol']), typ=o['Type'])) seen = set() uniq = [] for r in ours: if r['Name'] in seen: continue seen.add(r['Name']) uniq.append(r) emb = [r for r in uniq if r['d'] < 3.0 and np.isfinite(r['L']) and r['L'] >= 1e4] L_emb = np.array([r['L'] for r in emb]) L_emb_hii = np.array([r['L'] for r in emb if 'HII' in r['typ']]) print(f'embedded <3kpc N={len(L_emb)} (HII {len(L_emb_hii)})') # ---------- (b) WISE extended HII with self-calibrated L22 ---------- cal_f, cal_L, wise_rows = [], [], [] for r in csv.DictReader(open('wise22_photometry.csv')): if r['set'] == 'cal' and fnum(r['Lbol']) > 0 and fnum(r['f22rel']) > 0: cal_f.append(float(r['f22rel']) * float(r['d']) ** 2) cal_L.append(float(r['Lbol'])) elif r['set'] == 'wise' and fnum(r['f22rel']) > 0: wise_rows.append(r) cal_f, cal_L = np.array(cal_f), np.array(cal_L) # calibration: log L = log(f22*d^2) + C (slope fixed to 1); C = median offset C = np.median(np.log10(cal_L) - np.log10(cal_f)) scat = 1.4826 * np.median(np.abs((np.log10(cal_L) - np.log10(cal_f)) - C)) print(f'L22 calibration: N_cal={len(cal_f)}, offset C={C:.2f}, scatter {scat:.2f} dex') L_wise = np.array([10 ** (np.log10(float(r['f22rel']) * float(r['d']) ** 2) + C) for r in wise_rows]) R_wise = np.array([float(r['rad_as']) / 206265 * float(r['d']) * 1000 for r in wise_rows]) print(f'WISE known HII <3kpc with photometry: N={len(L_wise)}; ' f'L>1e4: {(L_wise>1e4).sum()}, L>1e5: {(L_wise>1e5).sum()}') # exclude WISE regions that ARE our embedded sources (compact, r<0.15pc) to avoid # double counting; extended = r_pc > 0.15 extended = R_wise > 0.15 print(f'extended (r>0.15pc): {extended.sum()}') # ---------- (c) ALS OB stars <3 kpc ---------- SPT_L = { # log L/Lsun, luminosity class V ZAMS-ish 'O3': 5.8, 'O4': 5.68, 'O5': 5.55, 'O5.5': 5.5, 'O6': 5.4, 'O6.5': 5.32, 'O7': 5.23, 'O7.5': 5.14, 'O8': 5.05, 'O8.5': 4.96, 'O9': 4.87, 'O9.5': 4.77, 'B0': 4.6, 'B0.5': 4.4, 'B1': 4.15, 'B1.5': 3.95, 'B2': 3.75, 'B2.5': 3.6, 'B3': 3.4} CLASS_BOOST = {'I': 0.7, 'II': 0.45, 'III': 0.25, 'IV': 0.12, 'V': 0.0} als = pickle.load(open('als_combined.pkl', 'rb')) Lob = [] for row in als: d = row['DistMode'] if not np.isfinite(d) or d <= 0 or d > 3000: continue sp = str(row['SpType']) m = re.match(r'^([OB])\s*(\d+(?:\.5)?)', sp) if not m: continue key = m.group(1) + m.group(2) if key not in SPT_L: continue lum = 0.0 mc = re.search(r'(I{1,3}|IV|V)\b', sp) if mc: lum = CLASS_BOOST.get(mc.group(1), 0.0) Lob.append(10 ** (SPT_L[key] + lum)) Lob = np.array(Lob) print(f'ALS OB <3kpc with usable SpType: N={len(Lob)}; L>1e4: {(Lob>1e4).sum()}, ' f'L>1e5: {(Lob>1e5).sum()}') # ---------- LFs + ratio ---------- bins = np.arange(4.0, 6.25, 0.25) cen = 0.5 * (bins[1:] + bins[:-1]) h_emb, _ = np.histogram(np.log10(L_emb), bins=bins) h_hii, _ = np.histogram(np.log10(L_emb_hii), bins=bins) h_wis, _ = np.histogram(np.log10(L_wise[extended & (L_wise > 1e4)]), bins=bins) h_ob, _ = np.histogram(np.log10(Lob), bins=bins) fig, axes = plt.subplots(1, 2, figsize=(12.5, 5)) ax = axes[0] ax.step(cen, h_emb, where='mid', lw=2.2, color='C0', label=f'embedded (N={h_emb.sum()})') ax.step(cen, h_wis, where='mid', lw=2.2, color='C2', label=f'WISE extended HII (N={h_wis.sum()})') ax.step(cen, h_emb + h_wis, where='mid', lw=1.5, ls='--', color='k', label='embedded+extended') ax.step(cen, h_ob, where='mid', lw=2.2, color='C3', label=f'ALS OB stars (N={h_ob.sum()})') ax.set_yscale('log') ax.set_xlabel(r'$\log L$ [L$_\odot$]') ax.set_ylabel('N per 0.25 dex (<3 kpc)') ax.legend(fontsize=8) ax.set_title('Volume-limited luminosity functions') ax = axes[1] with np.errstate(divide='ignore', invalid='ignore'): r_emb = h_emb / h_ob.astype(float) r_tot = (h_emb + h_wis) / h_ob.astype(float) e_emb = r_emb * np.sqrt(1 / np.maximum(h_emb, 1) + 1 / np.maximum(h_ob, 1)) e_tot = r_tot * np.sqrt(1 / np.maximum(h_emb + h_wis, 1) + 1 / np.maximum(h_ob, 1)) ax.errorbar(cen, r_emb, yerr=e_emb, fmt='o-', color='C0', label='embedded / OB') ax.errorbar(cen, r_tot, yerr=e_tot, fmt='s--', color='k', alpha=0.7, label='(embedded+extended HII) / OB') ax.set_yscale('log') ax.set_xlabel(r'$\log L$ [L$_\odot$]') ax.set_ylabel(r'$\phi_{\rm HII}/\phi_{\rm OB}$ ($\propto t_{\rm phase}/t_{\rm MS}$)') ax.legend(fontsize=8) ax.set_title('M11/P&G lifetime-ratio diagnostic') fig.tight_layout() fig.savefig('fig_lf_deficit_test.png', dpi=130) print('wrote fig_lf_deficit_test.png') for i, c in enumerate(cen): print(f'logL={c:.2f}: emb={h_emb[i]:3d} ext={h_wis[i]:3d} OB={h_ob[i]:4d} ' f'r_emb={r_emb[i]:.3f} r_tot={r_tot[i]:.3f}')