#!/usr/bin/env python3 """Independent public AC15 analysis. Python standard library only. The real timing diagnostic is distinct from synthetic migration records. This checks a short replay with independently authored RK4, not the entire IAS15 run. """ from pathlib import Path import argparse, hashlib, importlib.util, json, math from disk_torques import ( DEFAULT_DISK, EARTH_SOLAR_MASS_RATIO, torque_components, YEARS_PER_CODE_TIME, ) from migration_orbits import external_disk_force, barycentric_invariants def derivative(state, masses, disk): result = [0.0] * 14 for body in range(3): offset = body * 4 result[offset : offset + 2] = state[offset + 2 : offset + 4] for other in range(body + 1, 3): x = state[other * 4] - state[offset] y = state[other * 4 + 1] - state[offset + 1] inverse = math.hypot(x, y) ** -3 for component, distance in enumerate([x, y]): result[offset + 2 + component] += masses[other] * distance * inverse result[other * 4 + 2 + component] -= masses[body] * distance * inverse for body in [1, 2]: offset = body * 4 acceleration, _ = external_disk_force( masses[body], [state[offset] - state[0], state[offset + 1] - state[1]], [state[offset + 2] - state[2], state[offset + 3] - state[3]], disk, ) result[offset + 2] += acceleration[0] result[offset + 3] += acceleration[1] result[12] += masses[body] * ( state[offset + 2] * acceleration[0] + state[offset + 3] * acceleration[1] ) result[13] += masses[body] * ( state[offset] * acceleration[1] - state[offset + 1] * acceleration[0] ) return result def rk4(state, duration, maximum_step, masses, disk): state = state[:] steps = math.ceil(duration / maximum_step) step = duration / steps for _ in range(steps): first = derivative(state, masses, disk) second = derivative( [a + step * b / 2 for a, b in zip(state, first)], masses, disk ) third = derivative( [a + step * b / 2 for a, b in zip(state, second)], masses, disk ) fourth = derivative([a + step * b for a, b in zip(state, third)], masses, disk) state = [ a + step * (b + 2 * c + 2 * d + e) / 6 for a, b, c, d, e in zip(state, first, second, third, fourth) ] return state def invariants(state, masses): return barycentric_invariants( [ { "mass": mass, "position": state[index * 4 : index * 4 + 2], "velocity": state[index * 4 + 2 : index * 4 + 4], } for index, mass in enumerate(masses) ] ) def reproduce(directory): manifest_path = directory / "public-reproduction.json" if manifest_path.exists(): for entry in json.loads(manifest_path.read_text())["files"]: path = directory / entry["name"] if ( path.exists() and entry.get("sha256") and entry.get("role") != "expected-output" ): assert ( hashlib.sha256(path.read_bytes()).hexdigest() == entry["sha256"] ), ("Changed input: " + path.name) spec = importlib.util.spec_from_file_location( "timing", directory / "analyze-kepler223.py" ) timing = importlib.util.module_from_spec(spec) spec.loader.exec_module(timing) observed = timing.reconstruct(directory) torque_grid = [] for slope in [0.5, 1.2]: for viscosity in [1e-9, 1e-6, 1e-4]: disk = { **DEFAULT_DISK, "surfaceDensityAtOneAu": 1e-4, "temperatureSlope": slope, "viscosity": viscosity, } torque = torque_components(4 * EARTH_SOLAR_MASS_RATIO, 1, 0, disk) torque_grid.append( { "temperatureSlope": slope, "viscosity": viscosity, "waveCoefficient": torque["wave"] / torque["normalization"], "corotationCoefficient": torque["corotation"] / torque["normalization"], "totalCoefficient": torque["total"] / torque["normalization"], } ) captured = json.loads((directory / "capture-trajectory.json").read_text()) later = json.loads((directory / "fast-trajectory.json").read_text()) assert captured["settings"] == later["settings"] summaries = [] for record in [captured, later]: window = record["history"][len(record["history"]) // 2 :] spans = [ max(row[column] for row in window) - min(row[column] for row in window) for column in [6, 7] ] turns = [] for column in [6, 7]: differences = [ window[index + 1][column] - window[index][column] for index in range(len(window) - 1) ] turns.append( sum( differences[index] * differences[index + 1] < 0 for index in range(len(differences) - 1) ) ) summaries.append( { "id": record["id"], "endInnerOrbits": record["history"][-1][0], "finalSemimajorAxesAu": record["history"][-1][1:3], "lateHalfAngleSpansDegrees": spans, "lateHalfDirectionChanges": turns, } ) assert ( max(summaries[0]["lateHalfAngleSpansDegrees"]) < 25 and min(summaries[0]["lateHalfDirectionChanges"]) > 5 ) assert min(summaries[1]["lateHalfAngleSpansDegrees"]) > 9000 settings = captured["settings"] disk = { **DEFAULT_DISK, "surfaceDensityAtOneAu": settings["surfaceDensityAtOneAu"], "aspectRatioAtOneAu": settings["aspectRatioAtOneAu"], "temperatureSlope": settings["temperatureSlope"], "viscosity": 10 ** settings["logViscosity"], "temperatureTransition": settings["temperatureTransition"], } start = captured["clip"][0][1:] + [0.0, 0.0] duration = (captured["clip"][-1][0] - captured["clip"][0][0]) * captured[ "initialInnerPeriod" ] reference = captured["clip"][-1][1:] reference_checks = [] initial_invariants = invariants(start, captured["masses"]) for step in [0.04, 0.02, 0.01]: final = rk4(start, duration, step, captured["masses"], disk) end_invariants = invariants(final, captured["masses"]) reference_checks.append( { "maximumStep": step, "largestCartesianDifferenceFromIas15": max( abs(a - b) for a, b in zip(final, reference) ), "relativeEnergyMinusDiskWork": abs( end_invariants["energy"] - initial_invariants["energy"] - final[12] ) / abs(initial_invariants["energy"]), "relativeAngularMomentumMinusDiskTorque": abs( end_invariants["angularMomentum"] - initial_invariants["angularMomentum"] - final[13] ) / abs(initial_invariants["angularMomentum"]), } ) assert reference_checks[-1]["largestCartesianDifferenceFromIas15"] < 2e-6 assert ( reference_checks[-1]["largestCartesianDifferenceFromIas15"] < reference_checks[0]["largestCartesianDifferenceFromIas15"] / 50 ) local_torque = torque_components(4 * EARTH_SOLAR_MASS_RATIO, 1, 0, DEFAULT_DISK) local_drift = 2 * local_torque["total"] / (4 * EARTH_SOLAR_MASS_RATIO) return { "scope": "44 published measured TTVs; independent torque grid and RK4 check of one synthetic eight-orbit replay. This does not rerun the full12k/50k histories or fit a photodynamical posterior.", "torqueGrid": torque_grid, "oneAuExample": { "massRatio": 4 * EARTH_SOLAR_MASS_RATIO, "driftAuPerCodeTime": local_drift, "driftAuPerJulianYear": local_drift / YEARS_PER_CODE_TIME, "migrationTimeJulianYears": YEARS_PER_CODE_TIME / abs(local_drift), "eccentricityDampingJulianYears": local_torque["eccentricityTime"] * YEARS_PER_CODE_TIME, }, "observedAngles": { "rangesDegrees": observed["rangesDegrees"], "maximumQuarterAlignmentSpanDegrees": observed[ "maximumQuarterAlignmentSpanDegrees" ], "intervalMeaning": observed["intervalMeaning"], }, "syntheticWindows": summaries, "independentRk4Reference": reference_checks, } def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--directory", type=Path, default=Path(__file__).parent) parser.add_argument("--output", type=Path, default=Path("migration-results.json")) args = parser.parse_args() result = reproduce(args.directory) args.output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n") print(json.dumps(result, indent=2)) if __name__ == "__main__": main()