DELPHI - How to show chi-square distribution curve

Propability density function (chi-square distribution).

var
  i, iK : integer;
  iX, iY, iGamma, iNumerator, iDenominator : double;
  pSerie : TLineSeries;
  pIntegral : TIntegral;
begin
  Memo1.Clear;

  pSerie := TLineSeries( Chart1.Series[0] );
  pSerie.Clear;

  pIntegral := TIntegral.Create;
  pIntegral.iStep := 0.05;
  pIntegral.SetInterval( 0.001, 20 );

  { -- set param + calc GAMMA function }

  iK := 8;
  iGamma := GetGamma( iK/2 );

  for i := 0 to pIntegral.pX.Count - 1 do
    begin
      iX := pIntegral.pX.GetValue( i );

      try
        iNumerator := power( iX, iK/2-1 ) * exp( ( -1*iX ) /2 );
        iDenominator := power( 2, iK / 2 ) * iGamma;

        if iDenominator <> 0 then
          iY := iNumerator / iDenominator
        else
          iY := 0;
      except
        iY := 0;
      end;
      pIntegral.pY.AddValue( iY );

      { -- add to chart }

      pSerie.AddXY( iX, iY );
    end;
end;

And here is very important Gamma() function with k parameter.

function GammaStirF( X : double ) : double;
var
  y : double;
  w : double;
  v : double;
  stir : double;
begin
  w := 1 / x;
  stir := 7.87311395793093628397E-4;
  stir := -2.29549961613378126380E-4 + w * stir;
  stir := -2.68132617805781232825E-3 + w * stir;
  stir := 3.47222221605458667310E-3 + w * stir;
  stir := 8.33333333333482257126E-2 + w * stir;
  w := 1 + w * stir;
  y := exp(x);
  if x > 143.01608 then
    begin
      v := power( x, 0.5 * x -0.25 );
      y := v *( v / y );
    end
  else
    begin
      y := power( x, x -0.5 ) / y;
    end;
  result := 2.50662827463100050242 * y * w;
end;

function GetGamma( _x : double ) : double;
var
  p : double;
  PP : double;
  q : double;
  QQ : double;
  z : double;
  i : longint;
  SgnGam : double;
begin
  SgnGam := 1;
  q := abs( _x );
  if q > 33.0 then
  begin
    if _x < 0.0 then
    begin
      p := floor(q);
      i := round(p);

      if i mod 2 = 0 then
      begin
        SgnGam := -1;
      end;

      z := q - p;

      if z > 0.5 then
      begin
        p := p+1;
        z := q-p;
      end;

      z := q * sin( pi * z );
      z := Abs(z);
      z := pi / ( z * GammaStirF( q ) );
    end
    else
    begin
      z := GammaStirF( _x );
    end;
    result := SgnGam * z;
    exit;
  end;
  z := 1;
  while _x >= 3 do
  begin
    _x := _x - 1;
    z := z *_x;
  end;
  while _x < 0 do
  begin
    if _x > -0.000000001 then
    begin
      result := z / ( ( 1 + 0.5772156649015329 *_x ) * _x );
      exit;
    end;
    z := z /_x;
    _x := _x + 1;
  end;
  while _x < 2 do
  begin
    if _x < 0.000000001 then
    begin
      result := z / ( ( 1 + 0.5772156649015329 * _x ) *_x );
      exit;
    end;
    z := z /_x;
    _x := _x + 1.0;
  end;
  if _x = 2 then
  begin
    result := z;
    exit;
  end;

  _x := _x - 2.0;
  PP := 1.60119522476751861407E-4;
  PP := 1.19135147006586384913E-3 + _X * PP;
  PP := 1.04213797561761569935E-2 + _X * PP;
  PP := 4.76367800457137231464E-2 + _X * PP;
  PP := 2.07448227648435975150E-1 + _X * PP;
  PP := 4.94214826801497100753E-1 + _X * PP;
  PP := 9.99999999999999996796E-1 + _X * PP;
  QQ := -2.31581873324120129819E-5;
  QQ := 5.39605580493303397842E-4 + _X * QQ;
  QQ := -4.45641913851797240494E-3 + _X * QQ;
  QQ := 1.18139785222060435552E-2 + _X * QQ;
  QQ := 3.58236398605498653373E-2 + _X * QQ;
  QQ := -2.34591795718243348568E-1 + _X * QQ;
  QQ := 7.14304917030273074085E-2 + _X * QQ;
  QQ := 1.00000000000000000320 + _X * QQ;

  result := z * PP / QQ;
  exit;
end;

Output:

k=2

k=8