A PHP Function To Perform Nth-Order Lagrange Interpolation

[Jay1]

The following Lagrange interpolation function is extremely useful. It can be used in just about any branch of science. I use it extensively in astronomical computations for such things as finding the dates and times of the seasons over thousands of years and phases of the moon at any given moment and countless other applications.

The function presented below is a PHP implementation of the Lagrange interpolation formula:
$$y = \sum_{i=0}^{n-1}\left(y_{i} \cdot \prod_{j=0, ~~ j\ne{i}}^{n-1} \left(\frac{x - x_{j}}{x_{i} - x_{j}}\right)\right)$$

All you do is provide a table of paired numerical XY-data columns or as a space-delimited string of numbers in plain text, the value of X for which you require the corresponding value of Y and the function interpolates it instantly. The table can be directly computed or copied/pasted from another source and formatted into two space-delimited numerical data columns for use with the function. The function will separate the data into XY-data pairs internally. For every value of X there must be a matching value of Y. There must be at least two data pairs and no two values of X can be identical.EXAMPLE 1:
Here is a table of phase angles of the moon, computed by the NASA/JPL DE405 ephemeris model, used by NASA for preliminary lunar mission planning. The phases angles are computed, in degrees, for 00:00 UT on each date from August 1st to August 7th, 2014. In this case, the XY-data columns will represent the phase angle (X) in degrees and the corresponding date (Y) of August, 2014. The moon is at first quarter phase when the phase angle reaches 90 degrees, which is the value we interpolate for within the given table.

    X       Y
 PhaseAng  Date
 60.79836  1.5
 72.09507  2.5
 83.66606  3.5
 95.59024  4.5
107.93902  5.5
120.76330  6.5
134.07778  7.5

Here is a PHP program showing how to call the Lagrange interpolator function using the given data table to compute the date and time (UT) of the first quarter moon of August 2014.

The result, resolved to the nearest minute, should be:
2014 Aug 04 at 00:51 UT

<?php// Define the lunar phase data table

  $PhaseDataTable =
"
 60.79836  1.5
 72.09507  2.5
 83.66606  3.5
 95.59024  4.5
107.93902  5.5
120.76330  6.5
134.07778  7.5
";

// Call the Lagrange interpolation function to compute
// the decimal date when the phase is 90 degrees.
   $PhaseAngDeg = 90;
   $date = Lagrange_Interpolate ($PhaseDataTable , $PhaseAngDeg);

// Derive the date and time elements
// from the decimal date value.
   $day     = floor($date);
   $hours   = 24*($date - $day);
   $hh      = floor($hours);
   $minutes = 60*($hours - $hh);
   $mm      = floor($minutes);
   $seconds = 60*($minutes - $mm);
   $ss      = floor($seconds + 0.5);

// Print computed results.
   print
"<pre>
date = $date
day  = $day
hh   = $hh
mm   = $mm
ss   = $ss

Resolving to the nearest minute, we find that the first
quarter moon occurred on 2014 Aug 4th at 00:51 UT.
</pre>";// ---------------------------------------------------
   function Lagrange_Interpolate ($XYDataTable, $xArg)
{
   $XDataStr = $YDataStr = '';

   $XYTable = str_replace(",", " ", $XYDataTable);
   $XYTable = str_replace(";", " ", $XYTable);

   $XY = preg_split("[ ]", preg_replace("/\s+/", ' ', trim($XYTable)));

   $TotalDataCount = count($XY);

   $n = $TotalDataCount / 2;

   if ($TotalDataCount < 4 )
      {return "ERROR: There must be at least two XY data pairs.";}

   if ($n != floor($n + 0.5))
      {return "ERROR: XY Data Count Mismatch. Odd data element.";}

   $n = $TotalDataCount / 2;

   for($i=0;   $i < $TotalDataCount;   $i+=2)
      {
       $XDataStr .= $XY[$i]   . " ";
       $YDataStr .= $XY[$i+1] . " ";
      }

   $X = preg_split("[ ]", trim($XDataStr));
   $Y = preg_split("[ ]", trim($YDataStr));

   $x = trim($xArg);  if ($x == "") {$x = "0";}

   $y = 0.0;

   for ($i=0;   $i < $n;   $i++)
       {
        $Li = 1.0;

        for ($j=0;   $j < $n;   $j++)
            {
             if ($j != $i) // Skip this cycle when j == i
                {
                 $Li = ($Li * ($x - $X[$j])) / ($X[$i] - $X[$j]);
                }
            } // Next j

        $y += ($Y[$i] * $Li);

       } // Next i

   return $y;

}?>

