When doing field line integration from a magnetic field, would be useful to also get out the iota values. Multiple ways to do this but the easiest to implement initially is probbaly to just do the literal definition for iota,

$$ \iota = \lim_{n\rightarrow \infty} \frac{1}{n} \sum_1^n \Delta \theta_n$$

where $\Delta \theta_n$ is the change in poloidal angle in the $\phi=0$ plane from the $n-1$ toroidal transit to the $n$ th toroidal transit, which could be calculated as $arctan(\frac{Z-Z_{axis}}{R-R_{axis}})$

Where iota is calculated in FOCUS, as an example

https://github.com/PrincetonUniversity/FOCUS/blob/be28c104bebc1b4ea281547b03fc1e9fab34e8e2/sources/poinplot.f90#L113