DELPHI - How make Ferguson cubic interpolation curve

procedure TForm1.Ferguson;
var
  j, i, iIndex : integer;
  iNumSegmentsCount : integer;
  iT, i1_T, iP0x, iP0y, iP1x, iP1y : double;
  iF1, iF2, iF3, iF4 : double;
  iNewX, iNewY : double;
  pSerieLine : TLineSeries;
  iValue : extended;
begin

  pSerieLine := TLineSeries( Chart1.Series[1] );
  pSerieLine.Clear;

  { draw specified count of the cubic curves }
  for iIndex := 0 to px.Count - 2 do
    begin

       { every cube generate from I parts }
       for i := 0 to giDetail - 1 do
         begin
           iT := i / (giDetail-1);

           { -- calculate Hermid polynomials }

           iF1 := ( 2 * power( iT, 3 ) ) - ( 3 * power( iT, 2 ) ) + 1;
           iF2 := ( -2 * power( iT, 3 ) ) + ( 3 * power( iT, 2 ) );
           iF3 := power( iT, 3 ) - ( 2 * power( iT, 2 ) ) + iT;
           iF4 := power( iT, 3 ) - power( iT, 2 );

           { -- and vectors }

           iP0x := abs( pX.GetValue( iIndex + 1 ) - pX.GetValue( iIndex ) );
           iP0y := abs( pY.GetValue( iIndex + 1 ) - pY.GetValue( iIndex ) );

           iP1x := abs( pX.GetValue( iIndex ) - pX.GetValue( iIndex + 1 ) );
           iP1y := abs( pY.GetValue( iIndex ) - pY.GetValue( iIndex + 1 ) );

           iNewX := ( pX.GetValue( iIndex ) * iF1 ) +
                    ( pX.GetValue( iIndex + 1 ) * iF2 ) +
                    ( iP0x * iF3 ) +
                    ( iP1x * iF4 );
           iNewY := ( pY.GetValue( iIndex ) * iF1 ) +
                    ( pY.GetValue( iIndex + 1 ) * iF2 ) +
                    ( iP0y * iF3 ) +
                    ( iP1y * iF4 );

           pSerieLine.AddXY( iNewX, iNewY );
         end;

  end;
end;