Базовые функции (Pascal)
 

type 
     Complex = record
     Re,Im: real;
     end;

const
     Eps: real = 0.01;

 
  Функция определения знака вещественной величины

function Sign(R: real): integer;

begin
     Result:=1;
     if R<0 then Result:=-1;
end; // Sign

 
  Функция создания комплексного числа из вещественных компоннтов - действительной и мнимой частей

function MCompl(X, Y: real): сomplex;
var Z: complex;

begin
     if Abs(X)<Eps*Eps*Eps then X:=0;
     if Abs(Y)<Eps*Eps*Eps then Y:=0;
     Compl(Z,X,Y);
     MCompl:=Z;
end; // MComp

 

procedure Compl(var Z: complex;; X,Y: real);

begin
     if Abs(X)<1E-6 then X:=0;
     if Abs(Y)<1E-6 then Y:=0;
     Z.Re:=X;
     Z.Im:=Y;
end; // Compl

 
  Функция создания отрицательного комплексного числа

function CompNeg(X: complex;): complex;

begin
     Result:=MCompl(-X.Re,-X.Im);
end; // CompNeg

 
  Функция сложения двух комплексных чисел

function CompSum(X,Y: complex;): complex;;
var Z: complex;

begin
     Z.Re:=X.Re+Y.Re;
     Z.Im:=X.Im+Y.Im;
     Result:=Z;
end; // CompSum

 
  Функция разности двух комплексных чисел

function CompSub(X,Y: complex;): complex;;
var Z: complex;

begin
     Z.Re:=X.Re-Y.Re;
     Z.Im:=X.Im-Y.Im;
     Result:=Z;
end; // CompSub

 
  Функция произведения двух комплексных чисел

function CompMul(X,Y: complex): complex;
var Z: complex;

begin
     Z.Re:=X.Re*Y.Re-X.Im*Y.Im;
     Z.Im:=X.Re*Y.Im+Y.Re*X.Im;
     Result:=Z;
end; // CompMul

 
  Процедура деления двух комплексных чисел

procedure CompDiv(X,Y: complex; var Z: complex; var Success: boolean);

begin
     Success:=TRUE;
     if Abs (Sqr(Y.Re)+Sqr(Y.Im))<Eps then
     begin
          Success:=FALSE;
     end;
     Z.Re:=(X.Re*Y.Re+X.Im*Y.Im)/(Sqr(Y.Re)+Sqr(Y.Im));
     Z.Im:=(Y.Re*X.Im-X.Re*Y.Im)/(Sqr(Y.Re)+Sqr(Y.Im));
end; // CompDiv

 
  Функция деления двух комплексных чисел

function CompDiv1(X,Y: complex): complex;
var Z: complex;
    Success: boolean;
begin

     Success:=TRUE;
     if Abs (Sqr(Y.Re)+Sqr(Y.Im))<Eps then
     begin
          Success:=FALSE;
     end;
     Z.Re:=(X.Re*Y.Re+X.Im*Y.Im)/(Sqr(Y.Re)+Sqr(Y.Im));
     Z.Im:=(Y.Re*X.Im-X.Re*Y.Im)/(Sqr(Y.Re)+Sqr(Y.Im));
     Result:=Z;
end; // CompDiv1

 
  Функция возведения комплексного числа во вторую степень

function CompSqr(X: complex): complex;

begin
     Result:=CompMul(X,X);
end; // CompSqr

 
  Функция извлечения квадратного корня из комплексного числа

function CompSqrt(X: complex): complex;
var A,B: real;

begin
     A:=Sqrt((Sqrt(X.Re*X.Re+X.Im*X.Im)+X.Re)/2);
     B:=Sign(X.Im)*Sqrt((Sqrt(X.Re*X.Re+X.Im*X.Im)-X.Re)/2);
     Result:=MCompl(A,B);
end; // CompSqrt

 
  Функция вычисления вещественного дискриминанта третьего порядка

