i1 = np.abs(w0.v_z.value - 100).argmin()  # left of res
i2 = np.abs(w0.v_z.value - 118).argmin()  # in res
i3 = np.abs(w0.v_z.value - 155).argmin()  # right of res

for j, i in enumerate([i1, i2, i3]):
    orbits[:, i].cylindrical.plot(
        ['rho', 'z'], alpha=0.5, marker=',', axes=[axes[j]]
    )