What happens when we stop assuming what can't be assumed?
Astrophotography by GrandpaJunkrat ----- Python code classification
Cosmology from my armchair..............................
I do it metric .............................................................. Location:
bastion@protonmail.com The Netherlands
Bortle 7 - 9
This code is released under the CC0 1.0 Universal (CC0 1.0) Public Domain Dedication.
You may copy, modify, distribute, and use it for any purpose, even commercially, without asking permission.
"""
Version 4.6 - Universal, observables-based classifier
Strict matching rules:
- 0 or >1 matches → observable = unknown
- No nearest-neighbour
- Coordinates are always required and leading
- Gaia ID is optional (if provided, it is checked against coordinates)
- C-domain (compact remnants) has priority over Y-domain
- No object names (unreliable due to nomenclature ambiguity)
"""
# Auto-install required packages when running in Google Colab
try:
import google.colab
IN_COLAB = True
except ImportError:
IN_COLAB = False
if IN_COLAB:
print("Running in Google Colab – installing required packages...")
!pip install astroquery astropy numpy --quiet
import os
from datetime import datetime
from astropy.coordinates import SkyCoord
from astroquery.simbad import Simbad
from astroquery.vizier import Vizier
from astroquery.gaia import Gaia
import astropy.units as u
import numpy as np
# ============================================================================
# CONFIGURATION
# ============================================================================
WISE_SEARCH_RADIUS_ARCSEC = 30.0
SIMBAD_SEARCH_RADIUS_ARCSEC = 1.0
GAIA_SEARCH_RADIUS_ARCSEC = 1.0
LOG_DIR = "classification_logs"
# ============================================================================
# LOGGING
# ============================================================================
def save_log_to_file(identifier, log_lines, classification, search_info):
if not os.path.exists(LOG_DIR):
os.makedirs(LOG_DIR)
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
safe_name = identifier.replace(" ", "_").replace(":", "").replace("+", "p").replace(".", "_")
filename = os.path.join(LOG_DIR, f"{safe_name}_{timestamp}.log")
with open(filename, "w", encoding="utf-8") as f:
f.write("=" * 70 + "\n")
f.write(f"Classification Report - {identifier}\n")
f.write(f"Date/Time: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
f.write(f"Search info: {search_info}\n")
f.write("=" * 70 + "\n\n")
for line in log_lines:
f.write(line + "\n")
f.write("\n" + "=" * 70 + "\n")
f.write(f"CLASSIFICATION: {classification}\n")
f.write("=" * 70 + "\n")
return filename
# ============================================================================
# PRIMARY OBSERVABLES (STRICT MATCHING)
# ============================================================================
def get_gaia_by_coordinates(ra_str, dec_str):
"""
Strict Gaia matching by coordinates:
- 0 matches → unknown
- >1 matches → unknown
- Only 1 unique source accepted
"""
coord = SkyCoord(ra_str, dec_str, unit=(u.hourangle, u.deg))
radius = GAIA_SEARCH_RADIUS_ARCSEC * u.arcsec
query = f"""
SELECT
source_id, ra, dec,
phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag,
parallax, parallax_error,
ruwe, phot_bp_rp_excess_factor
FROM gaiadr3.gaia_source
WHERE 1=CONTAINS(
POINT('ICRS', ra, dec),
CIRCLE('ICRS', {coord.ra.deg}, {coord.dec.deg}, {radius.to(u.deg).value})
)
"""
try:
job = Gaia.launch_job_async(query, dump_to_file=False)
results = job.get_results()
if len(results) != 1:
return {"success": False}
row = results[0]
target_coord = SkyCoord(row["ra"], row["dec"], unit="deg")
separation = coord.separation(target_coord).arcsec
photometry = {
"G": row["phot_g_mean_mag"],
"BP": row["phot_bp_mean_mag"],
"RP": row["phot_rp_mean_mag"],
}
color = None
if row["phot_bp_mean_mag"] and row["phot_rp_mean_mag"]:
color = row["phot_bp_mean_mag"] - row["phot_rp_mean_mag"]
return {
"success": True,
"source_id": row["source_id"],
"ra": row["ra"],
"dec": row["dec"],
"coord_sep_arcsec": separation,
"photometry": photometry,
"color": color,
"parallax_mas": row["parallax"],
"parallax_error": row["parallax_error"],
"ruwe": row["ruwe"],
"bp_rp_excess": row["phot_bp_rp_excess_factor"],
}
except Exception as e:
print(f" DEBUG: Gaia query error: {e}")
return {"success": False}
def get_gaia_by_id(source_id):
"""
Direct Gaia query by source_id (100% reliable, no position uncertainty)
"""
query = f"""
SELECT
source_id, ra, dec,
phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag,
parallax, parallax_error,
ruwe, phot_bp_rp_excess_factor
FROM gaiadr3.gaia_source
WHERE source_id = {source_id}
"""
try:
job = Gaia.launch_job_async(query, dump_to_file=False)
results = job.get_results()
if len(results) != 1:
return {"success": False}
row = results[0]
photometry = {
"G": row["phot_g_mean_mag"],
"BP": row["phot_bp_mean_mag"],
"RP": row["phot_rp_mean_mag"],
}
color = None
if row["phot_bp_mean_mag"] and row["phot_rp_mean_mag"]:
color = row["phot_bp_mean_mag"] - row["phot_rp_mean_mag"]
return {
"success": True,
"source_id": row["source_id"],
"ra": row["ra"],
"dec": row["dec"],
"coord_sep_arcsec": 0.0,
"photometry": photometry,
"color": color,
"parallax_mas": row["parallax"],
"parallax_error": row["parallax_error"],
"ruwe": row["ruwe"],
"bp_rp_excess": row["phot_bp_rp_excess_factor"],
}
except Exception as e:
print(f" DEBUG: Gaia ID query error: {e}")
return {"success": False}
def get_wise_primary(ra_str, dec_str):
"""
Strict WISE matching:
- 0 matches → unknown
- >1 matches → unknown
"""
coord = SkyCoord(ra_str, dec_str, unit=(u.hourangle, u.deg))
vizier = Vizier(row_limit=-1)
radius = WISE_SEARCH_RADIUS_ARCSEC * u.arcsec
try:
result = vizier.query_region(coord, radius=radius, catalog="II/328/allwise")
if not result or len(result[0]) != 1:
return {"success": False}
row = result[0][0]
w1 = float(row["W1mag"])
w2 = float(row["W2mag"])
w3 = float(row["W3mag"])
w4 = float(row["W4mag"])
photometry = {"W1": w1, "W2": w2, "W3": w3, "W4": w4}
colors = {"W1W2": w1 - w2, "W2W3": w2 - w3, "W3W4": w3 - w4}
return {
"success": True,
"photometry": photometry,
"colors": colors,
"qual": row["qph"],
"errors": {
"W1": float(row["e_W1mag"]),
"W2": float(row["e_W2mag"]),
"W3": float(row["e_W3mag"]),
"W4": float(row["e_W4mag"]),
},
}
except Exception:
return {"success": False}
def get_simbad_primary(ra_str, dec_str):
"""
Strict SIMBAD matching:
- 0 matches → unknown
- >1 matches → unknown
"""
coord = SkyCoord(ra_str, dec_str, unit=(u.hourangle, u.deg))
Simbad.reset_votable_fields()
Simbad.add_votable_fields("sp_type", "otype")
try:
result = Simbad.query_region(coord, radius=SIMBAD_SEARCH_RADIUS_ARCSEC * u.arcsec)
if result is None or len(result) != 1:
print(f" DEBUG SIMBAD: No unique match (found {len(result) if result is not None else 0} sources)")
return {"success": False}
row = result[0]
sp_type = row["sp_type"].strip() if row["sp_type"] else None
otype = row["otype"].strip() if row["otype"] else None
print(f" DEBUG SIMBAD: sp_type='{sp_type}', otype='{otype}'")
return {
"success": True,
"sp_type": sp_type,
"otype": otype,
}
except Exception as e:
print(f" DEBUG SIMBAD error: {e}")
return {"success": False}
# ============================================================================
# DERIVED OBSERVABLES
# ============================================================================
def derive_ir_excess(gaia_primary, wise_primary):
if wise_primary and wise_primary.get("success"):
if wise_primary["colors"]["W2W3"] > 1.0:
return True
if gaia_primary and gaia_primary.get("success"):
if gaia_primary["bp_rp_excess"] and gaia_primary["bp_rp_excess"] > 0.2:
return True
return False
def derive_sed_peak(wise_primary):
if not wise_primary or not wise_primary.get("success"):
return None
phot = wise_primary["photometry"]
mags = [("W1", phot["W1"]), ("W2", phot["W2"]), ("W3", phot["W3"]), ("W4", phot["W4"])]
mags_sorted = sorted(mags, key=lambda x: x[1])
return mags_sorted[0][0]
def derive_tbol_placeholder():
return None
def derive_compactness(gaia_primary):
if gaia_primary and gaia_primary.get("success"):
return "compact"
return None
# ============================================================================
# DIAGNOSTIC OBSERVABLES
# ============================================================================
def diagnostic_accretion(gaia_primary):
if not gaia_primary or not gaia_primary.get("success"):
return False
if gaia_primary["ruwe"] and gaia_primary["ruwe"] > 1.4:
return True
if gaia_primary["bp_rp_excess"] and gaia_primary["bp_rp_excess"] > 0.3:
return True
return False
# ============================================================================
# COMPACT REMNANT DETECTION - EXTENDED VERSION
# ============================================================================
# Keywords for compact remnants (C-domain) - expanded
COMPACT_REMNANT_KEYWORDS = [
"WD", "white dwarf", "WhiteDwarf", "WD*", "D", "DA", "DB", "DQ", "DZ", "DC",
"DO", "PG", "sd", "subdwarf",
"neutron star", "NS", "pulsar", "PSR", "Crab", "CM Tau",
"black hole", "BH", "CV", "cataclysmic", "nova", "supernova", "SNR"
]
# Expanded list for SIMBAD otypes
COMPACT_REMNANT_OTYPES = [
"WD", "WD*", "WhiteDwarf", "D", "DA", "DB", "DQ", "DZ", "DC", "DO",
"Cataclysmic*", "CV*", "Nov*", "SN*", "SNR", "Pulsar", "Psr", "NeutronStar",
"BlackHole", "X", "gamma-ray", "High-mass*"
]
def is_compact_remnant(simbad_primary, gaia_primary=None):
"""
Detect if an object is a compact remnant (white dwarf, neutron star, black hole).
First checks SIMBAD, then falls back to Gaia RUWE + color as heuristic.
Returns: (is_remnant, type_string)
"""
# Check SIMBAD first
if simbad_primary and simbad_primary.get("success"):
sp_type = simbad_primary.get("sp_type", "")
otype = simbad_primary.get("otype", "")
if sp_type:
for keyword in COMPACT_REMNANT_KEYWORDS:
if keyword.upper() in sp_type.upper():
return True, f"SIMBAD: {sp_type}"
if otype:
for keyword in COMPACT_REMNANT_OTYPES:
if keyword.upper() in otype.upper():
return True, f"SIMBAD: {otype}"
# Fallback: Heuristic based on Gaia data for white dwarfs
if gaia_primary and gaia_primary.get("success"):
# White dwarfs typically have:
# - Very blue color (BP-RP < 0.5 or negative)
# - Small RUWE (< 1.4)
# - High proper motion (not implemented here)
color = gaia_primary.get("color")
ruwe = gaia_primary.get("ruwe")
if color is not None and ruwe is not None:
# Very blue color indicates potential white dwarf
if color < 0.5 and ruwe < 1.4:
# Additional check: absolute magnitude (would require parallax)
parallax = gaia_primary.get("parallax_mas")
if parallax and parallax > 0:
g_mag = gaia_primary["photometry"]["G"]
distance_pc = 1000.0 / parallax
abs_g = g_mag - 5 * np.log10(distance_pc) + 5
# White dwarfs typically have abs_g > 8 (faint)
if abs_g > 8:
return True, f"Gaia heuristic: color={color:.2f}, M_G={abs_g:.1f}"
return False, None
# ============================================================================
# MASS ESTIMATION
# ============================================================================
MASS_TABLE = {
"B5V": 5.5, "B6V": 4.8, "B7V": 4.2, "B8V": 3.6, "B9V": 2.8,
"A0V": 2.4, "A1V": 2.3, "A2V": 2.2, "A3V": 2.1, "A4V": 2.0,
"A5V": 1.9, "A6V": 1.8, "A7V": 1.7, "A8V": 1.6, "A9V": 1.5,
"F0V": 1.45, "F2V": 1.4, "F5V": 1.3, "F8V": 1.2, "F9V": 1.1,
"G0V": 1.08, "G2V": 1.0, "G5V": 0.95, "G8V": 0.85, "G9V": 0.80,
"K0V": 0.79, "K2V": 0.75, "K4V": 0.70, "K5V": 0.65, "K7V": 0.60,
"M0V": 0.55, "M2V": 0.45, "M4V": 0.25, "M5V": 0.18, "M6V": 0.12,
}
def mass_from_simbad(sp_type):
if not sp_type:
return None
for key, m in MASS_TABLE.items():
if key in sp_type:
return m
return None
def mass_from_gaia_photometry(gaia_primary):
if not gaia_primary or not gaia_primary.get("success"):
return None
g_mag = gaia_primary["photometry"]["G"]
parallax_mas = gaia_primary["parallax_mas"]
if parallax_mas is None or parallax_mas <= 0:
return None
distance_pc = 1000.0 / parallax_mas
abs_g = g_mag - 5 * np.log10(distance_pc) + 5
if abs_g < 2.0:
return 2.5
elif abs_g < 4.0:
return 1.2
elif abs_g < 6.0:
return 0.8
else:
return 0.5
def diagnostic_mass_estimate(gaia_primary, simbad_primary):
sp_type = simbad_primary["sp_type"] if simbad_primary and simbad_primary.get("success") else None
m = mass_from_simbad(sp_type)
if m is not None:
return m
m = mass_from_gaia_photometry(gaia_primary)
if m is not None:
return m
return 1.0
# ============================================================================
# M-AXIS (mass regime)
# ============================================================================
def get_m_label(mass):
if mass > 5.0:
return "M(high)"
elif mass > 1.5:
return "M(int)"
elif mass > 0.8:
return "M(low)"
else:
return "M(very_low)"
# ============================================================================
# MAPPING TO O, E, F
# ============================================================================
def map_to_o_label(ir_excess, is_compact_remnant=False):
if is_compact_remnant:
return "O(phot.point)"
if ir_excess:
return "O(phot.IRex.point)"
return "O(phot.noIRex.point)"
def map_to_e_label(sed_peak, is_compact_remnant=False):
if is_compact_remnant:
return "E(n.a.)"
if sed_peak in ["W3", "W4"]:
return "E(II)"
return "E(III)"
def map_to_f_label(ir_excess, accretion, is_compact_remnant=False):
if is_compact_remnant:
return "F(none)"
if ir_excess and accretion:
return "F(disc.full.acc)"
if ir_excess:
return "F(disc)"
return "F(none)"
# ============================================================================
# DOMAIN LOGIC (WITH C-DOMAIN PRIORITY)
# ============================================================================
def map_to_domain(mass, ir_excess, accretion, gaia_primary, simbad_primary):
# Check 1: Compact remnant (C-domain) - HIGHEST PRIORITY
is_remnant, remnant_type = is_compact_remnant(simbad_primary, gaia_primary)
if is_remnant:
print(f" DEBUG: Compact remnant detected: {remnant_type}")
return "C"
# Check 2: Evolved star (EV-domain)
sp_type = simbad_primary["sp_type"] if simbad_primary and simbad_primary.get("success") else None
if sp_type and any(cls in sp_type for cls in ["III", "II", "I", "giant", "supergiant"]):
return "EV"
# Check 3: Pre-equilibrium object (Y-domain)
if ir_excess and accretion:
return "Y"
# Check 4: Substellar post-contraction (XP-domain)
if mass < 0.075:
return "XP"
# Check 5: Default main sequence (MS-domain)
return "MS"
# ============================================================================
# MAIN CLASSIFICATION FUNCTION
# ============================================================================
def classify_object(ra_str, dec_str, gaia_id=None):
"""
Classify an object using coordinates (required) and optionally a Gaia ID.
Parameters:
-----------
ra_str : str
Right Ascension in hours, e.g., "21:42:46.04"
dec_str : str
Declination in degrees, e.g., "+66:05:13.80"
gaia_id : str, optional
Gaia DR3 source_id, e.g., "2218019861542693760"
Returns:
--------
str : Classification string, e.g., "C:(O(phot.point), M(very_low), E(n.a.); F(none))"
"""
log = []
# Validate: coordinates are required
if not ra_str or not dec_str:
print("ERROR: Coordinates (ra_str, dec_str) are required!")
return None
print("\n" + "=" * 70)
print(f"Coordinates: {ra_str} {dec_str}")
if gaia_id:
print(f"Provided Gaia ID: {gaia_id}")
print("-" * 70)
# Step 1: Always get Gaia data from coordinates first
gaia_from_coords = get_gaia_by_coordinates(ra_str, dec_str)
# Step 2: If a Gaia ID was provided, verify it matches the coordinates
gaia_primary = None
id_matches_coords = False
if gaia_id and gaia_from_coords.get("success"):
if str(gaia_from_coords["source_id"]) == str(gaia_id):
id_matches_coords = True
gaia_primary = gaia_from_coords
log.append(f"Provided Gaia ID matches coordinates: {gaia_id}")
else:
log.append(f"WARNING: Provided Gaia ID {gaia_id} does NOT match coordinates!")
log.append(f" Gaia from coordinates: {gaia_from_coords['source_id']}")
log.append(f" Proceeding with coordinates-only mode")
gaia_primary = gaia_from_coords
elif gaia_id and not gaia_from_coords.get("success"):
log.append(f"WARNING: No Gaia source found at coordinates, cannot verify Gaia ID {gaia_id}")
log.append(f" Proceeding with ID-only mode (no coordinate verification)")
gaia_primary = get_gaia_by_id(gaia_id)
else:
gaia_primary = gaia_from_coords
# Step 3: Determine search coordinates for WISE and SIMBAD
if gaia_primary and gaia_primary.get("success"):
search_ra_deg = gaia_primary["ra"]
search_dec_deg = gaia_primary["dec"]
coord = SkyCoord(search_ra_deg, search_dec_deg, unit="deg")
ra_hms = coord.to_string("hmsdms").split()[0]
dec_dms = coord.to_string("hmsdms").split()[1]
identifier = f"Gaia {gaia_primary['source_id']}"
else:
ra_hms = ra_str
dec_dms = dec_str
identifier = f"Coord_{ra_str.replace(':','')}_{dec_str.replace(':','')}"
# Step 4: Get WISE and SIMBAD data
wise_primary = get_wise_primary(ra_hms, dec_dms)
simbad_primary = get_simbad_primary(ra_hms, dec_dms)
# Step 5: Derived and diagnostic observables
ir_excess = derive_ir_excess(gaia_primary, wise_primary)
sed_peak = derive_sed_peak(wise_primary)
accretion = diagnostic_accretion(gaia_primary)
mass = diagnostic_mass_estimate(gaia_primary, simbad_primary)
# Step 6: Mapping to O, M, E, F (C-domain check happens inside map_to_domain)
m_label = get_m_label(mass)
# First get domain (which includes compact remnant detection)
domain = map_to_domain(mass, ir_excess, accretion, gaia_primary, simbad_primary)
# Check if we are dealing with a compact remnant
is_remnant = (domain == "C")
if is_remnant:
print(f" DEBUG: Setting is_remnant=True based on domain=C")
o_label = map_to_o_label(ir_excess, is_remnant)
e_label = map_to_e_label(sed_peak, is_remnant)
f_label = map_to_f_label(ir_excess, accretion, is_remnant)
classification = f"{domain}:({o_label}, {m_label}, {e_label}; {f_label})"
# Step 7: Logging
log.append(f"Gaia success: {gaia_primary.get('success') if gaia_primary else False}")
log.append(f"WISE success: {wise_primary.get('success')}")
log.append(f"SIMBAD success: {simbad_primary.get('success')}")
log.append(f"domain: {domain}")
log.append(f"ir_excess: {ir_excess}")
log.append(f"sed_peak: {sed_peak}")
log.append(f"accretion: {accretion}")
log.append(f"mass: {mass:.2f} M☉")
log.append(f"M-label: {m_label}")
log.append(f"O-label: {o_label}")
log.append(f"E-label: {e_label}")
log.append(f"F-label: {f_label}")
log.append(f"Final classification: {classification}")
for line in log:
print(" " + line)
search_info = f"Coordinates: {ra_str} {dec_str}" + (f", Gaia ID: {gaia_id}" if gaia_id else "")
logfile = save_log_to_file(identifier, log, classification, search_info)
print(f">>> CLASSIFICATION: {classification}")
print(f">>> Log saved: {logfile}")
return classification
"""
Version 4.6 - Universal, observables-based classifier
Strict matching rules:
- 0 or >1 matches → observable = unknown
- No nearest-neighbour
- Coordinates are always required and leading
- Gaia ID is optional (if provided, it is checked against coordinates)
- C-domain (compact remnants) has priority over Y-domain
- No object names (unreliable due to nomenclature ambiguity)
"""
import os
from datetime import datetime
from astropy.coordinates import SkyCoord
from astroquery.simbad import Simbad
from astroquery.vizier import Vizier
from astroquery.gaia import Gaia
import astropy.units as u
import numpy as np
# ============================================================================
# CONFIGURATION
# ============================================================================
WISE_SEARCH_RADIUS_ARCSEC = 30.0
SIMBAD_SEARCH_RADIUS_ARCSEC = 1.0
GAIA_SEARCH_RADIUS_ARCSEC = 1.0
LOG_DIR = "classification_logs"
# ============================================================================
# LOGGING
# ============================================================================
def save_log_to_file(identifier, log_lines, classification, search_info):
if not os.path.exists(LOG_DIR):
os.makedirs(LOG_DIR)
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
safe_name = identifier.replace(" ", "_").replace(":", "").replace("+", "p").replace(".", "_")
filename = os.path.join(LOG_DIR, f"{safe_name}_{timestamp}.log")
with open(filename, "w", encoding="utf-8") as f:
f.write("=" * 70 + "\n")
f.write(f"Classification Report - {identifier}\n")
f.write(f"Date/Time: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
f.write(f"Search info: {search_info}\n")
f.write("=" * 70 + "\n\n")
for line in log_lines:
f.write(line + "\n")
f.write("\n" + "=" * 70 + "\n")
f.write(f"CLASSIFICATION: {classification}\n")
f.write("=" * 70 + "\n")
return filename
# ============================================================================
# PRIMARY OBSERVABLES (STRICT MATCHING)
# ============================================================================
def get_gaia_by_coordinates(ra_str, dec_str):
"""
Strict Gaia matching by coordinates:
- 0 matches → unknown
- >1 matches → unknown
- Only 1 unique source accepted
"""
coord = SkyCoord(ra_str, dec_str, unit=(u.hourangle, u.deg))
radius = GAIA_SEARCH_RADIUS_ARCSEC * u.arcsec
query = f"""
SELECT
source_id, ra, dec,
phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag,
parallax, parallax_error,
ruwe, phot_bp_rp_excess_factor
FROM gaiadr3.gaia_source
WHERE 1=CONTAINS(
POINT('ICRS', ra, dec),
CIRCLE('ICRS', {coord.ra.deg}, {coord.dec.deg}, {radius.to(u.deg).value})
)
"""
try:
job = Gaia.launch_job_async(query, dump_to_file=False)
results = job.get_results()
if len(results) != 1:
return {"success": False}
row = results[0]
target_coord = SkyCoord(row["ra"], row["dec"], unit="deg")
separation = coord.separation(target_coord).arcsec
photometry = {
"G": row["phot_g_mean_mag"],
"BP": row["phot_bp_mean_mag"],
"RP": row["phot_rp_mean_mag"],
}
color = None
if row["phot_bp_mean_mag"] and row["phot_rp_mean_mag"]:
color = row["phot_bp_mean_mag"] - row["phot_rp_mean_mag"]
return {
"success": True,
"source_id": row["source_id"],
"ra": row["ra"],
"dec": row["dec"],
"coord_sep_arcsec": separation,
"photometry": photometry,
"color": color,
"parallax_mas": row["parallax"],
"parallax_error": row["parallax_error"],
"ruwe": row["ruwe"],
"bp_rp_excess": row["phot_bp_rp_excess_factor"],
}
except Exception as e:
print(f" DEBUG: Gaia query error: {e}")
return {"success": False}
def get_gaia_by_id(source_id):
"""
Direct Gaia query by source_id (100% reliable, no position uncertainty)
"""
query = f"""
SELECT
source_id, ra, dec,
phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag,
parallax, parallax_error,
ruwe, phot_bp_rp_excess_factor
FROM gaiadr3.gaia_source
WHERE source_id = {source_id}
"""
try:
job = Gaia.launch_job_async(query, dump_to_file=False)
results = job.get_results()
if len(results) != 1:
return {"success": False}
row = results[0]
photometry = {
"G": row["phot_g_mean_mag"],
"BP": row["phot_bp_mean_mag"],
"RP": row["phot_rp_mean_mag"],
}
color = None
if row["phot_bp_mean_mag"] and row["phot_rp_mean_mag"]:
color = row["phot_bp_mean_mag"] - row["phot_rp_mean_mag"]
return {
"success": True,
"source_id": row["source_id"],
"ra": row["ra"],
"dec": row["dec"],
"coord_sep_arcsec": 0.0,
"photometry": photometry,
"color": color,
"parallax_mas": row["parallax"],
"parallax_error": row["parallax_error"],
"ruwe": row["ruwe"],
"bp_rp_excess": row["phot_bp_rp_excess_factor"],
}
except Exception as e:
print(f" DEBUG: Gaia ID query error: {e}")
return {"success": False}
def get_wise_primary(ra_str, dec_str):
"""
Strict WISE matching:
- 0 matches → unknown
- >1 matches → unknown
"""
coord = SkyCoord(ra_str, dec_str, unit=(u.hourangle, u.deg))
vizier = Vizier(row_limit=-1)
radius = WISE_SEARCH_RADIUS_ARCSEC * u.arcsec
try:
result = vizier.query_region(coord, radius=radius, catalog="II/328/allwise")
if not result or len(result[0]) != 1:
return {"success": False}
row = result[0][0]
w1 = float(row["W1mag"])
w2 = float(row["W2mag"])
w3 = float(row["W3mag"])
w4 = float(row["W4mag"])
photometry = {"W1": w1, "W2": w2, "W3": w3, "W4": w4}
colors = {"W1W2": w1 - w2, "W2W3": w2 - w3, "W3W4": w3 - w4}
return {
"success": True,
"photometry": photometry,
"colors": colors,
"qual": row["qph"],
"errors": {
"W1": float(row["e_W1mag"]),
"W2": float(row["e_W2mag"]),
"W3": float(row["e_W3mag"]),
"W4": float(row["e_W4mag"]),
},
}
except Exception:
return {"success": False}
def get_simbad_primary(ra_str, dec_str):
"""
Strict SIMBAD matching:
- 0 matches → unknown
- >1 matches → unknown
"""
coord = SkyCoord(ra_str, dec_str, unit=(u.hourangle, u.deg))
Simbad.reset_votable_fields()
Simbad.add_votable_fields("sp_type", "otype")
try:
result = Simbad.query_region(coord, radius=SIMBAD_SEARCH_RADIUS_ARCSEC * u.arcsec)
if result is None or len(result) != 1:
print(f" DEBUG SIMBAD: No unique match (found {len(result) if result is not None else 0} sources)")
return {"success": False}
row = result[0]
sp_type = row["sp_type"].strip() if row["sp_type"] else None
otype = row["otype"].strip() if row["otype"] else None
print(f" DEBUG SIMBAD: sp_type='{sp_type}', otype='{otype}'")
return {
"success": True,
"sp_type": sp_type,
"otype": otype,
}
except Exception as e:
print(f" DEBUG SIMBAD error: {e}")
return {"success": False}
# ============================================================================
# DERIVED OBSERVABLES
# ============================================================================
def derive_ir_excess(gaia_primary, wise_primary):
if wise_primary and wise_primary.get("success"):
if wise_primary["colors"]["W2W3"] > 1.0:
return True
if gaia_primary and gaia_primary.get("success"):
if gaia_primary["bp_rp_excess"] and gaia_primary["bp_rp_excess"] > 0.2:
return True
return False
def derive_sed_peak(wise_primary):
if not wise_primary or not wise_primary.get("success"):
return None
phot = wise_primary["photometry"]
mags = [("W1", phot["W1"]), ("W2", phot["W2"]), ("W3", phot["W3"]), ("W4", phot["W4"])]
mags_sorted = sorted(mags, key=lambda x: x[1])
return mags_sorted[0][0]
def derive_tbol_placeholder():
return None
def derive_compactness(gaia_primary):
if gaia_primary and gaia_primary.get("success"):
return "compact"
return None
# ============================================================================
# DIAGNOSTIC OBSERVABLES
# ============================================================================
def diagnostic_accretion(gaia_primary):
if not gaia_primary or not gaia_primary.get("success"):
return False
if gaia_primary["ruwe"] and gaia_primary["ruwe"] > 1.4:
return True
if gaia_primary["bp_rp_excess"] and gaia_primary["bp_rp_excess"] > 0.3:
return True
return False
# ============================================================================
# COMPACT REMNANT DETECTION - EXTENDED VERSION
# ============================================================================
# Keywords for compact remnants (C-domain) - expanded
COMPACT_REMNANT_KEYWORDS = [
"WD", "white dwarf", "WhiteDwarf", "WD*", "D", "DA", "DB", "DQ", "DZ", "DC",
"DO", "PG", "sd", "subdwarf",
"neutron star", "NS", "pulsar", "PSR", "Crab", "CM Tau",
"black hole", "BH", "CV", "cataclysmic", "nova", "supernova", "SNR"
]
# Expanded list for SIMBAD otypes
COMPACT_REMNANT_OTYPES = [
"WD", "WD*", "WhiteDwarf", "D", "DA", "DB", "DQ", "DZ", "DC", "DO",
"Cataclysmic*", "CV*", "Nov*", "SN*", "SNR", "Pulsar", "Psr", "NeutronStar",
"BlackHole", "X", "gamma-ray", "High-mass*"
]
def is_compact_remnant(simbad_primary, gaia_primary=None):
"""
Detect if an object is a compact remnant (white dwarf, neutron star, black hole).
First checks SIMBAD, then falls back to Gaia RUWE + color as heuristic.
Returns: (is_remnant, type_string)
"""
# Check SIMBAD first
if simbad_primary and simbad_primary.get("success"):
sp_type = simbad_primary.get("sp_type", "")
otype = simbad_primary.get("otype", "")
if sp_type:
for keyword in COMPACT_REMNANT_KEYWORDS:
if keyword.upper() in sp_type.upper():
return True, f"SIMBAD: {sp_type}"
if otype:
for keyword in COMPACT_REMNANT_OTYPES:
if keyword.upper() in otype.upper():
return True, f"SIMBAD: {otype}"
# Fallback: Heuristic based on Gaia data for white dwarfs
if gaia_primary and gaia_primary.get("success"):
# White dwarfs typically have:
# - Very blue color (BP-RP < 0.5 or negative)
# - Small RUWE (< 1.4)
# - High proper motion (not implemented here)
color = gaia_primary.get("color")
ruwe = gaia_primary.get("ruwe")
if color is not None and ruwe is not None:
# Very blue color indicates potential white dwarf
if color < 0.5 and ruwe < 1.4:
# Additional check: absolute magnitude (would require parallax)
parallax = gaia_primary.get("parallax_mas")
if parallax and parallax > 0:
g_mag = gaia_primary["photometry"]["G"]
distance_pc = 1000.0 / parallax
abs_g = g_mag - 5 * np.log10(distance_pc) + 5
# White dwarfs typically have abs_g > 8 (faint)
if abs_g > 8:
return True, f"Gaia heuristic: color={color:.2f}, M_G={abs_g:.1f}"
return False, None
# ============================================================================
# MASS ESTIMATION
# ============================================================================
MASS_TABLE = {
"B5V": 5.5, "B6V": 4.8, "B7V": 4.2, "B8V": 3.6, "B9V": 2.8,
"A0V": 2.4, "A1V": 2.3, "A2V": 2.2, "A3V": 2.1, "A4V": 2.0,
"A5V": 1.9, "A6V": 1.8, "A7V": 1.7, "A8V": 1.6, "A9V": 1.5,
"F0V": 1.45, "F2V": 1.4, "F5V": 1.3, "F8V": 1.2, "F9V": 1.1,
"G0V": 1.08, "G2V": 1.0, "G5V": 0.95, "G8V": 0.85, "G9V": 0.80,
"K0V": 0.79, "K2V": 0.75, "K4V": 0.70, "K5V": 0.65, "K7V": 0.60,
"M0V": 0.55, "M2V": 0.45, "M4V": 0.25, "M5V": 0.18, "M6V": 0.12,
}
def mass_from_simbad(sp_type):
if not sp_type:
return None
for key, m in MASS_TABLE.items():
if key in sp_type:
return m
return None
def mass_from_gaia_photometry(gaia_primary):
if not gaia_primary or not gaia_primary.get("success"):
return None
g_mag = gaia_primary["photometry"]["G"]
parallax_mas = gaia_primary["parallax_mas"]
if parallax_mas is None or parallax_mas <= 0:
return None
distance_pc = 1000.0 / parallax_mas
abs_g = g_mag - 5 * np.log10(distance_pc) + 5
if abs_g < 2.0:
return 2.5
elif abs_g < 4.0:
return 1.2
elif abs_g < 6.0:
return 0.8
else:
return 0.5
def diagnostic_mass_estimate(gaia_primary, simbad_primary):
sp_type = simbad_primary["sp_type"] if simbad_primary and simbad_primary.get("success") else None
m = mass_from_simbad(sp_type)
if m is not None:
return m
m = mass_from_gaia_photometry(gaia_primary)
if m is not None:
return m
return 1.0
# ============================================================================
# M-AXIS (mass regime)
# ============================================================================
def get_m_label(mass):
if mass > 5.0:
return "M(high)"
elif mass > 1.5:
return "M(int)"
elif mass > 0.8:
return "M(low)"
else:
return "M(very_low)"
# ============================================================================
# MAPPING TO O, E, F
# ============================================================================
def map_to_o_label(ir_excess, is_compact_remnant=False):
if is_compact_remnant:
return "O(phot.point)"
if ir_excess:
return "O(phot.IRex.point)"
return "O(phot.noIRex.point)"
def map_to_e_label(sed_peak, is_compact_remnant=False):
if is_compact_remnant:
return "E(n.a.)"
if sed_peak in ["W3", "W4"]:
return "E(II)"
return "E(III)"
def map_to_f_label(ir_excess, accretion, is_compact_remnant=False):
if is_compact_remnant:
return "F(none)"
if ir_excess and accretion:
return "F(disc.full.acc)"
if ir_excess:
return "F(disc)"
return "F(none)"
# ============================================================================
# DOMAIN LOGIC (WITH C-DOMAIN PRIORITY)
# ============================================================================
def map_to_domain(mass, ir_excess, accretion, gaia_primary, simbad_primary):
# Check 1: Compact remnant (C-domain) - HIGHEST PRIORITY
is_remnant, remnant_type = is_compact_remnant(simbad_primary, gaia_primary)
if is_remnant:
print(f" DEBUG: Compact remnant detected: {remnant_type}")
return "C"
# Check 2: Evolved star (EV-domain)
sp_type = simbad_primary["sp_type"] if simbad_primary and simbad_primary.get("success") else None
if sp_type and any(cls in sp_type for cls in ["III", "II", "I", "giant", "supergiant"]):
return "EV"
# Check 3: Pre-equilibrium object (Y-domain)
if ir_excess and accretion:
return "Y"
# Check 4: Substellar post-contraction (XP-domain)
if mass < 0.075:
return "XP"
# Check 5: Default main sequence (MS-domain)
return "MS"
# ============================================================================
# MAIN CLASSIFICATION FUNCTION
# ============================================================================
def classify_object(ra_str, dec_str, gaia_id=None):
"""
Classify an object using coordinates (required) and optionally a Gaia ID.
Parameters:
-----------
ra_str : str
Right Ascension in hours, e.g., "21:42:46.04"
dec_str : str
Declination in degrees, e.g., "+66:05:13.80"
gaia_id : str, optional
Gaia DR3 source_id, e.g., "2218019861542693760"
Returns:
--------
str : Classification string, e.g., "C:(O(phot.point), M(very_low), E(n.a.); F(none))"
"""
log = []
# Validate: coordinates are required
if not ra_str or not dec_str:
print("ERROR: Coordinates (ra_str, dec_str) are required!")
return None
print("\n" + "=" * 70)
print(f"Coordinates: {ra_str} {dec_str}")
if gaia_id:
print(f"Provided Gaia ID: {gaia_id}")
print("-" * 70)
# Step 1: Always get Gaia data from coordinates first
gaia_from_coords = get_gaia_by_coordinates(ra_str, dec_str)
# Step 2: If a Gaia ID was provided, verify it matches the coordinates
gaia_primary = None
id_matches_coords = False
if gaia_id and gaia_from_coords.get("success"):
if str(gaia_from_coords["source_id"]) == str(gaia_id):
id_matches_coords = True
gaia_primary = gaia_from_coords
log.append(f"Provided Gaia ID matches coordinates: {gaia_id}")
else:
log.append(f"WARNING: Provided Gaia ID {gaia_id} does NOT match coordinates!")
log.append(f" Gaia from coordinates: {gaia_from_coords['source_id']}")
log.append(f" Proceeding with coordinates-only mode")
gaia_primary = gaia_from_coords
elif gaia_id and not gaia_from_coords.get("success"):
log.append(f"WARNING: No Gaia source found at coordinates, cannot verify Gaia ID {gaia_id}")
log.append(f" Proceeding with ID-only mode (no coordinate verification)")
gaia_primary = get_gaia_by_id(gaia_id)
else:
gaia_primary = gaia_from_coords
# Step 3: Determine search coordinates for WISE and SIMBAD
if gaia_primary and gaia_primary.get("success"):
search_ra_deg = gaia_primary["ra"]
search_dec_deg = gaia_primary["dec"]
coord = SkyCoord(search_ra_deg, search_dec_deg, unit="deg")
ra_hms = coord.to_string("hmsdms").split()[0]
dec_dms = coord.to_string("hmsdms").split()[1]
identifier = f"Gaia {gaia_primary['source_id']}"
else:
ra_hms = ra_str
dec_dms = dec_str
identifier = f"Coord_{ra_str.replace(':','')}_{dec_str.replace(':','')}"
# Step 4: Get WISE and SIMBAD data
wise_primary = get_wise_primary(ra_hms, dec_dms)
simbad_primary = get_simbad_primary(ra_hms, dec_dms)
# Step 5: Derived and diagnostic observables
ir_excess = derive_ir_excess(gaia_primary, wise_primary)
sed_peak = derive_sed_peak(wise_primary)
accretion = diagnostic_accretion(gaia_primary)
mass = diagnostic_mass_estimate(gaia_primary, simbad_primary)
# Step 6: Mapping to O, M, E, F (C-domain check happens inside map_to_domain)
m_label = get_m_label(mass)
# First get domain (which includes compact remnant detection)
domain = map_to_domain(mass, ir_excess, accretion, gaia_primary, simbad_primary)
# Check if we are dealing with a compact remnant
is_remnant = (domain == "C")
if is_remnant:
print(f" DEBUG: Setting is_remnant=True based on domain=C")
o_label = map_to_o_label(ir_excess, is_remnant)
e_label = map_to_e_label(sed_peak, is_remnant)
f_label = map_to_f_label(ir_excess, accretion, is_remnant)
classification = f"{domain}:({o_label}, {m_label}, {e_label}; {f_label})"
# Step 7: Logging
log.append(f"Gaia success: {gaia_primary.get('success') if gaia_primary else False}")
log.append(f"WISE success: {wise_primary.get('success')}")
log.append(f"SIMBAD success: {simbad_primary.get('success')}")
log.append(f"domain: {domain}")
log.append(f"ir_excess: {ir_excess}")
log.append(f"sed_peak: {sed_peak}")
log.append(f"accretion: {accretion}")
log.append(f"mass: {mass:.2f} M☉")
log.append(f"M-label: {m_label}")
log.append(f"O-label: {o_label}")
log.append(f"E-label: {e_label}")
log.append(f"F-label: {f_label}")
log.append(f"Final classification: {classification}")
for line in log:
print(" " + line)
search_info = f"Coordinates: {ra_str} {dec_str}" + (f", Gaia ID: {gaia_id}" if gaia_id else "")
logfile = save_log_to_file(identifier, log, classification, search_info)
print(f">>> CLASSIFICATION: {classification}")
print(f">>> Log saved: {logfile}")
return classification
# ============================================================================
# ENTRY POINT
# ============================================================================
if __name__ == "__main__":
print("=" * 70)
print("Version 4.5 – Universal, observables‑based classifier")
print("Coordinates are required. Gaia ID is optional.")
print("C‑domain (compact remnants) has priority over all other domains.")
print("=" * 70)
# ----------------------------------------------------------------------
# EXAMPLES
# ----------------------------------------------------------------------
# classify_object("21:42:46.04", "+66:05:13.80")
# classify_object("21:42:46.04", "+66:05:13.80", gaia_id="2218019861542693760")
# classify_object("05:34:31.947", "+22:00:52.15", gaia_id="3403818172572314624")
# ----------------------------------------------------------------------
# TEST OBJECTS
# ----------------------------------------------------------------------
print("\n" + "-" * 70)
print("TEST 1: White dwarf (expected: C‑domain)")
print("-" * 70)
classify_object("21:43:21.6078", "+66:16:03.0438")
print("\n" + "-" * 70)
print("TEST 2: Crab Pulsar (expected: C‑domain)")
print("-" * 70)
classify_object("05:34:31.947", "+22:00:52.15", gaia_id="3403818172572314624")
# ----------------------------------------------------------------------
# ADD YOUR OWN OBJECTS BELOW
# ----------------------------------------------------------------------
# classify_object("YOUR_RA", "YOUR_DEC")
# classify_object("YOUR_RA", "YOUR_DEC", gaia_id="YOUR_GAIA_ID")
O–M–E–F Classifier – User Guide (Version 4.5)
Author: H.N. van Roon
License: CC0 1.0 Universal (Public Domain)
What is this?
This Python script implements the O–M–E–F universal classifier, a fully observables‑based system that assigns any point‑like astronomical source to a physical domain (C, EV, Y, XP, MS) and a set of four orthogonal axes:
O — observational appearance
M — mass regime
E — evolutionary stage
F — flux origin / circumstellar environment
The classifier is coordinate‑based, strictly reproducible, and avoids ambiguous legacy labels such as “YSO”, “T Tauri”, “Herbig Ae/Be”, or “Class I/II/III”.
It uses only direct observables from Gaia, WISE, and SIMBAD — never object names.
Key features
Gaia DR3 as primary astrometric and photometric reference
WISE mid‑IR photometry for SED and IR‑excess diagnostics
SIMBAD spectral types and otypes for physical domain identification
Strict matching rules: 0 or >1 matches → observable = unknown
C‑domain priority: compact remnants override all other domains
Automatic logfile export to classification_logs/
No object names — only coordinates (prevents cross‑match ambiguity)
How to use it in Google Colab
Open Google Colab and create a new notebook.
Copy the entire script into the first code cell.
Run the cell (Shift + Enter).
Required packages (astroquery, astropy, numpy) will be installed automatically.
The classifier will run a few built‑in test objects (white dwarf, pulsar).
To classify your own objects, call:
python
classify_object("RA", "Dec")
classify_object("RA", "Dec", gaia_id="YOUR_GAIA_ID")
Example:
python
classify_object("21:42:46.04", "+66:05:13.80")
classify_object("05:34:31.947", "+22:00:52.15", gaia_id="3403818172572314624")
Important notes
Coordinates only
The classifier does not accept object names (e.g., “BD+65 1636”).
Names are ambiguous and cause cross‑match errors.
Use coordinates in:
HH:MM:SS.ss
±DD:MM:SS.ss
Search radii (strict matching)
Catalog Radius Purpose
Gaia 1 arcsec High‑precision astrometry; avoids false matches
WISE 30 arcsec Accounts for epoch differences and proper motion
SIMBAD 1 arcsec Reliable spectral type / otype lookup
Customising search radii
All search radii used by the classifier can be modified directly in the script.
They are defined at the top of the code in the configuration block:
WISE_SEARCH_RADIUS_ARCSEC = 30.0
SIMBAD_SEARCH_RADIUS_ARCSEC = 1.0
GAIA_SEARCH_RADIUS_ARCSEC = 1.0
Gaia ID is optional
If provided, it must match the coordinate‑based Gaia source.
If omitted, the classifier retrieves the Gaia source automatically.
Internet connection required
The script queries Gaia DR3, WISE, and SIMBAD online.
Logfile
A detailed classification log is saved to:
classification_logs/
