#!/usr/bin/env python3 """Offline Gaia DR3 tutorial. Python >=3.10, standard library only. A deliberately selected 24-row catalogue subset, not a Galactic population sample. Positions are ICRS degrees at J2016.0. Parallax and tangent position errors are mas; pmra means mu_alpha*cos(delta), in mas/Julian year. No zero-point correction is made. """ import argparse import csv import hashlib import json import math from pathlib import Path SCRIPT_DIRECTORY = Path(__file__).resolve().parent ROOT = SCRIPT_DIRECTORY.parents[2] if len(SCRIPT_DIRECTORY.parents) > 2 else SCRIPT_DIRECTORY DATA = ROOT / 'docs/astronomy/ac01/data' PARAMETERS = ['ra', 'dec', 'parallax', 'pmra', 'pmdec'] def covariance(row): errors = [float(row[name + '_error']) for name in PARAMETERS] result = [[0.0] * 5 for _ in range(5)] for i in range(5): result[i][i] = errors[i] ** 2 for j in range(i + 1, 5): result[i][j] = result[j][i] = float(row[PARAMETERS[i] + '_' + PARAMETERS[j] + '_corr']) * errors[i] * errors[j] return result def cholesky(matrix): factor = [[0.0] * len(matrix) for _ in matrix] for i in range(len(matrix)): for j in range(i + 1): remainder = matrix[i][j] - sum(factor[i][k] * factor[j][k] for k in range(j)) if i == j: if remainder <= 0: raise ValueError('Catalogue covariance is not positive definite') factor[i][j] = math.sqrt(remainder) else: factor[i][j] = remainder / factor[j][j] return factor def distance_quantiles(parallax_mas, sigma_mas, length_pc, intervals=32768): # Midpoint quadrature in log distance: dr = r d(log r). log_min, log_max = math.log(0.1), math.log(100000.0) step = (log_max - log_min) / intervals grid = [math.exp(log_min + (i + 0.5) * step) for i in range(intervals)] log_mass = [3 * math.log(r) - r / length_pc - 0.5 * ((parallax_mas - 1000 / r) / sigma_mas) ** 2 for r in grid] maximum = max(log_mass) mass = [math.exp(value - maximum) for value in log_mass] total = math.fsum(mass) quantiles = [] for target in (0.16, 0.5, 0.84): accumulated = 0.0 for i, weight in enumerate(mass): previous = accumulated accumulated += weight if accumulated >= target * total: fraction = (target * total - previous) / weight quantiles.append(math.exp(log_min + (i + fraction) * step)) break return quantiles def calculate(): with (DATA / 'gaia-dr3-sample.csv').open() as source: rows = list(csv.DictReader(source)) output = [] for row in rows: matrix = covariance(row) cholesky(matrix) parallax, sigma = float(row['parallax']), float(row['parallax_error']) posteriors = {str(length): distance_quantiles(parallax, sigma, length) for length in (500, 1350, 3000)} distance = posteriors['1350'][1] ra, dec = math.radians(float(row['ra'])), math.radians(float(row['dec'])) east_variance = matrix[0][0] + 40 * matrix[0][3] + 400 * matrix[3][3] output.append({ 'source_id': row['source_id'], 'parallax_mas': parallax, 'parallax_error_mas': sigma, 'inverse_distance_pc': 1000 / parallax if parallax > 0 else None, 'distance_quantiles_pc_by_prior_length': posteriors, 'icrs_position_pc_at_radial_median': [distance * math.cos(dec) * math.cos(ra), distance * math.cos(dec) * math.sin(ra), distance * math.sin(dec)], 'astrometric_covariance': matrix, 'propagated_east_sigma_2036_mas': math.sqrt(east_variance), 'propagated_east_sigma_2036_ignoring_correlation_mas': math.sqrt(matrix[0][0] + 400 * matrix[3][3]), 'ruwe': float(row['ruwe']), }) return {'release': 'Gaia DR3', 'reference_epoch': 'J2016.0', 'input_sha256': hashlib.sha256((DATA / 'gaia-dr3-sample.csv').read_bytes()).hexdigest(), 'measurement_provenance': 'Gaia DR3 catalogue parameters; not Gaia epoch data', 'prior': 'Teaching EDSD prior r^2 exp(-r/L); not the Bailer-Jones EDR3 catalogue prior', 'zero_point_correction': 'none; illustrative single-source inference, not calibrated population science', 'row_count': len(rows), 'negative_parallax_count': sum(float(row['parallax']) < 0 for row in rows), 'rows': output} def validate(result): assert result['row_count'] == 24 and result['negative_parallax_count'] == 3 for row in result['rows']: for quantiles in row['distance_quantiles_pc_by_prior_length'].values(): assert 0 < quantiles[0] < quantiles[1] < quantiles[2] < 100000 coarse = distance_quantiles(row['parallax_mas'], row['parallax_error_mas'], 1350, 16384) fine = row['distance_quantiles_pc_by_prior_length']['1350'] assert max(abs(a / b - 1) for a, b in zip(coarse, fine)) < 0.0003 # When the likelihood is constant in distance, this prior is Gamma(3,L). # Its analytic CDF independently checks the grid quantiles. no_information = distance_quantiles(0, 1e10, 1000) for distance, probability in zip(no_information, (0.16, 0.5, 0.84)): x = distance / 1000 assert abs(1 - math.exp(-x) * (1 + x + x*x/2) - probability) < 1e-6 if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('--check', action='store_true', help='Check committed output without writing it') parser.add_argument('--data-dir', type=Path, help='Folder containing the downloaded CSV, provenance.json and saved distance-results.json') args = parser.parse_args() args_data = args.data_dir if args_data is not None: DATA = args_data.resolve() provenance = json.loads((DATA / 'provenance.json').read_text()) assert hashlib.sha256((DATA / 'gaia-dr3-sample.csv').read_bytes()).hexdigest() == provenance['sha256'], 'Input checksum differs from frozen provenance' result = calculate() validate(result) path = DATA / 'distance-results.json' serialized = json.dumps(result, indent=2, ensure_ascii=False) + '\n' if args.check: assert path.read_text() == serialized, 'Saved result differs; rerun without --check and inspect the changes' else: path.write_text(serialized) print('PASS: 24 real Gaia DR3 rows; 3 negative parallaxes retained; covariance positive definite; quadrature convergence and analytic prior-limit checks passed.')