geom_model.py

import numpy as np
import matplotlib.pyplot as plt

def compensatory_radius(estimate_dispersion, R):
    cosarg = estimate_dispersion/R**2
    # Cos will go negative when arg > pi/2. arg < 0 makes no sense, so no need to worry about that.
    # Trigonometry breaks down for this at pi/2? Could probably have some sort of analytical
    # continuation if needed, however pi/2 dispersion leading to center seems implausible for a
    # model.
    r_2 = R*np.cos(cosarg)
    r_2[cosarg > np.pi/2] = np.nan
    
    return r_2

# Just some radius. I think I (JP) got an analytical result that the r_2 is independent
# from R when minimizing (squared?) distance when assuming Gaussian/Von-Mises distribution
# for the phase uncertainty. Seems off that this is dependent on units (R^2 in the formula makes
# this unit-dependent I think).
R = 13

ts = np.arange(0, 30, 0.01)
alphas = np.linspace(0, 10, 3)

for alpha in alphas:
    phase_dispersion = ts*alpha
    r_2 = compensatory_radius(phase_dispersion, R)
    plt.plot(ts, r_2)

plt.show()