Constructing magnetic surface coordinate system from discrete $\Psi (R, Z)$ data

Given an axisymmetric tokamak equilibrium in $(R, \phi, Z)$ coordinates (e.g., 2D data $\Psi (R, Z)$ on a rectangular grids $(R, Z)$ in G-file), we can construct a magnetic surface coordinates $(\psi , \theta , \phi )$ by the following two steps. (1) Find out a series of magnetic surfaces on $(R, Z)$ plane and select radial coordinates for each magnetic surface (e.g. the poloidal flux within each magnetic surface). (2) Specify the Jacobian or some property you want the poloidal angle to have. Then calculate the poloidal angle of each point on each flux surface (on the $\phi = \ensuremath{\operatorname{const}}$ plane) by using Eq. (188) (if the Jacobian is specified) or some method specified by you to achieve some property you prefer for the poloidal angle (if you do not specify a Jacobian). Then we obtain the magnetic surface coordinates system $(\psi , \theta , \phi )$.