Given a fct F(x, y) , return the ns values of (xs, ys) where
f(xs, ys) = f(xp, yp). Ttot is the total time variable after a close contour
(some codes need ttot to compute the rotational transform along the contour).
Method: Hamilton's equation dx/dt=x*dF/dy; dy/dt=-x*dF/fx are integrated
starting from (xp, yp) until the cummulative angle between (x-x0, y-y0)
and (xnew-x0, ynew-y0) reaches 2*pi, where (xnew, ynew) are new coordinates
along the trajectory.
As a second step, the (xnew, ynew) coordinates are interpolated onto
the equidistant angle grid to yield the (xs, ys).
integer, parameter :: r8 = selected_real_kind(12,100)
integer, intent(in) :: nx, ny
real(r8), intent(in) :: x(nx), y(ny)
real(r8), intent(in) :: f(nx, ny)
real(r8), intent(in) :: x0, y0 approximate contour centre
real(r8), intent(in) :: xp, yp starting point
integer, intent(in) :: ns
real(r8), intent(out) :: ttot, xs(ns), ys(ns)
integer, intent(out) :: ier
integer, parameter :: r8 = selected_real_kind(12,100)