#!/usr/bin/env python3 """Independent pinned Mandel-Agol reference through installed batman; no Lab imports.""" import argparse from importlib.metadata import version import json from pathlib import Path import warnings # batman 2.5.3 imports its tests; one upstream docstring has an invalid escape. # Limit this workaround to that import-time spelling warning, not numerical warnings. with warnings.catch_warnings(): warnings.filterwarnings("ignore", message="invalid escape sequence.*", category=DeprecationWarning) warnings.filterwarnings("ignore", message="invalid escape sequence.*", category=SyntaxWarning) import batman import numpy as np from scipy.special import ellipk, ellipe from scipy.integrate import quad def central_edge_analytic(radius, linear, quadratic): # Mandel-Agol Table1 caseV / Eq7: planet edge through stellar center, radius<1/2. modulus_squared = 4 * radius * radius lam = 1 / 3 + 2 / (9 * np.pi) * (4 * (2 * radius * radius - 1) * ellipe(modulus_squared) + (1 - 4 * radius * radius) * ellipk(modulus_squared)) eta = 1.5 * radius ** 4 return 1 - ((1 - linear - 2 * quadratic) * radius ** 2 + (linear + 2 * quadratic) * lam + quadratic * eta) / (1 - linear / 3 - quadratic / 6) def cartesian_integral(radius, distance, linear, quadratic): # Independent nested Cartesian area integral; the browser uses annular arcs. lower, upper = max(-1., distance - radius), min(1., distance + radius) if lower >= upper: return 1. split = (1 - radius * radius + distance * distance) / (2 * distance) if distance else upper points = [lower] + ([split] if lower < split < upper else []) + [upper] def strip(horizontal): limit = min(np.sqrt(max(0., 1 - horizontal * horizontal)), np.sqrt(max(0., radius * radius - (horizontal - distance) ** 2))) def vertical(angle): ordinate = limit * np.sin(angle) missing = 1 - np.sqrt(max(0., 1 - horizontal * horizontal - ordinate * ordinate)) return (1 - linear * missing - quadratic * missing ** 2) * limit * np.cos(angle) return 2 * quad(vertical, 0, np.pi / 2, epsabs=2e-12, epsrel=2e-12)[0] blocked = sum(quad(lambda angle: strip(left + (right - left) * np.sin(angle) ** 2) * 2 * (right - left) * np.sin(angle) * np.cos(angle), 0, np.pi / 2, epsabs=2e-12, epsrel=2e-12)[0] for left, right in zip(points[:-1], points[1:])) return float(1 - blocked / (np.pi * (1 - linear / 3 - quadratic / 6))) def calculate(): rows = [] for radius in [0.02, 0.07, 0.12, 0.18, 0.35]: for linear, quadratic in [(0., 0.), (.45, .2), (.8, .15), (.05, .85)]: separations = sorted(set([0., radius, max(0., radius - 1e-7), radius + 1e-7, 1 - radius - 1e-7, 1 - radius, 1 - radius + 1e-7, 1., 1 + radius - 1e-7, 1 + radius, 1 + radius + 1e-7, .1, .5, 1.25])) params = batman.TransitParams() params.t0, params.per, params.rp, params.a = 0., 3., radius, 10. params.inc, params.ecc, params.w = 90., 0., 90. params.limb_dark, params.u = 'quadratic', [linear, quadratic] times = np.arcsin(np.array(separations) / params.a) * params.per / (2 * np.pi) model = batman.TransitModel(params, times) fluxes = model.light_curve(params) for separation, flux in zip(separations, fluxes): rows.append({'radius_ratio': radius, 'separation': separation, 'limb_linear': linear, 'limb_quadratic': quadratic, 'host_flux': float(central_edge_analytic(radius, linear, quadratic) if separation == radius else flux), 'batman_raw_host_flux': float(flux), 'reference_method': 'Mandel-Agol Table1 caseV/Eq7 with SciPy elliptic integrals' if separation == radius else 'batman analytic quadratic', **({'cartesian_host_flux': cartesian_integral(radius, separation, linear, quadratic)} if separation == radius or (separation == 1 and linear == .45) else {})}) return {'schema': 'ac17-mandel-agol-reference-v1', 'kind': 'independent analytic synthetic reference, not observed data', 'distribution': 'batman-package', 'distribution_version': version('batman-package'), 'module_reported_version': batman.__version__, 'numpy_version': np.__version__, 'source': 'Mandel and Agol astro-ph/0210099v1 Sections2/4, Eq1 and quadratic Table1', 'convention': 'circular orbit, front half, per3 days, a/Rstar10, impact0; times chosen to give declared separation', 'reference_failure_scope': 'Raw batman z=p/roundoff-near-z=p values are retained; exact z=p uses separately transcribed published special-case elliptic formula, checked by independent Cartesian area integration. No coefficients or browser quadrature are tuned to the failed library values.', 'rows': rows} def exposure_references(inputs): results = [] for case in inputs['cases']: settings = case['settings'] radius, impact, axis = [settings[key] for key in ['radiusRatio','impact','scaledAxis']] period, seconds = settings['periodDays'], settings['exposureSeconds'] contact = period/(2*np.pi)*np.arcsin(np.sqrt(max(0.,(1+radius)**2-impact**2)/(axis**2-impact**2))) centers = np.array([0., contact, contact-seconds/172800, contact+seconds/172800, contact*.5]) params = batman.TransitParams() params.t0, params.per, params.rp, params.a = 0., period, radius, axis params.inc, params.ecc, params.w = float(np.degrees(np.arccos(impact/axis))), 0., 90. params.limb_dark, params.u = 'quadratic',[settings['limbLinear'],settings['limbQuadratic']] refined = [] for count in [4096,8192]: offsets = ((np.arange(count)+.5)/count-.5)*seconds/86400 times = (centers[:,None]+offsets[None,:]).ravel() model = batman.TransitModel(params,times,nthreads=1) values = model.light_curve(params).reshape(len(centers),count) refined.append(np.mean(values,axis=1)) results.append({'name':case['name'],'settings':settings, 'kind':'independent analytic kernel plus composite-midpoint exposure integration', 'coarse_time_samples':4096,'fine_time_samples':8192, 'maximum_time_refinement_difference':float(np.max(abs(refined[0]-refined[1]))), 'rows':[{'time_days':float(time),'host_averaged':float(flux), 'blended_averaged':float((flux+settings['contamination'])/(1+settings['contamination']))} for time,flux in zip(centers,refined[1])]}) return results def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--out', type=Path, default=Path('analytic-reference-rerun.json')) parser.add_argument('--input', type=Path, default=Path('synthetic-input.json')) args = parser.parse_args() result = calculate() if args.input.exists(): raw = args.input.read_bytes() import hashlib result['synthetic_input_sha256'] = hashlib.sha256(raw).hexdigest() result['exposure_references'] = exposure_references(json.loads(raw)) args.out.parent.mkdir(parents=True, exist_ok=True) args.out.write_text(json.dumps(result, indent=2, allow_nan=False) + '\n') print('Analytic reference rows:', len(result['rows']), 'batman distribution:', result['distribution_version']) if __name__ == '__main__': main()