Given an axisymmetric tokamak equilibrium in (R,ϕ,Z) coordinates (e.g., 2D data Ψ(R,Z) on a rectangular grids (R,Z) in G-file), we can construct a magnetic surface coordinates (ψ,𝜃,ϕ) 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 that we want the poloidal angle to have. Then calculate the poloidal angle of each point on each flux surface (on the ϕ = const plane) by using Eq. (168) (if the Jacobian is specified) or some method specified by us to achieve some property we prefer for the poloidal angle (if a Jacobian is not directly specified). Then we obtain the magnetic surface coordinates system (ψ,𝜃,ϕ).