(********************************************
  Version: 7.1
  Date: 2007/02/01
  This is shareware.
  Copyright (C) Michal Pohanka
  All rights reserved.

  Internet: http://www.volny.cz/profipohanka
 ********************************************)

unit UPohMath;

interface

uses
  Math, Types, UPohTypes;

type
  TOptimParab = record
    a, b: double; // intital range
    x, fx: double; // zde se uklada vysledek
    d, u, v, w: double;
    fu, fv, fw: double; // function values
    e: double;
    AccuracyFrac, ZeroTreshold: double; // max odchylka je AccuracyFrac*Abs(x)+ZeroTreshold, pokud je mensi je optimalizace u konce
  end;

procedure SwapDouble(var a, b: double);
procedure SwapInt(var a, b: integer);

function PohZero(a: single): boolean; overload;
function PohZero(a: double): boolean; overload;
function PohZero(a: extended): boolean; overload;

function LessInt(x: double): integer;

function FitImageWithinRect(var ARect: TRect; ImageWidth, ImageHeight, RectWidth, RectHeight: integer): boolean; overload; // centers the image in the rect
function FitImageWithinRect(var left, top, width, height: integer; ImageWidth, ImageHeight, RectWidth, RectHeight: integer): boolean; overload; // centers the image in the rect
function FitImageWithinRect(var left, top, width, height: double; ImageWidth, ImageHeight, RectWidth, RectHeight: double): boolean; overload; // centers the image in the rect
function FixRectWtoH(r: TRect; WtoH: double): TRect; // shrinks the rect so that the returned rect has width/height=WtoH

function PohSign(a: double): TValueSign; overload; // returns -1, 0, 1
function PohSign(a: extended): TValueSign; overload; // returns -1, 0, 1

function IsEqualSingle(a, b: single; dig: integer = 6): boolean; // dig is number of left most digits to compare
function IsEqual(a, b: double; dig: integer = 11): boolean; overload;// dig is number of left most digits to compare
function IsEqualExt(a, b: extended; dig: integer = 15): boolean; // nemuze byt overloaded
function IsEqual(a, b: TXYSingle; dig: integer = 6): boolean; overload;// dig is number of left most digits to compare
function IsEqual(a, b: TXYDouble; dig: integer = 11): boolean; overload;// dig is number of left most digits to compare
function IsEqual(a, b: TXYZDouble; dig: integer = 11): boolean; overload;// dig is number of left most digits to compare

function PohDist(x1, y1: double): double; overload;
function PohDist(x1, y1: extended): extended; overload;
function PohDist(a, b: TXYSingle): single; overload;
function PohDist(a, b: TXYDouble): double; overload;
function PohDist(a, b: TXYZDouble): double; overload;
function PohDist(x1, y1, x2, y2: double): double; overload;
function PohDist(x1, y1, x2, y2: extended): extended; overload;
function PohDist(x1, y1, z1, x2, y2, z2: double): double; overload;

procedure GetEllipseRange(var wx, wy, hx, hy: double; a, b, alphaRad: double); // [wx, wy] the most right point for Y down, [hx, hy] the most top point for Y down; axis radius (a, b) must be >0; alpha is ellipse CCW rotation for Y down
function GetEllipseVolumePart(a, b, x1, x2: double): double; // axis radius (a, b) must be >0; x1 a x2 is any number
procedure RotatePointDeg(var DestX, DestY: double; OrigX, OrigY, Alpha: double); overload; // clockwise for Y up, CCW for Y down
procedure RotatePointRad(var DestX, DestY: double; OrigX, OrigY, Alpha: double); overload; // clockwise for Y up, CCW for Y down
procedure RotatePointDeg(var DestX, DestY: double; OrigX, OrigY, Alpha, CenterX, CenterY: double); overload; // clockwise for Y up, CCW for Y down
procedure RotatePointRad(var DestX, DestY: double; OrigX, OrigY, Alpha, CenterX, CenterY: double); overload; // clockwise for Y up, CCW for Y down

function CircleLineIntersect(var x1, y1, x2, y2: double; r: double; y0, q: double): boolean; // circle: x^2+y^2=r; line:=y=y0+q*x;
function CircleLineTanIntersect(var x1, y1, x2, y2: double; r: double; x0, y0, alphaRad: double): boolean; // circle: x^2+y^2=r; line: primka prochazejici bodem [x0, y0] pod uhlem alfa od osy X; CCW for Y up
function CircleLineCotIntersect(var x1, y1, x2, y2: double; r: double; x0, y0, alphaRad: double): boolean; // circle: x^2+y^2=r; line: primka prochazejici bodem [x0, y0] pod uhlem alfa od osy Y; CCW for Y up
function Circle2Intersect(var x1, y1, x2, y2: double; cx1, cy1, r1, cx2, cy2, r2: double): boolean; // circle: (x-cx1)^2+(y-cy)^2=r

function LineIntersect(var x, y: double; ax1, ay1, ax2, ay2: double; bx1, by1, bx2, by2: double): boolean;
procedure PointOnLine(var p: TXYDouble; p1, p2: TXYDouble; GivenPoint: char); overload; // GivenPoint='X', 'Y'  and is stored in p
procedure PointOnLine(var p: TXYZDouble; p1, p2: TXYZDouble; GivenPoint: string); overload; // GivenPoint='X', 'Y', 'Z', 'XY', 'XZ', 'YZ' and is stored in p
function PointInLineSegment(p1x, p1y, p2x, p2y, testPointx, testPointy: double; precision: double=1e-10): boolean; overload; // tests ranges only, not position on line
function PointInLineSegment(p1, p2, testPoint: TXYZDouble; precision: double=1e-10): boolean; overload; // tests ranges only, not position on line
function PointInRect(x, y, x1, y1, x2, y2: double): boolean;

function LinInterp(x, x0, y0, x1, y1: double): double; // !!!x0<>x1!!! ; returns y=(y1-y0)/(x1-x0)*(x-x0)+y0;
function LinInterpSafe(x, x0, y0, x1, y1: double): double; // returns (y1-y0)/(x1-x0)*(x-x0)+y0; if x1=x0 then average is returned (y0+y1)/2
function ParabInterp(x, x0, y0, x1, y1, x2, y2: double): double; overload; // !!!
function ParabInterp(var a, b, c: double; x0, y0, x1, y1, x2, y2: double): double; overload; // !!!
function ParabInterp(x, a, b, c: double): double; overload; // !!!
function QuadrInterp(x, y, x1, y1, x2, y2, v11, v21, v12, v22: double): double; // vXY
function QuadrInterpSafe(x, y, x1, y1, x2, y2, v11, v21, v12, v22: double): double; // vXY if x1=x2 or y1=y2 then average is returned
function GetAngleRad(x1, y1, x2, y2: double): double; // 0rad = 0 clock, Pi=6, 2Pi rad = 12 clock; CW for x and y like display coordinates, CCW for Y up
function GetAngleDeg(x1, y1, x2, y2: double): double; // 0o = 0 clock, 90=6, 360o = 12 clock; CW for x and y like display coordinates, CCW for Y up
function HalfCycleRad(alpha: double): double; // returns alpha within range -Pi..Pi
function HalfCycleDeg(alpha: double): double; // returns alpha within range -180..180
function OneCycleRad(alpha: double): double; // returns alpha within range -Pi..Pi
function OneCycleDeg(alpha: double): double; // returns alpha within range -180..180

function Tridiag(var a, b, c, x, r: TDoubleArray): boolean; // Solves for a vector x[0..n-1] the tridiagonal linear set given by a[0..n-1], b[0..n-1], c[0..n-1], and "right side"=r[0..n-1] are input vectors and are not modied.
// b1 c1 0  0 ... 0   x1   r1
// a2 b2 c1 0 ... 0 * x2 = r2

procedure FFT(var data: TArrXYDouble; inverse: boolean=false); // Number of samples must be positive integer power of 2
procedure CreateCubicSpline(var points, line: TArrXYDouble; Xmin, Xmax: double; throughPoints: boolean=true; lineStep: double=1; minDist: double=1); overload; // Xmin=points[0].x; Xmax=points[High{points)].x; minDist=lineStep
procedure CreateCubicSpline(var points, line: TArrXYInteger; Xmin, Xmax: double; throughPoints: boolean=true; lineStep: double=1; minDist: double=1); overload; // Xmin=points[0].x; Xmax=points[High{points)].x; minDist=lineStep

