Since the routine already traces the closed constant flux contours and can compute the B-field map, it would be very useful to additionally pre-compute and provide flux-surface averaged values for ease of use in downstream applications.
While not an exhaustive list, these quantities from an equivalent OMFIT implementation could be a good starting point:
