I've been working on a program which calculates, given a point and 4 surrounding points, the [Lagrange polynomial][1], in order to interpolate a value. Consider that I'm not a mathematician and I better understand code than formulas.

I've came up with the following code, which works, but I really don't think is general (and it's pretty ugly to me).

    using System;
    using System.Collections.Generic;
    using System.Linq;

    namespace LagrangeInterpolation
    {
        public class Point : ICloneable
        {
            public double X { get; set; }

            public double Y { get; set; }

            public double Value { get; set; }   

            public Point Clone()
            {
                return (Point)this.MemberwiseClone();
            }           
        }

        public static class Lagrange
        {
            public static Point Interpolate(Point[] controlPoints, Point point)
            {           
                var A = -controlPoints[0].X - controlPoints[1].X + controlPoints[2].X + controlPoints[3].X;
                var B = -controlPoints[0].X + controlPoints[1].X + controlPoints[2].X - controlPoints[3].X;
                var C = +controlPoints[0].X - controlPoints[1].X + controlPoints[2].X - controlPoints[3].X;
                var X = 4 * point.X - controlPoints.Sum(x => x.X);

                var D = -controlPoints[0].Y - controlPoints[1].Y + controlPoints[2].Y + controlPoints[3].Y;
                var E = -controlPoints[0].Y + controlPoints[1].Y + controlPoints[2].Y - controlPoints[3].Y;
                var F = +controlPoints[0].Y - controlPoints[1].Y + controlPoints[2].Y - controlPoints[3].Y;
                var Y = 4 * point.Y - controlPoints.Sum(x => x.Y);

                var r = (X / 4 - B * Y / 4 * E) / (1 - D / 4 * E);
                var s = (Y - D * r) / E;

                var prevR = 0d;
                var prevS = 0d;
                const double precision = 0.00000001;
                while (!(prevR - r < precision && prevS - s < precision))
                {
                    prevR = r;
                    prevS = s;
                    r = (X - B * s) / (A + C * s);
                    s = (Y - D * r) / (E + F * r);
                }

                // Interpolate value
                var result = point.Clone();
                result.Value = ((1 - r) * (1 - s) * controlPoints[0].Value + (1 - r) * (1 + s) * controlPoints[1].Value + (1 + r) * (1 + s) * controlPoints[2].Value + (1 + r) * (1 - s) * controlPoints[3].Value) / 4;
                            
                return result;
            }
        }
    }

The input parameters are:
 
 - `controlPoints`: the 4 points, each one with its `Value`. 
 - `point`: the point for which we want to calculate the interpolated value

The returned `Point` is a clone of the Input point with the Value property set.
Every instance of `Point` have `X` and `Y` normalised within range -1..1 (I subtract the quadrilateral center from each point).

Example:

	controlPoints = new [] {
		new Point()	{ X = -0.033675000000000566, Y = -0.02564999999999884, Value = 1.2787 },
		new Point()	{ X = -0.035524999999999807, Y = 0.024329999999999075, Value = 1.329 },
		new Point()	{ X = 0.03370499999999943, Y = 0.02564999999999884, Value = 1.3376 },
		new Point()	{ X = 0.035494999999999166, Y = -0.024329999999999075, Value = 1.302 }
	}

	point = new Point() { X = 0.018148174616284152, Y = -0.014201699949808244 }


Expected result is `Point.Value = 1.3044829106888913`

Can someone suggest a better way (formally and mathematically) to perform this calculation?


  [1]: https://en.wikipedia.org/wiki/Lagrange_polynomial