procedure GaussFiltr(var dataIn, dataOut: TSingleArray; sigma: single); overload; // filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
procedure GaussFiltr(var data: TSingleArray; sigma: single); overload; // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))
procedure GaussFiltr(var dataIn, dataOut: TDoubleArray; sigma: double); overload; // filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
procedure GaussFiltr(var data: TDoubleArray; sigma: double); overload; // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))
procedure GaussFiltrAdapt(var dataIn, dataOut, sigma: TSingleArray); overload; // filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
procedure GaussFiltrAdapt(var data, sigma: TSingleArray); overload; // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))
procedure GaussFiltrAdapt(var dataIn, dataOut, sigma: TDoubleArray); overload; // filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
procedure GaussFiltrAdapt(var data, sigma: TDoubleArray); overload; // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))

function T_K(mV, T0: double): double;

procedure OptimParabInit(var values: TOptimParab; a, b, c, fc: double; AccuracyFrac: double = 1e-6; ZeroTreshold: double = 1e-6); // c must be within a and b, fc is function value in c; max odchylka je AccuracyFrac*Abs(x)+ZeroTreshold, pokud je mensi je optimalizace u konce
function OptimParab1step(var values: TOptimParab): boolean; // true if done; after this function compute fu for u
function OptimParab2step(var values: TOptimParab): boolean; // true if done

procedure IdentityMatrix(var a: TDoubleMatrixAlgn; n: integer); overload; // allocates 'a' if necessary

function isPow2(x: cardinal): boolean;
function Pow2Low(x: cardinal): cardinal;

implementation

uses
  SysUtils;

function ParabInterp(x, x0, y0, x1, y1, x2, y2: double): double; overload; // !!!
var
  a, b, c, d, x02, x12, x22: double;