function Diskr(P1,P2,P3,P4,P5,P6,P7,P8,P9: real): real;

begin
     Diskr:=P1*P5*P9+P4*P8*P3+P7*P2*P6-P7*P5*P3-P1*P8*P6-P4*P2*P9;
end; // Diskr

 
  Функция вычисления комплексного дискриминанта третьего порядка

function DiskrComp(P1,P2,P3,P4,P5,P6,P7,P8,P9: complex): complex;
var S1,S2,S3,S4,S5,S6,D1,D2: complex;

begin
     S1:=CompMul(P1,CompMul(P5,P9));
     S2:=CompMul(P4,CompMul(P8,P3));
     S3:=CompMul(P7,CompMul(P2,P6));
     S4:=CompMul(P7,CompMul(P5,P3));
     S5:=CompMul(P1,CompMul(P8,P6));
     S6:=CompMul(P4,CompMul(P2,P9));
     D1:=CompSum(S1,CompSum(S2,S3));
     D2:=CompSum(S4,CompSum(S5,S6));
     Result:=CompSub(D1,D2);
end; // DiskrComp

 
  Функция вычисления синуса и косинуса направленного отрезка

procedure SС(X1,Y1,X2,Y2: real; var S1,C1,Dl: real; var Prizn: boolean);

var Dx,Dy: real;
begin
     Prizn:=FALSE;
     Dx:=X2-X1; Dy:=Y2-Y1; Dl:=Sqrt(Sqr(Dx)+Sqr(Dy));
     if Dl>0 then
     begin
          S1:=Dy/Dl;
          C1:=Dx/Dl;
          Prizn:=TRUE;
     end;
end; // Sc

 
  Функция вычисления арксинуса вещественного числа

function ArcSin(X: extended): extended;

var A: extended;
begin
     if(X>-1) and (X<1) then A:=ArcTan(X/Sqrt(1-Sqr(X)));
     if(X>1-Sqr(Eps*Eps))then A:=PI/2;
     if(X<-1+Sqr(Eps*Eps))then A:=-PI/2;
     ArcSin:=A;
end; // ArcSin

 
  Функция вычисления арккосинуса вещественного числа

function ArcCos(X: extended): extended;

var A: real;
begin
     ArcCos:=Pi/2-ArcSin(X);
end; // ArcCos

 
  Функция вычисления величины угла между двумя прямыми линиями, заданными направляющими синусами и косинусами

procedure Fi(Sa,Ca,Sb,Cb: real; var Df: real);
var Da1,Da2,Fi1,Fi2: real;

begin
     Fi2:=ArcSin(Sb);
     if (Sa>=0) and (Ca>=0) then Fi1:=ArcSin(Sa);
     if (Sa>=0) and (Ca<0) then Fi1:=Pi-ArcSin(Sa);
     if (Sa<0) and (Ca>=0) then Fi1:=ArcSin(Sa);
     if (Sa<0) and (Ca<0) then Fi1:=-Pi-ArcSin(Sa);
     if (Sb>=0) and (Cb>=0) then Fi2:=ArcSin(Sb);
     if (Sb>=0) and (Cb<0) then Fi2:=Pi-ArcSin(Sb);
     if (Sb<0) and (Cb>=0) then Fi2:=ArcSin(Sb);
     if (Sb<0) and (Cb<0) then Fi2:=-Pi-ArcSin(Sb);
     Df:=Fi2-Fi1;
end; // Fi

 
  Функция определения координат точки пересечения двух вещественных прямых линий, заданных координатами пар несовпадающих точек

function LinLin(X1,Y1,X2,Y2,X3,Y3,X4,Y4: real; var XX,YY,ZZ: real): boolean;
var A1,B1,C1,A2,B2,C2: real;
      Prizn: boolean;

