I have a [data file][1] containing the four orbital element $(x,y,p_x,p_y)$ of a two-dimensional closed loop orbit. Note that the exact implicit equation of the orbit is not known. 

The plot of the 2D orbit on the $(x,y)$ plane

    data = Import["lyap_4d.dat", "Table"];
    d00 = data[[All, {1, 2}]];
     
    C0 = ListLinePlot[d00, AspectRatio -> 1, PlotRange -> All, 
         Frame -> True, Axes -> False, PlotStyle -> {Black, Thick}]

[![enter image description here][2]][2]

Now let's define a rectangular grid of initial conditions around this orbit

    data = Flatten[Table[{i, j}, {i, 0.785, 1.135, 0.005}, 
           {j, -0.45, 0.45, 0.005}], 1];
    nic = Length[data]
    L0 = ListPlot[data, PlotStyle -> {Blue, PointSize[0.001]}];
    plot1 = Show[{L0, C0}]

[![enter image description here][3]][3]

Finally, we determine which of the initial conditions are inside the orbit

    poly = Cases[Normal@C0, Line[x_] :> x, Infinity];
    inPolyQ[poly_, pt_] := Graphics`Mesh`PointWindingNumber[poly, pt] =!= 0
    data2 = Select[data, inPolyQ[poly[[1]], #] &];
    nic2 = Length[data2]
    L1 = ListPlot[data2, PlotStyle -> {Red, PointSize[0.001]}];
    plot2 = Show[{L1, C0}, Frame -> True, Axes -> False]

[![enter image description here][4]][4]

The list `data2` contains all the $(x_0,y_0)$ initial conditions inside the orbit. 

My question is the following: For the boundary orbit the $p_x$ and $p_y$ values are known. How can I use these values so as to predict in a way the $(p_{x0},p_{y0})$ of the orbits with initial conditions inside the orbit? In other words, how can I interpolate the $p_x$ and $p_y$ of the orbit by taking into account the values of the boundary?

I use version 9.0 of Mathematica in Win XP SP3. 

Many thanks in advance! 

  [1]: http://www.mediafire.com/download/cf1f0yz0jmdu1da/lyap_4d.dat
  [2]: https://i.sstatic.net/lMcQA.png
  [3]: https://i.sstatic.net/r5YHx.png
  [4]: https://i.sstatic.net/XhEWO.png