begin
  x02:=sqr(x0);
  x12:=sqr(x1);
  x22:=sqr(x2);
  a:=-((x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2)/((x0-x1)*x22+(x1-x2)*x02+(x2-x0)*x12);
  d:=((x0-x1)*(x22-x0*x2+(x0-x2)*x1));
  b:=((y1-y2)*x02+(y0-y1)*x22+(y2-y0)*x12)/d;
  c:=((x1*y2-x2*y1)*x02+(x2*y0-x0*y2)*x12+(y1*x0-y0*x1)*x22)/d;
end;

function ParabInterp(var a, b, c: double; x0, y0, x1, y1, x2, y2: double): double; overload; // !!!
var
  d, x02, x12, x22: double;
begin
  x02:=sqr(x0);
  x12:=sqr(x1);
  x22:=sqr(x2);
  a:=-((x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2)/((x0-x1)*x22+(x1-x2)*x02+(x2-x0)*x12);
  d:=((x0-x1)*(x22-x0*x2+(x0-x2)*x1));
  b:=((y1-y2)*x02+(y0-y1)*x22+(y2-y0)*x12)/d;
  c:=((x1*y2-x2*y1)*x02+(x2*y0-x0*y2)*x12+(y1*x0-y0*x1)*x22)/d;
end;

function ParabInterp(x, a, b, c: double): double; overload; // !!!
begin
  result:=a*sqr(x)+b*x+c;
end;

function isPow2(x: cardinal): boolean;
begin
  result:=(x and (x-1))=0;
end;

function Pow2Low(x: cardinal): cardinal;
var
  i, n: integer;
begin
  n:=0;
  while (x>1) do begin
    Inc(n);
    x:=x shr 1;
  end;
  result:=1;
  for i:=1 to n do result:=result shl 1;
end;

procedure CreateCubicSpline(var points, line: TArrXYInteger; Xmin, Xmax: double; throughPoints: boolean=true; lineStep: double=1; minDist: double=1); overload; // Xmin=points[0].x; Xmax=points[High{points)].x; minDist=lineStep
var
  i, n: integer;
  apoints, aline: TArrXYDouble;
begin
  n:=Length(points);
  SetLength(apoints, n);
  for i:=0 to n-1 do with apoints[i] do begin
    x:=points[i].x;
    y:=points[i].y;
  end;
  CreateCubicSpline(apoints, aline, Xmin, Xmax, throughPoints, lineStep, minDist);
  n:=Length(aline);
  SetLength(line, n);
  for i:=0 to n-1 do with line[i] do begin
    x:=Round(aline[i].x);
    y:=Round(aline[i].y);
  end;
  SetLength(apoints, 0);
  SetLength(aline, 0);
end;

procedure CreateCubicSpline(var points, line: TArrXYDouble; Xmin, Xmax: double; throughPoints: boolean=true; lineStep: double=1; minDist: double=1); // Xmin=points[0].x; Xmax=points[High{points)].x; minDist=lineStep
var
  i, n, nl: integer;
  k, kk: integer;
  xk, xk1: double;
  h: TDoubleArray;
  a, b, c, M, r: TDoubleArray;
  x, d, x1, x2, x3, x4, y1, y2, y3, y4, aa, bb, cc, dd: double;
  linearne: boolean;
begin
  try
    n:=Length(points);
    nl:=Round((Xmax-Xmin)/lineStep)+100;
    SetLength(line, nl);
    if Xmin>Xmax then begin
      d:=Xmin;
      Xmin:=Xmax;
      Xmax:=d;
    end;
    if n<1 then begin // nic=0
      // vytvorit body
      i:=0;
      x:=Xmin;
      while x<=Xmax do begin
        if i>=nl then begin
          nl:=nl+100;
          SetLength(line, nl);
        end;
        line[i].x:=x;
        line[i].y:=0;
        Inc(i);
        x:=x+lineStep;
        if throughPoints and (x>Xmax) and (i>0) and not IsEqual(line[i-1].x, Xmax) then x:=Xmax;
      end;
      SetLength(line, i);
    end
    else if n=1 then begin // konstanta
      // vytvorit body
      i:=0;
      x:=Xmin;
      k:=0;
      while x<=Xmax do begin
        if i>=nl then begin
          nl:=nl+100;
          SetLength(line, nl);
        end;
        line[i].x:=x;
        line[i].y:=points[0].y;
        Inc(i);
        x:=x+lineStep;
        if throughPoints and (k=0) and (points[k].x<x) then begin
          x:=points[0].x;
          Inc(k);
        end;
        if throughPoints and (x>Xmax) and (i>0) and not IsEqual(line[i-1].x, Xmax) then x:=Xmax;
      end;
      SetLength(line, i);
    end
    else if n=2 then begin // linearni interpolace
      x1:=points[0].x;
      y1:=points[0].y;
      x2:=points[1].x;
      y2:=points[1].y;
      // vytvorit body
      i:=0;
      x:=Xmin;
      k:=0;
      while (k<n) and (points[k].x<x) do Inc(k);
      while x<=Xmax do begin
        if i>=nl then begin
          nl:=nl+100;
          SetLength(line, nl);
        end;
        line[i].x:=x;
        line[i].y:=LinInterpSafe(x, x1, y1, x2, y2);
        Inc(i);
        x:=x+lineStep;
        while throughPoints and (k<n) and (points[k].x<x) do begin
          x:=points[k].x;
          Inc(k);
        end;
        if throughPoints and (x>Xmax) and (i>0) and not IsEqual(line[i-1].x, Xmax) then x:=Xmax;
      end;
      SetLength(line, i);
    end
    else if n=3 then begin // kvadraticka interpolace
      x1:=points[0].x;
      y1:=points[0].y;
      x2:=points[1].x;
      y2:=points[1].y;
      x3:=points[2].x;
      y3:=points[2].y;
      if (Abs(x1-x3)<1e-100) or (Abs((x1-x3)*(-x1*x3+x1*x2+x2*x3-Sqr(x2)))<1e-100)
      then linearne:=true
      else begin
        linearne:=false;
        cc:=(x1*y3-x3*y1-x2*y3+x2*y1+x3*y2-x1*y2)/((x1-x3)*(-x1*x3+x1*x2+x2*x3-Sqr(x2)));
        bb:=(y1-y3+cc*Sqr(x3)-cc*Sqr(x1))/(x1-x3);
        aa:=y3-bb*x3-cc*Sqr(x3);
      end;
      // vytvorit body
      i:=0;
      x:=Xmin;
      k:=0;
      while (k<n) and (points[k].x<x) do Inc(k);
      while x<=Xmax do begin
        if i>=nl then begin
          nl:=nl+100;
          SetLength(line, nl);
        end;
        line[i].x:=x;
        if linearne
        then line[i].y:=LinInterpSafe(x, x1, y1, x3, y3)
        else line[i].y:=aa+bb*x+cc*Sqr(x);
        Inc(i);
        x:=x+lineStep;
        while throughPoints and (k<n) and (points[k].x<x) do begin
          x:=points[k].x;
          Inc(k);
        end;
        if throughPoints and (x>Xmax) and (i>0) and not IsEqual(line[i-1].x, Xmax) then x:=Xmax;
      end;
      SetLength(line, i);
    end
    else if n=4 then begin // polynom 3 interpolace
      x1:=points[0].x;
      y1:=points[0].y;
      x2:=points[1].x;
      y2:=points[1].y;
      x3:=points[2].x;
      y3:=points[2].y;
      x4:=points[3].x;
      y4:=points[3].y;
      if (Abs(x4-x3)<1e-100) or (Abs(x4*x2-x4*x3+x2*x3-Sqr(x2))<1e-100) or (Abs(x3*x4*x2-x4*x1*x2+x4*Sqr(x1)-x4*x1*x3+Sqr(x1)*x2-x2*x1*x3-Sqr(x1)*x1+Sqr(x1)*x3)<1e-100)
         or (Abs((x4-x3)*(x4*x2-x4*x3+x2*x3-Sqr(x2)))<1e-100)
         or (Abs((x4-x3)*(x4*x2-x4*x3+x2*x3-Sqr(x2))*(x3*x4*x2-x4*x1*x2+x4*Sqr(x1)-x4*x1*x3+Sqr(x1)*x2-x2*x1*x3-Sqr(x1)*x1+Sqr(x1)*x3))<1e-100)
      then linearne:=true
      else begin
        linearne:=false;
        dd:=(Sqr(x4)*x1*y2-Sqr(x4)*x2*y1-Sqr(x4)*x3*y2-Sqr(x4)*x1*y3+Sqr(x4)*x3*y1+y3*x2*Sqr(x4)-y3*Sqr(x2)*x4-x4*Sqr(x1)*y2+x4*Sqr(x1)*y3+Sqr(x3)*y2*x4-x4*y1*Sqr(x3)+x4*y1*Sqr(x2)-Sqr(x1)*y4*x3-y4*x2*Sqr(x3)+Sqr(x1)*x3*y2-y1*x3*Sqr(x2)+y3*x1*Sqr(x2)-Sqr(x3)*x1*y2-x1*y4*Sqr(x2)+y4*x3*Sqr(x2)+Sqr(x1)*x2*y4-Sqr(x1)*x2*y3+x1*y4*Sqr(x3)+Sqr(x3)*x2*y1)/
            ((x4-x3)*(x4*x2-x4*x3+x2*x3-Sqr(x2))*(x3*x4*x2-x4*x1*x2+x4*Sqr(x1)-x4*x1*x3+Sqr(x1)*x2-x2*x1*x3-Sqr(x1)*x1+Sqr(x1)*x3));
        cc:=(x3*y2-y2*x4-y4*x3+x4*y3-x4*dd*Sqr(x3)*x3+dd*Sqr(x4)*x4*x3-x2*y3+x2*y4-x2*dd*Sqr(x4)*x4+x2*dd*Sqr(x3)*x3-dd*Sqr(x2)*x2*x3+dd*Sqr(x2)*x2*x4)/((x4-x3)*(x4*x2-x4*x3+x2*x3-Sqr(x2)));
        bb:=(y4-y3+cc*Sqr(x3)-cc*Sqr(x4)-dd*Sqr(x4)*x4+dd*Sqr(x3)*x3)/(x4-x3);
        aa:=y4-bb*x4-cc*Sqr(x4)-dd*Sqr(x4)*x4;
      end;
      // vytvorit body
      i:=0;
      x:=Xmin;
      k:=0;
      while (k<n) and (points[k].x<x) do Inc(k);
      while x<=Xmax do begin
        if i>=nl then begin
          nl:=nl+100;
          SetLength(line, nl);
        end;
        line[i].x:=x;
        if linearne
        then line[i].y:=LinInterpSafe(x, x1, y1, x4, y4)
        else line[i].y:=aa+bb*x+cc*Sqr(x)+dd*Sqr(x)*x;
        Inc(i);
        x:=x+lineStep;
        while throughPoints and (k<n) and (points[k].x<x) do begin
          x:=points[k].x;
          Inc(k);
        end;
        if throughPoints and (x>Xmax) and (i>0) and not IsEqual(line[i-1].x, Xmax) then x:=Xmax;
      end;
      SetLength(line, i);
    end
    else begin // kubicky splajn
      SetLength(h, n);
      h[0]:=0;
      for i:=1 to n-1 do h[i]:=Max(minDist, points[i].x-points[i-1].x);

      SetLength(a, n-2);
      SetLength(b, n-2);
      SetLength(c, n-2);
      SetLength(M, n-2);
      SetLength(r, n-2);
      for i:=0 to n-3 do begin
        a[i]:=h[i+1];
        b[i]:=2*(h[i+1]+h[i+2]);
        c[i]:=h[i+2];
        r[i]:=6*(1/h[i+1]*points[i].y-(1/h[i+1]+1/h[i+2])*points[i+1].y+1/h[i+2]*points[i+2].y);
      end;
      a[0]:=0;
      c[n-3]:=0;
      Tridiag(a, b, c, M, r);
      SetLength(M, n);
      for i:=n-3 downto 0 do M[i+1]:=M[i];
      M[0]:=0;
      M[n-1]:=0;
      // vytvorit body
      i:=0;
      x:=Xmin;
      kk:=0;
      while (kk<n) and (points[kk].x<x) do Inc(kk);
      k:=1;
      while (k<n-1) and (points[k].x<x) do Inc(k);
      while x<=Xmax do begin
        if i>=nl then begin
          nl:=nl+100;
          SetLength(line, nl);
        end;
        line[i].x:=x;
        xk:=points[k].x-x;
        xk1:=x-points[k-1].x;
        line[i].y:=(M[k-1]*xk*sqr(xk)+M[k]*xk1*sqr(xk1))/(6*h[k])+((points[k-1].y-M[k-1]*sqr(h[k])/6)*xk+(points[k].y-M[k]*sqr(h[k])/6)*xk1)/h[k];
        Inc(i);
        x:=x+lineStep;
        while throughPoints and (kk<n) and (points[kk].x<x) do begin
          x:=points[kk].x;
          Inc(kk);
          k:=Min(kk, n-1);
        end;
        while (k<n-1) and (points[k].x<x) do begin
          if throughPoints then x:=points[k].x;
          Inc(k);
        end;
        if throughPoints and (x>Xmax) and (i>0) and not IsEqual(line[i-1].x, Xmax) then x:=Xmax;
      end;
      SetLength(line, i);
    end;
  except
    SetLength(line, 0);
  end;
  SetLength(h, 0);
  SetLength(a, 0);
  SetLength(b, 0);
  SetLength(c, 0);
  SetLength(M, 0);
  SetLength(r, 0);
end;

procedure FFT(var data: TArrXYDouble; inverse: boolean=false);
var
  numBits, i, j, k, n, blockSize, blockEnd, num, m: word;
  delta_angle, delta_ar: double;
  alpha, beta: double;
  tr, ti, ar, ai: double;
  angleNumerator: double;
  tempXY: TXYDouble;
  function IsPowerOfTwo(x: longword): boolean;
  var
    i, y:  longword;
  begin
    y:=2;
    for i:=1 to 31 do begin
      if x=y then begin
        IsPowerOfTwo:=true;
        exit;
      end;
      y:=y shl 1;
    end;
    IsPowerOfTwo:=false;
  end;
  function NumberOfBitsNeeded(powerOfTwo: word): word;
  var
    i: word;
  begin
    NumberOfBitsNeeded:=0;
    for i:=0 to 16 do if (powerOfTwo and (1 shl i)) <> 0 then begin
      NumberOfBitsNeeded:=i;
      exit;
    end;
  end;
  function ReverseBits(index, numBits: word): word;
  var
    i, rev: word;
  begin
    rev:=0;
    for i:=0 to numBits-1 do begin
      rev:=(rev shl 1) or (index and 1);
      index:=index shr 1;
    end;
    ReverseBits := rev;
  end;
begin
  num:=Length(data);
  if not IsPowerOfTwo(num) then exit; // Number of samples is not a positive integer power of 2

  numBits := NumberOfBitsNeeded(num);

  j:=0;
  for i:=0 to num-1 do begin
    if j>i then begin
      tempXY:=data[i];
      data[i]:=data[j];
      data[j]:=tempXY;
    end;
    m:=num shr 1;
    while (m>=2) and (j>=m) do begin
      j:=j-m;
      m:=m shr 1;
    end;
    j:=j+m;
  end;

  angleNumerator:=6.283185307179586;
  if inverse then angleNumerator:=-angleNumerator;
  blockEnd:=1;
  blockSize:=2;
  while blockSize<=num do begin
      delta_angle:=angleNumerator/blockSize;
      alpha:=sin(0.5*delta_angle);
      alpha:=2.0*alpha*alpha;
      beta:=sin(delta_angle);
      i:=0;
      while i<num do begin
          ar:=1.0;    (* cos(0) *)
          ai:=0.0;    (* sin(0) *)
          j:=i;
          for n:=0 to blockEnd-1 do begin
              k:=j+blockEnd;
              tr:=ar*data[k].X-ai*data[k].Y;
              ti:=ar*data[k].Y+ai*data[k].X;
              data[k].X:=data[j].X-tr;
              data[k].Y:=data[j].Y-ti;
              data[j].X:=data[j].X+tr;
              data[j].Y:=data[j].Y+ti;
              delta_ar:=alpha*ar+beta*ai;
              ai:=ai-(alpha*ai-beta*ar);
              ar:=ar-delta_ar;
              Inc(j);
          end;
          i:=i+blockSize;
      end;
      blockEnd:=blockSize;
      blockSize:=blockSize shl 1;
  end;
end;

function FitImageWithinRect(var ARect: TRect; ImageWidth, ImageHeight, RectWidth, RectHeight: integer): boolean;
var
  width, height: integer;
begin
  result:=FitImageWithinRect(ARect.Left, ARect.Top, width, height, ImageWidth, ImageHeight, RectWidth, RectHeight);
  ARect.Right:=ARect.Left+width;
  ARect.Bottom:=ARect.Top+height;
end;

function FitImageWithinRect(var left, top, width, height: integer; ImageWidth, ImageHeight, RectWidth, RectHeight: integer): boolean;
var
  kx, ky, k: double;
begin
  result:=true;
  try
    kx:=RectWidth/ImageWidth;
    ky:=RectHeight/ImageHeight;
    if kx<ky
    then k:=kx
    else k:=ky;
    width:=Round(k*ImageWidth);
    height:=Round(k*ImageHeight);
    left:=(RectWidth-width) div 2;
    top:=(RectHeight-height) div 2;
  except
    left:=0;
    top:=0;
    width:=0;
    height:=0;
    result:=false;
  end;
end;

function FitImageWithinRect(var left, top, width, height: double; ImageWidth, ImageHeight, RectWidth, RectHeight: double): boolean;
var
  kx, ky, k: double;
begin
  result:=true;
  try
    kx:=RectWidth/ImageWidth;
    ky:=RectHeight/ImageHeight;
    if kx<ky
    then k:=kx
    else k:=ky;
    width:=k*ImageWidth;
    height:=k*ImageHeight;
    left:=(RectWidth-ImageWidth) / 2;
    top:=(RectHeight-ImageHeight) / 2;
  except
    left:=0;
    top:=0;
    width:=0;
    height:=0;
    result:=false;
  end;
end;

function FixRectWtoH(r: TRect; WtoH: double): TRect; // shrinks the rect so that the returned rect has width/height=WtoH
var
  k: double;
  newDim, dim2: integer;
begin
  result:=r;
  if (r.Right>r.Left) and (r.Bottom>r.Top) then try
    k:=(r.Right-r.Left)/(r.Bottom-r.Top);
    if k>WtoH
    then begin // ramecek je sirsi
      newDim:=Trunc(WtoH*(r.Bottom-r.Top));
      dim2:=((r.Right-r.Left)-newDim) div 2;
      r.Left:=r.Left+dim2;
      r.Right:=r.Right-dim2;
    end
    else begin // ramecek je vyssi
      newDim:=Trunc((r.Right-r.Left)/WtoH);
      dim2:=((r.Bottom-r.Top)-newDim) div 2;
      r.Top:=r.Top+dim2;
      r.Bottom:=r.Bottom-dim2;
    end;
  except
  end;
  result:=r;
end;

function HalfCycleRad(alpha: double): double; // returns alpha within range -Pi..Pi
begin
  result:=Frac(alpha/(2*Pi))*2*Pi;
  if result<-Pi then result:=result+2*Pi;
  if result>Pi then result:=result-2*Pi;
end;

function HalfCycleDeg(alpha: double): double; // returns alpha within range -180..180
begin
  result:=HalfCycleRad(alpha*Pi/180)*180/Pi;
end;

function OneCycleRad(alpha: double): double; // returns alpha within range -Pi..Pi
begin
  result:=Frac(alpha/(2*Pi))*2*Pi;
  if result<0 then result:=result+2*Pi;
end;

function OneCycleDeg(alpha: double): double; // returns alpha within range -180..180
begin
  result:=Frac(alpha/360)*360;
  if result<0 then result:=result+360;
end;

function LineIntersect(var x, y: double; ax1, ay1, ax2, ay2: double; bx1, by1, bx2, by2: double): boolean;
var
  del: double;
  ta: double;
begin
  result:=true;
  try
    del:=ax2*by2-ax2*by1-ax1*by2+ax1*by1+ay2*bx1-ay2*bx2-ay1*bx1+ay1*bx2;
    if PohZero(del) then raise Exception.Create('Parallel lines');
    ta:=(-ax1*by2+ax1*by1+ay1*bx2-ay1*bx1+bx1*by2-bx2*by1)/del;
    x:=ax1+ta*(ax2-ax1);
    y:=ay1+ta*(ay2-ay1);
  except
    result:=false;
    x:=0;
    y:=0;
  end;
end;

function Circle2Intersect(var x1, y1, x2, y2: double; cx1, cy1, r1, cx2, cy2, r2: double): boolean; // circle: (x-cx1)^2+(y-cy)^2=r
var
  offx, offy, cdist, a: double;
  del, det: double;
  AA: double;
begin
  result:=true;
  cdist:=PohDist(cx1, cy1, cx2, cy2);
  try
    if IsEqual(r1, r2) and IsEqual(cx1, cx2) and IsEqual(cy1, cy2) then begin // kruznice splyvaji
      x1:=-r1+cx1;
      y1:=cy1;
      x2:=r1+cx1;
      y2:=cy1;
      exit;
    end;
    if ((cdist>r1+r2) and not IsEqual(cdist, r1+r2)) or
       ((cdist+r1<r2) and not IsEqual(cdist+r1, r2)) or
       ((cdist+r2<r1) and not IsEqual(cdist+r2, r1))
    then raise Exception.Create('No intersection');
    // posunuti kruznice1 do pocatku
    offx:=cx1;
    offy:=cy1;
    // cx1:=0; // nebude pouzito, protoze je to nula
    // cy1:=0; // nebude pouzito, protoze je to nula
    cx2:=cx2-offx;
    cy2:=cy2-offy;
    if IsEqual(cdist, r1+r2) or IsEqual(cdist+r2, r1)
    then begin // kruznice se dotykaji
      a:=GetAngleRad(0, 0, cx2, cy2);
      RotatePointRad(x1, y1, 0, -r1, -a);
      x2:=x1;
      y2:=y1;
    end
    else if IsEqual(cdist+r1, r2)
    then begin
      a:=GetAngleRad(0, 0, cx2, cy2);
      RotatePointRad(x1, y1, 0, -r1, -(a+Pi));
      x2:=x1;
      y2:=y1;
    end
    else begin // normalni dva pruseciky
      del:=2*(cx2*cx2+cy2*cy2);
      det:=-sqr(cx2)*(
             sqr(sqr(r2))-2*sqr(cx2*r2)+2*sqr(cx2*cy2)+sqr(sqr(cx2))
             -2*sqr(r1*cy2)-2*sqr(r1*r2)-2*sqr(cy2*r2)+sqr(sqr(cy2))-2*sqr(r1*cx2)+sqr(sqr(r1))
            );
      AA:=sqr(r1)*cy2+cy2*sqr(cx2)+sqr(cy2)*cy2-cy2*sqr(r2);
      if det<0 then raise Exception.Create('No intersection');
      // y1, y2
      y1:=(AA+sqrt(det))/del;
      y2:=(AA-sqrt(det))/del;
      // x1
      det:=sqr(r1)-sqr(y1);
      if det<0 then raise Exception.Create('No intersection');
      x1:=sqrt(det);
      if not IsEqual(sqr(r2), sqr(x1-cx2)+sqr(y1-cy2), 8) then x1:=-x1;
      // x2
      det:=sqr(r1)-sqr(y2);
      if det<0 then raise Exception.Create('No intersection');
      x2:=sqrt(det);
      if not IsEqual(sqr(r2), sqr(x2-cx2)+sqr(y2-cy2), 8) then x2:=-x2;
      if IsEqual(x1, x2, 8) and IsEqual(y1, y2, 8) then x2:=-x2;
    end;

    // posunoti vysledku zpet
    x1:=x1+offx;
    y1:=y1+offy;
    x2:=x2+offx;
    y2:=y2+offy;
  except
    x1:=0;
    y1:=0;
    x2:=0;
    y2:=0;
    result:=false;
  end;
end;

procedure PointOnLine(var p: TXYDouble; p1, p2: TXYDouble; GivenPoint: char);
var
  v: double;
begin
  if UpperCase(GivenPoint)='X'
  then begin
    v:=p2.x-p1.x;
    if PohZero(v)
    then p.y:=p1.y
    else p.y:=(p.x-p1.x)/v*(p2.y-p1.y)+p1.y;
  end
  else if UpperCase(GivenPoint)='Y'
  then begin
    v:=p2.y-p1.y;
    if PohZero(v)
    then p.x:=p1.x
    else p.x:=(p.y-p1.y)/v*(p2.x-p1.x)+p1.x;
  end;
end;

procedure PointOnLine(var p: TXYZDouble; p1, p2: TXYZDouble; GivenPoint: string);
var
  v, v1, v2: double;
begin
  if UpperCase(GivenPoint)='X'
  then begin
    v:=p2.x-p1.x;
    if PohZero(v)
    then begin
      p.y:=p1.y;
      p.z:=p1.z;
    end
    else begin
      p.y:=(p.x-p1.x)/v*(p2.y-p1.y)+p1.y;
      p.z:=(p.x-p1.x)/v*(p2.z-p1.z)+p1.z;
    end;
  end
  else if UpperCase(GivenPoint)='Y'
  then begin
    v:=p2.y-p1.y;
    if PohZero(v)
    then begin
      p.x:=p1.x;
      p.z:=p1.z;
    end
    else begin
      p.x:=(p.y-p1.y)/v*(p2.x-p1.x)+p1.x;
      p.z:=(p.y-p1.y)/v*(p2.z-p1.z)+p1.z;
    end;
  end
  else if UpperCase(GivenPoint)='Z'
  then begin
    v:=p2.z-p1.z;
    if PohZero(v)
    then begin
      p.x:=p1.x;
      p.y:=p1.y;
    end
    else begin
      p.x:=(p.z-p1.z)/v*(p2.x-p1.x)+p1.x;
      p.y:=(p.z-p1.z)/v*(p2.y-p1.y)+p1.y;
    end;
  end;

  if UpperCase(GivenPoint)='XY'
  then begin
    v1:=p2.x-p1.x;
    v2:=p2.y-p1.y;
    if PohZero(v1) and PohZero(v2)
    then p.z:=p1.z
    else begin
      if Abs(v1)>Abs(v2)
      then p.z:=(p.x-p1.x)/v1*(p2.z-p1.z)+p1.z
      else p.z:=(p.y-p1.y)/v2*(p2.z-p1.z)+p1.z;
    end;
  end
  else if UpperCase(GivenPoint)='XZ'
  then begin
    v1:=p2.x-p1.x;
    v2:=p2.z-p1.z;
    if PohZero(v1) and PohZero(v2)
    then p.y:=p1.y
    else begin
      if Abs(v1)>Abs(v2)
      then p.y:=(p.x-p1.x)/v1*(p2.y-p1.y)+p1.y
      else p.y:=(p.z-p1.z)/v2*(p2.y-p1.y)+p1.y;
    end;
  end
  else if UpperCase(GivenPoint)='YZ'
  then begin
    v1:=p2.y-p1.y;
    v2:=p2.z-p1.z;
    if PohZero(v1) and PohZero(v2)
    then p.x:=p1.x
    else begin
      if Abs(v1)>Abs(v2)
      then p.x:=(p.y-p1.y)/v1*(p2.x-p1.x)+p1.x
      else p.x:=(p.z-p1.z)/v2*(p2.x-p1.x)+p1.x;
    end;
  end;
end;

function PointInLineSegment(p1x, p1y, p2x, p2y, testPointx, testPointy: double; precision: double): boolean; // tests ranges only, not position on line
begin
  if p1x>p2x then SwapDouble(p1x, p2x);
  if p1y>p2y then SwapDouble(p1y, p2y);
  result:=((p1x-precision<=testPointx) and (testPointx<=p2x+precision) and
           (p1y-precision<=testPointy) and (testPointy<=p2y+precision));
end;

function PointInLineSegment(p1, p2, testPoint: TXYZDouble; precision: double): boolean; // tests ranges only, not position on line
begin
  if p1.x>p2.x then SwapDouble(p1.x, p2.x);
  if p1.y>p2.y then SwapDouble(p1.y, p2.y);
  if p1.z>p2.z then SwapDouble(p1.z, p2.z);
  result:=((p1.x-precision<=testPoint.x) and (testPoint.x<=p2.x+precision) and
           (p1.y-precision<=testPoint.y) and (testPoint.y<=p2.y+precision) and
           (p1.z-precision<=testPoint.z) and (testPoint.z<=p2.z+precision));
end;

function PointInRect(x, y, x1, y1, x2, y2: double): boolean;
var
  d: double;
begin
  if x1>x2 then begin
    d:=x1;
    x1:=x2;
    x2:=d;
  end;
  if y1>y2 then begin
    d:=y1;
    y1:=y2;
    y2:=d;
  end;
  result:=(x1<=x) and (x<=x2) and (y1<=y) and (y<=y2);
end;

function PohDist(x1, y1: double): double;
begin
  result:=sqrt(max(0, sqr(x1)+sqr(y1)));
end;

function PohDist(x1, y1: extended): extended;
begin
  result:=sqrt(max(0, sqr(x1)+sqr(y1)));
end;

function PohDist(a, b: TXYSingle): single;
begin
  result:=sqrt(max(0, sqr(b.x-a.x)+sqr(b.y-a.y)));
end;

function PohDist(a, b: TXYDouble): double;
begin
  result:=sqrt(max(0, sqr(b.x-a.x)+sqr(b.y-a.y)));
end;

function PohDist(a, b: TXYZDouble): double;
begin
  result:=sqrt(max(0, sqr(b.x-a.x)+sqr(b.y-a.y)+sqr(b.z-a.z)));
end;

function PohDist(x1, y1, x2, y2: double): double;
begin
  result:=sqrt(max(0, sqr(x2-x1)+sqr(y2-y1)));
end;

function PohDist(x1, y1, x2, y2: extended): extended;
begin
  result:=sqrt(max(0, sqr(x2-x1)+sqr(y2-y1)));
end;

function PohDist(x1, y1, z1, x2, y2, z2: double): double;
begin
  result:=sqrt(max(0, sqr(x2-x1)+sqr(y2-y1)+sqr(z2-z1)));
end;

function IsEqualSingle(a, b: single; dig: integer = 6): boolean; // dig is number of left most digits to compare
var
  nda, ndb: integer;
begin
  if (Abs(a)<1e-100) and (Abs(b)<1e-100) then begin
    result:=true;
    exit;
  end;
  try
    if (Abs(a)<1e-100)
    then nda:=-101
    else nda:=LessInt(Log10(Abs(a))+0.00001);
    if (Abs(b)<1e-100)
    then ndb:=-101
    else ndb:=LessInt(Log10(Abs(b))+0.00001);
    if Abs(nda-ndb)>1 then result:=false
    else result:=Abs(a-b)<Power(10, Max(nda, ndb)-dig);
  except
    result:=false;
  end;
end;


function IsEqual(a, b: double; dig: integer): boolean;
var
  nda, ndb: integer;
begin
  if (Abs(a)<1e-100) and (Abs(b)<1e-100) then begin
    result:=true;
    exit;
  end;
  try
    if (Abs(a)<1e-100)
    then nda:=-101
    else nda:=LessInt(Log10(Abs(a))+0.0000001);
    if (Abs(b)<1e-100)
    then ndb:=-101
    else ndb:=LessInt(Log10(Abs(b))+0.0000001);
    if Abs(nda-ndb)>1 then result:=false
    else result:=Abs(a-b)<Power(10, Max(nda, ndb)-dig);
  except
    result:=false;
  end;
end;

function IsEqualExt(a, b: extended; dig: integer): boolean;
var
  nda, ndb: integer;
begin
  if (Abs(a)<1e-1930) and (Abs(b)<1e-1930) then begin
    result:=true;
    exit;
  end;
  try
    if (Abs(a)<1e-1930)
    then nda:=-1931
    else nda:=LessInt(Log10(Abs(a))+0.0000001);
    if (Abs(b)<1e-1930)
    then ndb:=-1931
    else ndb:=LessInt(Log10(Abs(b))+0.0000001);
    if Abs(nda-ndb)>1 then result:=false
    else result:=Abs(a-b)<Power(10, Max(nda, ndb)-dig);
  except
    result:=false;
  end;
end;

function IsEqual(a, b: TXYSingle; dig: integer = 6): boolean;
begin
  result:=IsEqualSingle(a.x, b.x, dig) and IsEqualSingle(a.y, b.y, dig);
end;

function IsEqual(a, b: TXYDouble; dig: integer): boolean;
begin
  result:=IsEqual(a.x, b.x, dig) and IsEqual(a.y, b.y, dig);
end;

function IsEqual(a, b: TXYZDouble; dig: integer): boolean;
begin
  result:=IsEqual(a.x, b.x, dig) and IsEqual(a.y, b.y, dig) and IsEqual(a.z, b.z, dig);
end;

function PohSign(a: double): TValueSign; // returns -1, 0, 1
begin
  if PohZero(a)
  then result:=0
  else if a>=0 then result:=1
  else result:=-1;
end;

function PohSign(a: extended): TValueSign; // returns -1, 0, 1
begin
  if PohZero(a)
  then result:=0
  else if a>=0 then result:=1
  else result:=-1;
end;

function PohZero(a: single): boolean;
begin
  result:=Abs(a)<1e-30;
end;

function PohZero(a: double): boolean;
begin
  result:=Abs(a)<1e-100;
end;

function PohZero(a: extended): boolean;
begin
  result:=Abs(a)<1e-1930;
end;

function LessInt(x: double): integer;
begin
  if x>=0
  then result:=Trunc(x)
  else result:=Trunc(x)-1;
end;

procedure GetEllipseRange(var wx, wy, hx, hy: double; a, b, alphaRad: double);
var
  q, alpha: double;
  det: double;
begin
  if a<0 then a:=0;
  if b<0 then b:=0;
  try
    alphaRad:=Frac(alphaRad/Pi)*Pi;
    if alphaRad<-Pi/2 then alphaRad:=-Pi-alphaRad;
    if alphaRad>+Pi/2 then alphaRad:=Pi-alphaRad;
    if Abs(alphaRad)<0.01*Pi/180
    then begin// nenatoceno
      wx:=a;
      wy:=0;
    end
    else begin
      q:=Cot(alphaRad);
      det:=sqrt(b*b+q*q*a*a);
      wx:=q*a*a/det;
      wy:=-b*b/det;
      RotatePointRad(wx, wy, wx, wy, -alphaRad);
    end;
    alpha:=alphaRad;
    alphaRad:=Frac((alphaRad+Pi/2)/Pi)*Pi;
    if Abs(alphaRad)<0.01*Pi/180
    then begin// nenatoceno
      hx:=0;
      hy:=-a;
    end
    else begin
      q:=Cot(alphaRad);
      det:=sqrt(b*b+q*q*a*a);
      hx:=-q*a*a/det;
      hy:=b*b/det;
      RotatePointRad(hx, hy, hx, hy, -alpha);
    end;
    if wx<0
    then wx:=-wx
    else wy:=-wy;
    if hy>0
    then hy:=-hy
    else hx:=-hx;
  except
    wx:=a;
    wy:=0;
    hx:=0;
    hy:=-b;
  end;
end;

function GetEllipseVolumePart(a, b, x1, x2: double): double;
var
  det1, det2: double;
begin
  x1:=Min(a, Max(-a, x1));
  x2:=Min(a, Max(-a, x2));
  if x1>x2 then SwapDouble(x1, x2);
  try
    det1:=sqrt(Max(0, a*a-x1*x1));
    if abs(det1/x1)<1e-30
    then det1:=(b*x1*det1/a+b*a*Pi/2)/2
    else det1:=(b*x1*det1/a+b*a*arctan(x1/det1))/2;
    det2:=sqrt(Max(0, a*a-x2*x2));
    if abs(det2/x2)<1e-30
    then det2:=(b*x2*det2/a+b*a*Pi/2)/2
    else det2:=(b*x2*det2/a+b*a*arctan(x2/det2))/2;
    result:=2*(det2-det1); // nad x i pod x osou
  except
    result:=0;
  end;
end;

procedure RotatePointDeg(var DestX, DestY: double; OrigX, OrigY, Alpha: double);
begin
  RotatePointRad(DestX, DestY, OrigX, OrigY, Alpha*Pi/180);
end;

procedure RotatePointRad(var DestX, DestY: double; OrigX, OrigY, Alpha: double);
var
  ca, sa: Extended;
begin
  SinCos(Alpha, sa, ca);
  DestX:=OrigX*ca+OrigY*sa;
  DestY:=OrigY*ca-OrigX*sa;
end;

procedure RotatePointDeg(var DestX, DestY: double; OrigX, OrigY, Alpha, CenterX, CenterY: double);
begin
  RotatePointRad(DestX, DestY, OrigX, OrigY, Alpha*Pi/180, CenterX, CenterY);
end;

procedure RotatePointRad(var DestX, DestY: double; OrigX, OrigY, Alpha, CenterX, CenterY: double);
var
  ca, sa: Extended;
begin
  SinCos(Alpha, sa, ca);
  DestX:=CenterX+(OrigX-CenterX)*ca+(OrigY-CenterY)*sa;
  DestY:=CenterY-(OrigX-CenterX)*sa+(OrigY-CenterY)*ca;
end;

function CircleLineIntersect(var x1, y1, x2, y2: double; r: double; y0, q: double): boolean; // circle: x^2+y^2=r; line:=y=y0+q*x;
var
  det, A, B: double;
begin
  try
    result:=true;
    if (r<0) then r:=0;
    det:=r*r-y0*y0+q*q*r*r;
    if det<0 then begin
      x1:=0;
      y1:=0;
      x2:=0;
      y2:=0;
      result:=false;
      exit;
    end;
    A:=1+q*q;
    B:=-y0*q;
    det:=sqrt(det);
    x1:=(B-det)/A;
    x2:=(B+det)/A;
    y1:=y0+q*x1;
    y2:=y0+q*x2;
  except
    x1:=0;
    y1:=0;
    x2:=0;
    y2:=0;
    result:=false;
  end;
end;

function CircleLineTanIntersect(var x1, y1, x2, y2: double; r: double; x0, y0, alphaRad: double): boolean; // circle: x^2+y^2=r; line: primka prochazejici bodem [x0, y0] pod uhlem alfa od osy X;
var
  q: double;
  chyba: boolean;
begin
  if r<0 then r:=0;
  try
    chyba:=false;
    q:=tan(alphaRad);
  except
    chyba:=true;
    q:=0;
  end;
  if chyba or (Abs(q)>1e10)
  then begin // primka je rovnobezna s Y osou
    if UPohMath.IsEqual(-r, x0)
    then begin // tecna v -r
      x1:=-r;
      y1:=0;
      x2:=-r;
      y2:=0;
      result:=true;
    end
    else if UPohMath.IsEqual(r, x0)
    then begin // tecna v r
      x1:=r;
      y1:=0;
      x2:=r;
      y2:=0;
      result:=true;
    end
    else if (x0<-r) or (r<x0)
    then begin // neprotina
      x1:=0;
      y1:=0;
      x2:=0;
      y2:=0;
      result:=false;
    end
    else begin // protina
      x1:=x0;
      y1:=-sqrt(r*r-x0*x0);
      x2:=x0;
      y2:=-y1;
      result:=true;
    end;
  end
  else try
    result:=CircleLineIntersect(x1, y1, x2, y2, r, y0-q*x0, q);
  except
    x1:=0;
    y1:=0;
    x2:=0;
    y2:=0;
    result:=false;
  end;
end;

function CircleLineCotIntersect(var x1, y1, x2, y2: double; r: double; x0, y0, alphaRad: double): boolean; // circle: x^2+y^2=r; line: primka prochazejici bodem [x0, y0] pod uhlem alfa od osy Y;
begin
  result:=CircleLineTanIntersect(x1, y1, x2, y2, r, x0, y0, Pi/2-alphaRad);
end;

procedure SwapDouble(var a, b: double);
var
  c: double;
begin
  c:=a;
  a:=b;
  b:=c;
end;

procedure SwapInt(var a, b: integer);
var
  c: integer;
begin
  c:=a;
  a:=b;
  b:=c;
end;

function LinInterp(x, x0, y0, x1, y1: double): double; // returns (y1-y0)/(x1-x0)*(x-x0)+y0; x1<>x0
begin
  result:=(y1-y0)/(x1-x0)*(x-x0)+y0;
end;

function LinInterpSafe(x, x0, y0, x1, y1: double): double; // returns (y1-y0)/(x1-x0)*(x-x0)+y0; if x1=x0 then average is returned (y0+y1)/2
begin
  if IsEqual(x0,x1) then result:=(y1+y0)/2
  else result:=(y1-y0)/(x1-x0)*(x-x0)+y0;
end;

function QuadrInterp(x, y, x1, y1, x2, y2, v11, v21, v12, v22: double): double; // vXY
begin
  result:=LinInterp(y, y1, LinInterp(x, x1, v11, x2, v21), y2, LinInterp(x, x1, v12, x2, v22));
end;

function QuadrInterpSafe(x, y, x1, y1, x2, y2, v11, v21, v12, v22: double): double; // vXY
begin
  result:=LinInterpSafe(y, y1, LinInterpSafe(x, x1, v11, x2, v21), y2, LinInterpSafe(x, x1, v12, x2, v22));
end;

function GetAngleRad(x1, y1, x2, y2: double): double; // 0rad = 12 clock, 2Pi rad = 24 clock; x and y like display coordinates
var
  l, a: double;
begin
  try
    l:=sqrt(sqr(x2-x1)+sqr(y2-y1));
    a:=arcsin(Max(-1,Min(1,(y2-y1)/l)))+Pi/2;
    if (x2-x1)<0 then a:=2*Pi-a;
    result:=a;
  except
    result:=0;
  end;
end;

function GetAngleDeg(x1, y1, x2, y2: double): double; // 0o = 12 clock, 360o = 24 clock; x and y like display coordinates
begin
  result:=GetAngleRad(x1, y1, x2, y2)*180/Pi;
end;

function Tridiag(var a, b, c, x, r: TDoubleArray): boolean; // Solves for a vector x[0..n-1] the tridiagonal linear set given by a[0..n-1], b[0..n-1], c[0..n-1], and "right side"=r[0..n-1] are input vectors and are not modied. X vector must be allocated.
var
  i, n, n1: integer;
  bet: double;
  gam: TDoubleArray;
begin
  result:=false;
  n:=Length(a);
  n1:=n-1;
  if (Length(a)<>Length(b)) or (Length(b)<>Length(c)) or (Length(c)<>Length(r)) or PohZero(b[0]) or (n<3) then exit;

  SetLength(gam, n);
  bet:=b[0];
  x[0]:=r[0]/bet;

  for i:=1 to n1 do begin // Decomposition and forward substitution.
    gam[i]:=c[i-1]/bet;
    bet:=b[i]-a[i]*gam[i];
    if Abs(bet)<1e-300 then begin // Algorithm fails
      SetLength(gam, 0);
      exit;
    end;
    x[i]:=(r[i]-a[i]*x[i-1])/bet;
  end;
  for i:=n-2 downto 0 do x[i]:=x[i]-gam[i+1]*x[i+1]; // Backsubstitution.

  SetLength(gam, 0);
  result:=true;
end;

procedure OptimParabInit(var values: TOptimParab; a, b, c, fc: double; AccuracyFrac: double = 1e-6; ZeroTreshold: double = 1e-6); // c must be within a and b, fc is function value in c; max odchylka je AccuracyFrac*Abs(x)+ZeroTreshold, pokud je mensi je optimalizace u konce
var
  temp: double;
begin
  values.a:=a;
  values.b:=b;
  values.AccuracyFrac:=AccuracyFrac;
  values.ZeroTreshold:=ZeroTreshold;
  with values do begin
    if a>b then begin
      temp:=a;
      a:=b;
      b:=temp;
    end;
    e:=0;
    d:=0;
    x:=c;
    v:=c;
    w:=c;
    fx:=fc;
    fv:=fc;
    fw:=fc;
  end;
end;

function OptimParab1step(var values: TOptimParab): boolean; // true if done
const
  CGOLD = 0.381966;
var
  iter: integer;
  etemp, p, q, r, xm, tol1, tol2: double;
  function SIGN(a, b: double): double;
  begin
    if b>=0
    then result:=abs(a)
    else result:=-abs(a);
  end;
begin
  result:=false;
  with values do begin
    xm:=0.5*(a+b);
    tol1:=AccuracyFrac*abs(x)+ZeroTreshold;
    tol2:=2*tol1;
    if abs(x-xm)<=(tol2-0.5*(b-a)) then begin
      result:=true;
      exit;
    end;
    if abs(e)>tol1 then begin
      r:=(x-w)*(fx-fv);
      q:=(x-v)*(fx-fw);
      p:=(x-v)*q-(x-w)*r;
      q:=2.0*(q-r);
      if q>0 then p:=-p;
      q:=abs(q);
      etemp:=e;
      e:=d;
      if (abs(p)>=abs(0.5*q*etemp)) or (p<=q*(a-x)) or (p>=q*(b-x))
      then begin
        if x>=xm
        then e:=(a-x)
        else e:=(b-x);
        d:=CGOLD*e;
      end
      // The above conditions determine the acceptability of the parabolic fit. Here we take the golden section step into the larger of the two segments.
      else begin
        d:=p/q; // Take the parabolic step.
        u:=x+d;
        if (u-a<tol2) or  (b-u<tol2) then d:=SIGN(tol1, xm-x);
      end;
    end
    else begin
      if x>=xm
      then e:=a-x
      else e:=b-x;
      d:=CGOLD*e;
    end;

    if abs(d)>=tol1
    then u:=x+d
    else u:=x+SIGN(tol1,d);
    // Now is the one function evaluation per iteration.
    // compute fu for u
  end;
end;

function OptimParab2step(var values: TOptimParab): boolean; // true if done; after this function compute fu for u
var
  xm, tol1, tol2: double;
begin
  result:=false;
  with values do begin
    if fu<fx
    then begin
      if u>=x
      then a:=x
      else b:=x;
      // SHFT(v,w,x,u);
      v:=w;
      w:=x;
      x:=u;
      // SHFT(fv,fw,fx,fu);
      fv:=fw;
      fw:=fx;
      fx:=fu;
    end
    else begin
      if u<x
      then a:=u
      else b:=u;
      if (fu<=fw) or (abs(w-x)<ZeroTreshold)
      then begin
        v:=w;
        w:=u;
        fv:=fw;
        fw:=fu;
      end
      else begin
        if (fu<=fv) or (abs(v-x)<ZeroTreshold) or (abs(v-w)<ZeroTreshold)
        then begin
          v:=u;
          fv:=fu;
        end;
      end;
    end;
    // test end of optimization
    xm:=0.5*(a+b);
    tol1:=AccuracyFrac*abs(x)+ZeroTreshold;
    tol2:=2*tol1;
    if abs(x-xm)<=(tol2-0.5*(b-a)) then begin
      result:=true;
      exit;
    end;
  end;
end;

function T_K(mV, T0: double): double;
begin
  if mV<0 then result:=2.5173462E+01*mV-1.1662878E+00*power(mV,2)-1.0833638E+00*power(mV,3)-8.9773540E-01*power(mV,4)-3.7342377E-01*power(mV,5)-8.6632643E-02*power(mV,6)-1.0450598E-02*power(mV,7)-5.1920577E-04*power(mV,8)+T0
  else if mV<20.819 then result:=2.5083550E+01*mV+7.8601060E-02*power(mV,2)-2.5031310E-01*power(mV,3)+8.3152700E-02*power(mV,4)-1.2280340E-02*power(mV,5)+9.8040360E-04*power(mV,6)-4.4130300E-05*power(mV,7)+1.0577340E-06*power(mV,8)-1.0527550E-08*power(mV,9)+T0
  else result:=-1.3180580E+02+4.8302220E+01*mV-1.6460310E+00*power(mV,2)+5.4647310E-02*power(mV,3)-9.6507150E-04*power(mV,4)+8.8021930E-06*power(mV,5)-3.1108100E-08*power(mV,6)+T0;
end;

procedure IdentityMatrix(var a: TDoubleMatrixAlgn; n: integer); // allocates 'a' if necessary
var
  c, r: integer;
begin
  SetLength(a, n);
  for c:=0 to n-1 do begin
    SetLength(a[c], n);
    for r:=0 to n-1 do a[c,r]:=0;
    a[c,c]:=1;
  end;
end;

procedure GaussFiltr(var dataIn, dataOut: TSingleArray; sigma: single);// filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
var
  i, k: integer;
  f: TSingleArray;
  n, n1, nf: integer;
  s2: single;
  sum, sumk: single;
begin
  nf:=Round(sigma*3+2);
  s2:=2*sqr(sigma);
  SetLength(f, nf+1);
  for k:=0 to nf do f[k]:=exp(-k*k/s2);
  n:=Length(dataIn);
  n1:=n-1;
  SetLength(dataOut, n);
  sumk:=0;
  for k:=-nf to -1 do sumk:=sumk+f[-k];
  sumk:=sumk+1;
  for k:=1 to nf do sumk:=sumk+f[k];
  if 2*nf<n1
  then begin
    for i:=0 to nf-1 do try
      sum:=0;
      for k:=-nf to -1 do sum:=sum+dataIn[Max(0, i+k)]*f[-k];
      sum:=sum+dataIn[i];
      for k:=1 to nf do sum:=sum+dataIn[i+k]*f[k];
      dataOut[i]:=sum/sumk;
    except
    end;
    for i:=nf to n1-nf do try
      sum:=0;
      for k:=-nf to -1 do sum:=sum+dataIn[i+k]*f[-k];
      sum:=sum+dataIn[i];
      for k:=1 to nf do sum:=sum+dataIn[i+k]*f[k];
      dataOut[i]:=sum/sumk;
    except
    end;
    for i:=n1-nf+1 to n1 do try
      sum:=0;
      for k:=-nf to -1 do sum:=sum+dataIn[i+k]*f[-k];
      sum:=sum+dataIn[i];
      for k:=1 to nf do sum:=sum+dataIn[Min(i+k, n1)]*f[k];
      dataOut[i]:=sum/sumk;
    except
    end;
  end
  else for i:=0 to n1 do try
    sum:=0;
    for k:=-nf to -1 do sum:=sum+dataIn[Max(0, i+k)]*f[-k];
    sum:=sum+dataIn[i];
    for k:=1 to nf do sum:=sum+dataIn[Min(i+k, n1)]*f[k];
    dataOut[i]:=sum/sumk;
  except
  end;
  SetLength(f, 0);
end;

procedure GaussFiltr(var data: TSingleArray; sigma: single); // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))
var
  dataNew: TSingleArray;
begin
  GaussFiltr(data, dataNew, sigma);
  data:=Copy(dataNew);
  SetLength(dataNew, 0);
end;

procedure GaussFiltr(var dataIn, dataOut: TDoubleArray; sigma: double);// filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
var
  i, k: integer;
  f: TDoubleArray;
  n, n1, nf: integer;
  s2: double;
  sum, sumk: double;
begin
  nf:=Round(sigma*3+2);
  s2:=2*sqr(sigma);
  SetLength(f, nf+1);
  for k:=0 to nf do f[k]:=exp(-k*k/s2);
  n:=Length(dataIn);
  n1:=n-1;
  SetLength(dataOut, n);
  sumk:=0;
  for k:=-nf to -1 do sumk:=sumk+f[-k];
  sumk:=sumk+1;
  for k:=1 to nf do sumk:=sumk+f[k];
  if 2*nf<n1
  then begin
    for i:=0 to nf-1 do try
      sum:=0;
      for k:=-nf to -1 do sum:=sum+dataIn[Max(0, i+k)]*f[-k];
      sum:=sum+dataIn[i];
      for k:=1 to nf do sum:=sum+dataIn[i+k]*f[k];
      dataOut[i]:=sum/sumk;
    except
    end;
    for i:=nf to n1-nf do try
      sum:=0;
      for k:=-nf to -1 do sum:=sum+dataIn[i+k]*f[-k];
      sum:=sum+dataIn[i];
      for k:=1 to nf do sum:=sum+dataIn[i+k]*f[k];
      dataOut[i]:=sum/sumk;
    except
    end;
    for i:=n1-nf+1 to n1 do try
      sum:=0;
      for k:=-nf to -1 do sum:=sum+dataIn[i+k]*f[-k];
      sum:=sum+dataIn[i];
      for k:=1 to nf do sum:=sum+dataIn[Min(i+k, n1)]*f[k];
      dataOut[i]:=sum/sumk;
    except
    end;
  end
  else for i:=0 to n1 do try
    sum:=0;
    for k:=-nf to -1 do sum:=sum+dataIn[Max(0, i+k)]*f[-k];
    sum:=sum+dataIn[i];
    for k:=1 to nf do sum:=sum+dataIn[Min(i+k, n1)]*f[k];
    dataOut[i]:=sum/sumk;
  except
  end;
  SetLength(f, 0);
end;

procedure GaussFiltr(var data: TDoubleArray; sigma: double); // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))
var
  dataNew: TDoubleArray;
begin
  GaussFiltr(data, dataNew, sigma);
  data:=Copy(dataNew);
  SetLength(dataNew, 0);
end;

procedure GaussFiltrAdapt(var dataIn, dataOut, sigma: TSingleArray); // filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
var
  i, k: integer;
  f: TSingleArray;
  n, n1, nf: integer;
  s2: single;
  sum, sumk: single;
begin
  n:=Length(dataIn);
  n1:=n-1;
  SetLength(dataOut, n);
  for i:=0 to n1 do try
    nf:=Round(sigma[i]*3+2);
    s2:=2*sqr(sigma[i]);
    SetLength(f, nf+1);
    sumk:=0;
    for k:=0 to nf do begin
      f[k]:=exp(-k*k/s2);
      sumk:=sumk+f[k];
    end;
    sumk:=2*sumk-1; // jednu jednicku v nule odecist (aby tam nebyla dvakrat)

    sum:=0;
    for k:=-nf to nf do sum:=sum+dataIn[EnsureRange(i+k, 0, n1)]*f[Abs(k)];
    dataOut[i]:=sum/sumk;
  except
  end;
  SetLength(f, 0);
end;

procedure GaussFiltrAdapt(var data, sigma: TSingleArray); // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))
var
  dataNew: TSingleArray;
begin
  GaussFiltrAdapt(data, dataNew, sigma);
  data:=Copy(dataNew);
  SetLength(dataNew, 0);
end;

procedure GaussFiltrAdapt(var dataIn, dataOut, sigma: TDoubleArray); // filter dataIn using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2)) and save result to dataOut
var
  i, k: integer;
  f: TDoubleArray;
  n, n1, nf: integer;
  s2: double;
  sum, sumk: double;
begin
  n:=Length(dataIn);
  n1:=n-1;
  SetLength(dataOut, n);
  for i:=0 to n1 do try
    nf:=Round(sigma[i]*3+2);
    s2:=2*sqr(sigma[i]);
    SetLength(f, nf+1);
    sumk:=0;
    for k:=0 to nf do begin
      f[k]:=exp(-k*k/s2);
      sumk:=sumk+f[k];
    end;
    sumk:=2*sumk-1; // jednu jednicku v nule odecist (aby tam nebyla dvakrat)

    sum:=0;
    for k:=-nf to nf do sum:=sum+dataIn[EnsureRange(i+k, 0, n1)]*f[Abs(k)];
    dataOut[i]:=sum/sumk;
  except
  end;
  SetLength(f, 0);
end;

procedure GaussFiltrAdapt(var data, sigma: TDoubleArray); // filter data using convolution with Gauss kernel e^(- dx^2 / (2*sigma^2))
var
  dataNew: TDoubleArray;
begin
  GaussFiltrAdapt(data, dataNew, sigma);
  data:=Copy(dataNew);
  SetLength(dataNew, 0);
end;

end.