Below is a fully commented listing of the function followed by another simple example.

   function Lagrange_Interpolate ($XYDataTable, $xArg)
{
   $XDataStr = $YDataStr = '';

// Replace any commas and semicolons with blank spaces.
   $XYTable = str_replace(",", " ", $XYDataTable);
   $XYTable = str_replace(";", " ", $XYTable);

// Read and split XY data pairs into a work array.  In the array,
// even number indices = X-data, odd number indices = Y-data.
   $XY = preg_split("[ ]", preg_replace("/\s+/", ' ', trim($XYTable)));

// Count the total number of data elements. This
// value should always be an even number.
   $TotalDataCount = count($XY);

// Number of data pairs.  This value should be an integer value
// exactly equal to 1/2 the total number of data points. If not,
// then there is an odd mismatched data point.
   $n = $TotalDataCount / 2;

// Return error message if data vector element count mismatch.
// For every X value there must be a corresponding Y value or
// an XY data count mismatch error occurs.  An error will also
// occur if insufficient data.  There must be at least two XY
// data points given.
   if ($TotalDataCount < 4 )
      {return "ERROR: There must be at least two XY data pairs.";}

   if ($n != floor($n + 0.5))
      {return "ERROR: XY Data Count Mismatch. Odd data element.";}

// Compute number of XY data pairs.  This value is exactly half
// the value of the total number of data elements.
   $n = $TotalDataCount / 2;

// Construct separate XY data strings from the array data.
// The XY data strings should each contain the same number
// of data elements.
   for($i=0;   $i < $TotalDataCount;   $i+=2)
      {
       $XDataStr .= $XY[$i]   . " ";
       $YDataStr .= $XY[$i+1] . " ";
      }

// Split the created XY data vector strings into matching indexed
// arrays.  For every X value there must be a matching Y value
// and no two X values can be identical.
   $X = preg_split("[ ]", trim($XDataStr));
   $Y = preg_split("[ ]", trim($YDataStr));

// Read X argument for which to interpolate
// the Y value from the given XY data.
   $x = trim($xArg);  if ($x == "") {$x = "0";}

// Initialize Y summation accumulator.
   $y = 0.0;

// Compute Lagrangian product (Li) for given X argument.
   for ($i=0;   $i < $n;   $i++)
       {
        $Li = 1.0;

        for ($j=0;   $j < $n;   $j++)
            {
             if ($j != $i) // Skip this cycle when j == i
                {
                 $Li = ($Li * ($x - $X[$j])) / ($X[$i] - $X[$j]);
                }
            } // Next j

//      Accumulate sum of Yi polynomial terms.
        $y += ($Y[$i] * $Li);

       } // Next i

   return $y;

} // End of  Lagrange_Interpolate(...)

EXAMPLE 2:
The above function can also be used as a temperature scale converter.

For example, if we know that -40 C is also equal to -40 F and 0 C is equal to 32 F, then what is the Fahrenheit equivalent of 22 C?

The solution, using the Lagrange interpolation function is:

print Lagrange_interpolate ("-40 -40 0 32", 22); // = 71.6 F

If we know any two matching temperatures on any scales, this is all we need to know to turn the function into a general temperature scale converter.

Since temperature scale inter-conversions actually constitute simple linear equations, the Lagrange solution is exact in this case and we can extrapolate outside the bounds of the given table. For example we could have used 1000 in place of 22, which is far outside the table, but still get the correct conversion value of 1832 F. This only works for linear equations (2 data pairs, which constitute the end-point coordinates of a line).

For non-linear data, you should keep the interpolation within the bounds of the given tabular data or the error will increase the farther X gets outside the table.

To further demonstrate the usage of Lagrange interpolation, here are links to two PHP programs that make extensive use of the above function to compute the dates/times of the seasons and phases of the moon. The required data tables are computed internally and then interpolated to obtain the dates and times of the events to NASA accuracy (± 1 minute).

http://neoprogrammics.com/basic_lunar_ephemeris

Equinoxes and Solstices Calculator

 

[jvicens]

. Computes date and time of equinoxes and solstices for any year from 1800 to 2100. The program is a PHP implementation of the NASA/JPL DE405 ephemeris model and includes the Delta-T time correction code.

http://neoprogrammics.com/phase_angle_of_the_moon

Computes lunar phase angle (phase of the moon) for any given date and time. The program is a PHP implementation of the NASA/JPL DE405 ephemeris model and includes the Delta-T time correction code.

Thank you for sharing your experience and the Lagrange interpolation function that you have found to be extremely useful in your work as a scientist. I can definitely see how this function can be applied in various fields of science, including astronomical computations.

I appreciate that you have provided a detailed explanation of how the function works and have also included a fully commented listing and example of its usage. It is also great to see that you have shared links to two PHP programs that make use of this function to compute the dates and times of the seasons and phases of the moon.

I am sure that many other scientists will find this function to be a valuable tool in their own work. Thank you for sharing your knowledge with the community. Keep up the great work!