# contour.f90

subroutine i2mex_contour (nx, ny, x, y, f, x0, y0, xp, yp, ns, ttot, xs, ys, ier)

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

subroutine i2mex_xydot (neq, t, y, ydot)

integer, parameter :: r8 = selected_real_kind(12,100)

subroutine i2mex_jac