unit GaussUnit;

interface

uses TypesUnit;

type

TGauss = class
  private
    Dis: Double;
    Sigma: Double;
    px: TFloatArray;
    py: TFloatArray;
  public
    constructor Create(Ds, Sigm: Double);
    property x: TFloatArray read px;
    property y: TFloatArray read py;
    destructor Destroy; override;
  end;

TGauss1D = class (TGauss)
  private
    pIntegralSum: TFloatArray;
    pIntegralNorm: Double;
  public
    constructor Create(Ds, Sigm: Double);
    property IntegralSum: TFloatArray read pIntegralSum;
    property IntegralNorm: Double read pIntegralNorm;
    destructor Destroy; override;
  end;



implementation

constructor TGauss.Create(Ds, Sigm: Double);
var i, n: Integer;
    N2Sigma2, h_2, Mid: Double;
begin
  inherited Create;
  Dis:=Ds;
  Sigma:=Sigm;

  n:=RoundUp(Dis/cGaussStep);
  SetLength(px, n+1);
  SetLength(py, n+1);

  N2Sigma2:=2*Sigma*Sigma;
  h_2:=cGaussStep/2;

  px[0]:=0;
  py[0]:=0;
  for i:=1 to n-1 do begin
    px[i]:=px[i-1]+cGaussStep;
    Mid:=px[i]-h_2;
    py[i]:=exp(-Mid*Mid/N2Sigma2);
  end;
  px[n]:=Dis;
  Mid:=(px[n]+px[n-1])/2;
  py[n]:=exp(-Mid*Mid/N2Sigma2);
end;



destructor TGauss.Destroy;
begin
  px:=nil;
  py:=nil;
  inherited
end;



constructor TGauss1D.Create(Ds, Sigm: Double);
var i, n: Integer;
begin
  inherited Create(Ds, Sigm);
  n:=RoundUp(Dis/cGaussStep);
  SetLength(pIntegralSum, n+1);

  pIntegralSum[0]:=0;
  for i:=1 to n-1 do
    pIntegralSum[i]:=pIntegralSum[i-1]+cGaussStep*py[i];
  pIntegralSum[n]:=pIntegralSum[n-1]+(px[n]-px[n-1])*py[n];
  pIntegralNorm:=pIntegralSum[n];
end;



destructor TGauss1D.Destroy;
begin
  px:=nil;
  py:=nil;
  pIntegralSum:=nil;
  inherited;
end;

end.