begin
     Prizn:=FALSE;
     A1:=Y1-Y2; B1:=X2-X1;
     A2:=Y3-Y4; B2:=X4-X3;
     ZZ:=A1*B2-B1*A2;

     if Abs(ZZ)>1E-2 then
     begin
          C1:=X1*Y2-Y1*X2;
          C2:=X3*Y4-Y3*X4;
          XX:=(B1*C2-C1*B2)/ZZ;
          YY:=(C1*A2-A1*C2)/ZZ;
          ZZ:=1;
          Prizn:=TRUE
     end;

     if (Abs(ZZ)<=1E-2) or (Abs(XX)>1000000000) or (Abs(YY)>1000000000) then
     begin
          XX:=B1;
          YY:=-A1;
          ZZ:=0;
     end;

     Result:=Prizn;
end; // LinLin

 

 

function ProcA(X1,Y1,X2,Y2: real; var Xt,Yt: real; V: real; var XD: real):  boolean;
var S,C,Dl,Dx,Xt1,Yt1,Df,X3,Y3,X4,Y4,A,B: real;
    LDummy,Prizn: boolean;

begin
     ProcA:=FALSE;
     SC(X1,Y1,X2,Y2,S,C,Dl,LDummy);
     Fi(0,1,S,C,Df);
     Dx:=90*Pi/180+Df;
     X3:=Xt+200*Cos(Dx);
     Y3:=Yt+200*Sin(Dx);
     X4:=Xt-200*Cos(Dx);
     Y4:=Yt-200*Sin(Dx);
     Prizn:=LinLin(X1,Y1,X2,Y2,X3,Y3,X4,Y4,Xt1,Yt1,C);
     A:=Sqr(Xt-Xt1); B:=Sqr(Yt-Yt1);
     XD:=Sqrt(A+B);

     if XD< V then
     begin
          ProcA:=TRUE;
          Xt:=Xt1;Yt:=Yt1;
     end;
end; // ProcA

 
  Функция определения координат точки пересечения двук комплексных прямых линий

function LinLinComp(X1,Y1,X2,Y2,X3,Y3,X4,Y4: complex; var XX,YY,ZZ:  complex): boolean;
var A1,B1,C1,A2,B2,C2: complex;
    B,BB,BBB,Prizn: boolean;
    V,XD: real;
label fin;

begin

     Prizn:=FALSE;

     if (X1.Im=0) and (X2.Im=0) and (X3.Im=0) and (X4.Im=0) and (Y1.Im=0) and (Y2.Im=0) and (Y3.Im=0) and (Y4.Im=0) then
     BEGIN
          V:=0.0001;
          BB:=ProcA(X1.Re,Y1.Re,X2.Re,Y2.Re,X3.Re,Y3.Re,V,XD);
          BBB:=ProcA(X1.Re,Y1.Re,X2.Re,Y2.Re,X4.Re,Y4.Re,V,XD);
          if BB and BBB then
          begin
               B1:=CompSub(X2,X1);A1:=CompSub(Y1,Y2);
               XX:=B1;
               YY:=CompSub(MCompl(0,0),A1);
               ZZ:=MCompl(0,0);
               goto fin;
          end;
     END;

     A1:=CompSub(Y1,Y2); B1:=CompSub(X2,X1);
     A2:=CompSub(Y3,Y4); B2:=CompSub(X4,X3);
     ZZ:=CompSub(CompMul(A1,B2),CompMul(B1,A2));

     if (Abs(ZZ.Re)>1E-25) or (Abs(ZZ.Im)>1E-25) then
     begin
          C1:=CompSub(CompMul(X1,Y2),CompMul(Y1,X2));
          C2:=CompSub(CompMul(X3,Y4),CompMul(Y3,X4));
          CompDiv(CompSub(CompMul(B1,C2),CompMul(C1,B2)),ZZ,XX,B);
          CompDiv(CompSub(CompMul(C1,A2),CompMul(A1,C2)),ZZ,YY,B);
          ZZ:=MCompl(1,0);
          Prizn:=TRUE
     end;

     if (Abs(ZZ.Re)<=1E-25) or (Abs(XX.Re)>1000000000) or (Abs (YY.Re)>1000000000) then
     begin
          XX:=B1;
          YY:=CompSub(MCompl(0,0),A1);
          ZZ:=MCompl(0,0);
     end;

     fin:
     Result:=Prizn;

