hdf5_to_reb.py

import h5py
import numpy as np
import rebound

def hdf5_to_rebound(row):
    row = row.reshape((-1,5))
    Pmin = np.min(row[:,1])
    TCmin = np.min(row[:,4])
    epoch = TCmin - 0.1 * Pmin
    
    sim = rebound.Simulation()
    sim.add(m=1)
    
    for m, P, h, k, Tc in row:
        e, pomega = np.sqrt(h*h + k*k), np.arctan2(k,h)
        l0 = np.mod(2.0*np.pi*(epoch - Tc)/P, 2*np.pi) + 0.5*np.pi
        inc = 1e-3*np.random.normal() # small random inclination
        P = P/365.2569 * np.pi*2  # days to year/2pi
        sim.add(m=m, P=P, e=e, pomega=pomega, M=l0-pomega, inc=inc, Omega="uniform", primary=sim.particles[0])

    sim.move_to_com()
    return sim

h5data = h5py.File("NBody_MCMC_Posteriors.hdf5")
mcmc_data = h5data["Kepler-11"]['DefaultPriors']['PosteriorSample']

mcmc_sample = mcmc_data[np.random.randint(0,len(mcmc_data))] # random sample

sim = hdf5_to_rebound(mcmc_sample)
sim.status()