// -----------------------------------------------------------------
// Hipercomplex functions
// for Quaternion Formula Compiler of the Fractal Explorer
// <c> 2004, Sirotinsky Arthur.
// =================================================================
// Please, do not modify this code without authors permission !
// =================================================================
Unit ComplexQ;

Interface
   uses Math;

type
   TQuat = array [0..3] of Single;


const
   QuatZ : TQuat = (0,0,0,0);
   Quat1 : TQuat = (1,0,0,0);
   SmalT = 1E-24;


Function FuncDisp(Func: Integer; qq: TQuat): TQuat;
Function QAdd(q1,q2: TQuat): TQuat;
Function QSub(q1,q2: TQuat): TQuat;
Function QSqr(q1: TQuat): TQuat;
Function QMul(q1,q2: TQuat): TQuat;
Function QRev(q1: TQuat): TQuat;
Function QDiv(q1,q2: TQuat): TQuat;
Function QDivR(q1: TQuat; a: Single): TQuat;
Function QSubR(q1: TQuat; a: Single): TQuat;
Function QMulR(q1: TQuat; a: Single): TQuat;
Function QExp(q: TQuat): TQuat;
Function QLog(q: TQuat): TQuat;
Function QPow(a,b: TQuat): TQuat;          // a^b
Function QSqrt(q: TQuat): TQuat;
Function QSqrt3(q: TQuat): TQuat;

Function QModule(q: TQuat): Single;



{ =================================================================== }
Implementation


Function QAdd(q1,q2: TQuat): TQuat;
Begin
  Result[0]:=q1[0]+q2[0];
  Result[1]:=q1[1]+q2[1];
  Result[2]:=q1[2]+q2[2];
  Result[3]:=q1[3]+q2[3];
End;

Function QSub(q1,q2: TQuat): TQuat;
Begin
  Result[0]:=q1[0]-q2[0];
  Result[1]:=q1[1]-q2[1];
  Result[2]:=q1[2]-q2[2];
  Result[3]:=q1[3]-q2[3];
End;

Function QSqr(q1: TQuat): TQuat;
  var tt: Single;
Begin
  tt:= q1[0]+q1[0];
  Result[0]:=q1[0]*q1[0] - q1[1]*q1[1] - q1[2]*q1[2] - q1[3]*q1[3];
  Result[1]:=q1[1]*tt;
  Result[2]:=q1[2]*tt;
  Result[3]:=q1[3]*tt;
End;

Function QMul(q1,q2: TQuat): TQuat;
Begin
  Result[0]:= q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3];   // real part
  Result[1]:= q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2];   // i part
  Result[2]:= q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1];   // j part
  Result[3]:= q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0];   // k part
End;

Function QRev(q1: TQuat): TQuat;
  var tt: Double;
Begin
  tt :=q1[0]*q1[0] + q1[1]*q1[1] + q1[2]*q1[2] + q1[3]*q1[3] + SmalT;         //  1/q1
  Result[0]:=q1[0]/tt;
  Result[1]:=q1[1]/tt;
  Result[2]:=q1[2]/tt;
  Result[3]:=q1[3]/tt;
End;

Function QDiv(q1,q2: TQuat): TQuat;
Begin
  Result:=QMul(q1, QRev(q2));
End;

Function QDivR(q1: TQuat; a: Single): TQuat;
Begin
  Result[0]:=q1[0]/a;
  Result[1]:=q1[1]/a;
  Result[2]:=q1[2]/a;
  Result[3]:=q1[3]/a;
End;


Function QSubR(q1: TQuat; a: Single): TQuat;
Begin
  Result:=q1;
  Result[0]:=q1[0]-a;
End;

Function QMulR(q1: TQuat; a: Single): TQuat;
Begin
  Result[0]:=a*q1[0];
  Result[1]:=a*q1[1];
  Result[2]:=a*q1[2];
  Result[3]:=a*q1[3];
End;

Function QExp(q: TQuat): TQuat;
  var n,e,f: Double;