end; // LinLinComp

 

procedure ComDuga(Dx,Dy,Xc1,Yc1,Z,R1: real; var X1,X2,Y1,Y2: complex; var  ComplOut: boolean);
var A,B,C,D,Sqd,Dy2,Dx2: real;
label fin;

begin
     ComplOut:=FALSE;
     if Abs(Dy)>Abs(Dx) then
     begin
          Dy2:=2*Dy;Dx2:=2*Dx;
          A:=(1+Sqr(Dx)/Sqr(Dy))*2;
          B:=-2*Xc1+Z*Dx/Sqr(Dy)-Dx2*Yc1/Dy;
          C:=Sqr(Xc1)+Sqr(Z)/4/Sqr(Dy)-Z*Yc1/Dy+Sqr(Yc1)-Sqr(R1);
          D:=Sqr(B)-2*A*C;
          if Abs(D)<0.00000001 then D:=0;
          if D>=0 then
          begin
               Sqd:=Sqrt(D);
               X1.Re:=(-B+Sqd)/A;
               X1.Im:=0;
               X2.Re:=(-B-Sqd)/A;
               X2.Im:=0;
               Y1.Re:=(Z+X1.Re*Dx2)/Dy2;
               Y1.Im:=0;
               Y2.Re:=(Z+X2.Re*Dx2)/Dy2;
               Y2.Im:=0;
          end else
          begin
               X1.Re:=-B/A; X2.Re:=-B/A;
               X1.Im:=Sqrt(Abs(D))/A; X2.Im:=-Sqrt(Abs(D))/A;
               Y1.Re:=(Z+(-B/A)*Dx2)/Dy2; Y1.Im:=Sqrt(Abs(D))/A*Dx2/Dy2;
               Y2.Re:=(Z+(-B/A)*Dx2)/Dy2; Y2.Im:=-Sqrt(Abs(D))/A*Dx2/Dy2;
               ComplOut:=TRUE;
          end;
     end
     else
     begin
          Dx:=-Dx; Dy:=-Dy;
          Dy2:=2*Dy;Dx2:=2*Dx;
          A:=(1+Sqr(Dy)/Sqr(Dx))*2;
          B:=-2*Yc1+Z*Dy/Sqr(Dx)-Dy2*Xc1/Dx;
          C:=Sqr(Yc1)+Sqr(Z)/4/Sqr(Dx)-Z*Xc1/Dx+Sqr(Xc1)-Sqr(R1);
          D:=Sqr(B)-2*A*C;
          if D>=0 then
          begin
               Sqd:=Sqrt(D);
               Y1.Re:=(-B+Sqd)/A;
               Y1.Im:=0;
               Y2.Re:=(-B-Sqd)/A;
               Y2.Im:=0;
               X1.Re:=(Z+Y1.Re*Dy2)/(Dx2);
               X1.Im:=0;
               X2.Re:=(Z+Y2.Re*Dy2)/(Dx2);
               X2.Im:=0;
          end else
          begin
               Y1.Re:=-B/A; Y2.Re:=-B/A;
               Y1.Im:=Sqrt(Abs(D))/A; Y2.Im:=-Sqrt(Abs(D))/A;
               X1.Re:=(Z+(-B/A)*Dy2)/Dx2; X1.Im:=Sqrt(Abs(D))/A*Dy2/Dx2;
               X2.Re:=(Z+(-B/A)*Dy2)/Dx2; X2.Im:=-Sqrt(Abs(D))/A*Dy2/Dx2;
               ComplOut:=TRUE;
          end;
     end;


fin:
end; // ComDuga