2D irregular mesh interpolation fortran

Question

How can I implement two-dimensional interpolation in FORTRAN, where the data looks as shown below. x and y are two coordinates, z is a value that depends on them x uniformly, but y is not evenly distributed and the maximum value of y corresponding to uniform values ​​of x continues to increase. Without loss of accuracy -

• What is the easiest way to get the value of z based on a given x and y?
• What is the fastest way to get the value of z based on given x and y?

Thanks SM

xyz ----------- 0 0 - 0 0.014 - 0 0.02 - ..... .... 0.1 0 - 0.1 0.02 - 0.1 0.03 - ....... ..... 1.0 0 - 1.0 0.05 - 1.0 0.08 - ....... ....... 

Answer 1

I'm going to assume that you already read your data in an N x 3 array, following the format you gave. I assume that you do not know in advance what the X-interval is - you definitely do not know the Y-interval, since it is changing. Therefore, I would recommend the following strategy:

• Find out the interval X: start from the first line and go through the elements of X until you see a change in value. Now you know XSTART and XSTEP — you'll need it later.
• Do a binary search in your array for an X value until you find an XFOUND value such that XFOUND <X <XFOUND + XSTART
• Assuming you specify “somewhere in the list”, you find the corresponding value of Y - depending on whether it is greater or less than this value, you increase or decrease the array until you find the first record <Y. The corresponding values ​​are X11, Y11, Z11. The next line of the array has X12 Y12 and Z12.
• You need two more points before you can interpolate - repeat this process, looking for the “next higher value of X”. This will give you XYZ21 and XYZ22
• Now you can think about calculating the interpolated value of Z. In general, there are different methods with different accuracy:
• "Nearest neighbor": find the nearest point and use its Z value (simplest, least accurate)
• linear interpolation: find the three closest points and linearly interpolate values ​​based on relative distances
• “higher order estimates”: for this, you usually need to create a full mapping of network grid points, so you can perform interline spline and get smooth interpolation, which will usually be more accurate at points between grid points (assuming that the function, described by samples is actually a smooth function!)

My FORTRAN is a little rusty - hope this helps.

PS - perhaps a simpler approach is to use the fact that the values ​​of X are already evenly distributed. This allows for more accurate interpolation. See this image:

Answer 2

Before finding the fastest way to solve this problem, I would suggest finding a way to solve this problem. In short, I would suggest:

1) Find the Delaunay triangulation for points (x, y). The Fortran code for this, for example, is in GEOMPACK .

2) For interpolation, given the Delaunay triangulation, find the triangle that contains the point where the interpolation is to be performed, and then interpolate the z value based on the location of the point relative to each of the vertices of the triangle.

(EDIT) I forgot the name of the method (if I ever knew), but thanks to @Floris, a good approach to interpolating in a triangle is called barycentric interpolation , which finds the interpolated value based on the ratio of the areas in the three smaller triangles into which you can split the big triangle by drawing lines from a point inside a large triangle to each of the corners of the triangle. The area of ​​each small triangle can be found along the length of each side of the triangle according to the Heron formula .

If speed improvements were needed, I think that they could be obtained mainly by finding a quick way to find a triangle that contains an interpolation point, but I would like to profile some test runs before choosing which bits of code to optimize.