Begin
  n:=Sqrt(q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
  e:=Exp(q[0]);
  f:=e*Sin(n);
  If n<>0 Then f:=f/n
          Else f:=0;

//	static double n, f, ex;
//	n = sqrt(q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
//	ex = exp(q[0]);
//	f = ex*sin(n);
//	if (n != 0.0) f /= n;

  Result[0]:=e*Cos(n);
  Result[1]:=f*q[1];
  Result[2]:=f*q[2];
  Result[3]:=f*q[3];
//      c[0] = ex*cos(n);
//	c[1] = f*q[1];
//	c[2] = f*q[2];
//	c[3] = f*q[3];
End;


Function QLog(q: TQuat): TQuat;
  var n,f: Double;
Begin
  n:=Sqrt(q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
  If n<>0 Then begin
     f := arctan2(n, q[0]);
     f := f/n;
     Result[0]:=0.5*Ln(q[0]*q[0] + n*n);
     Result[1]:=f*q[1];
     Result[2]:=f*q[2];
     Result[3]:=f*q[3];
  end
          Else begin
     Result[0]:=0.5*Ln(q[0]*q[0]);
     Result[1]:=arctan2(0, q[0]);
     Result[2]:=0;
     Result[3]:=0;
  end;
End;


Function QPow(a,b: TQuat): TQuat;          // a^b
  var an,bnp: Double;
      p: TQuat;
Begin
  an :=a[0]*a[0] + a[1]*a[1] + a[2]*a[2] + a[3]*a[3];
  bnp:=b[1]*b[1] + b[2]*b[2] + b[3]*b[3];
  If (Abs(an)<1E-100) and ( (b[0]>0) or (bnp<>0) ) Then begin
     Result[0]:= 0;
     Result[1]:= 0;
     Result[2]:= 0;
     Result[3]:= 0;
  end
          Else begin
     p:=QLog(a);
     p:=QMul(p,b);
     Result:=QExp(p);
  end;
End;


Function QSqrt(q: TQuat): TQuat;
Begin
  Result:=QLog(q);
  Result[0]:=Result[0]*0.5;
  Result[1]:=Result[1]*0.5;
  Result[2]:=Result[2]*0.5;
  Result[3]:=Result[3]*0.5;
  Result:=QExp(q);
End;

Function QSqrt3(q: TQuat): TQuat;
Begin
  Result:=QLog(q);
  Result[0]:=Result[0]*0.3333333;
  Result[1]:=Result[1]*0.3333333;
  Result[2]:=Result[2]*0.3333333;
  Result[3]:=Result[3]*0.3333333;
  Result:=QExp(q);
End;


{ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ }
Function FuncDisp(Func: Integer; qq: TQuat): TQuat;
Begin
  If Abs(qq[0])>1E6 Then qq[0]:=1E6;    // -do not change or delete these lines !!!
  If Abs(qq[1])>1E6 Then qq[1]:=1E6;    //
  If Abs(qq[2])>1E6 Then qq[2]:=1E6;    //
  If Abs(qq[3])>1E6 Then qq[3]:=1E6;    //
  Case func of
    0: Result:=qq;           // QIdent(real,imag);
    1: Result:=QuatZ;        // QZero
    2: Result:=Quat1;        // QOne
    3: begin                 // QAbs1
         Result:=qq;
         Result[0]:=Abs(Result[0]);
       end;
    4: begin                 // QAbs2
         Result[0]:=Abs(qq[0]);
         Result[1]:=Abs(qq[1]);
         Result[2]:=Abs(qq[2]);
         Result[3]:=Abs(qq[3]);
       end;
    5: Result:=QRev(qq);     // QRecip
    6: Result:=QSqrt(qq);    // QSqrt
    7: Result:=QSqrt3(qq);   // QSqrt3
    8: Result:=QSqr(qq);     // QSqr
    9: begin                 // QTriple
         Result:=QSqr(qq);
         Result:=QMul(Result,qq);
       end;
   10: Result:=QSqr(QSqr(qq));  // QFour
   11: Result:=QExp(qq);     // QExp
   12: Result:=QLog(qq);     // QLog
  end;
  If Abs(Result[0])>1E6 Then Result[0]:=1E6;    // -do not change or delete these lines !!!
  If Abs(Result[1])>1E6 Then Result[1]:=1E6;    //
  If Abs(Result[2])>1E6 Then Result[2]:=1E6;    //
  If Abs(Result[3])>1E6 Then Result[3]:=1E6;    //
End;



Function QModule(q: TQuat): Single;
Begin
  Result:= q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3];
End;

Begin
END.
{ =================================================================== }
{ =================================================================== }
