# Rotation period against diameter
# (c) SsODNet, B.Carry & J.Berthier, Lagrange/OCA & LTE/Obs. Paris

import rocks
import numpy as np
import matplotlib.pyplot as plt

# Read SsODNet BFT: only children, diameter and rotation period (for speedup)
bft = rocks.load_bft(columns=["children", "diameter.value", "spins.period.value"])

# Some columns contain `arrays` in the BFT, to access all values you can simply `explode` the DataFrame
df = bft.explode("spins.period.value")


# All Solar System Objects (SSOs)
fig, ax = plt.subplots()
ax.scatter(
    df["diameter.value"],
    df["spins.period.value"],
    s=1,
    alpha=0.2,
    color="grey",
    label="SSOs",
)

# SSO with satellites
cond = df["children"] > 0
ax.scatter(
    df.loc[cond, "diameter.value"],
    df.loc[cond, "spins.period.value"],
    s=4,
    alpha=0.9,
    color="C0",
    label="Binaries",
)

# Valley between Principal-Axis Rotators and Tumblers:
# https://ui.adsabs.harvard.edu/abs/2025NatAs...9..493Z/abstract
xx = np.logspace(-1, 3, num=10)
yy = 13.38 * xx**0.783
ax.plot(xx, yy, color="black", linestyle="--", label="Zhou+2025")

# Spin limit
ax.axhline( 2.2, linestyle='dotted', color='red', label='Spin limit')

# Axes
ax.set_xlabel("Diameter / km")
ax.set_ylabel("Period / h")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(0.1, 3e3)
ax.set_ylim(0.1, 1e3)
ax.invert_yaxis()

leg = ax.legend(loc="upper right")
for lh in leg.legend_handles:
    lh.set_alpha(1)

# To disk
fig.tight_layout()
fig.savefig( 'period-diameter.png', transparent=True, dpi=180)