File : minifuns-ops.adb

with Ints, Messages;
use Ints, Messages;

pragma Elaborate (Ints);

package body MiniFuns.Ops is

  T3: constant Integer := 3*Int(M_Trunc);

  function Det(B: Mat3) return Scalar is
  begin
    return B.U.A*B.V.B-B.U.B*B.V.A;
  end Det;

  procedure Normalize(B: in out Bas3) is
--- want B.U.A=1 and B.V.B=1
  begin
    B.U.B := B.U.B/B.U.A;
    B.U.C := B.U.C/B.U.A;
    B.U.A := Scal_One;
    B.V.A := B.V.A/B.V.B;
    B.V.C := B.V.C/B.V.B;
    B.V.B := Scal_One;
  end Normalize;

  procedure Translate(R0,S0: in Scalar; F: in out Vec3) is
  begin
    F.C := F.C-R0*F.A-S0*F.B;
  end Translate;

  procedure Prod(B: in Mat3; W1: in Vec3; W2: out Vec3) is
  begin
    W2.A := W1.A*B.U.A+W1.B*B.V.A;
    W2.B := W1.A*B.U.B+W1.B*B.V.B;
    W2.C := W1.A*B.U.C+W1.B*B.V.C+W1.C;
  end Prod;

  procedure Prod(B1,B2: in Mat3; B3: out Mat3) is
  begin
    Prod(B1,B2.U,B3.U);
    Prod(B1,B2.V,B3.V);
  end Prod;

  procedure Inv(B: in out Mat3) is
    IDet: constant Scalar := Inv(Det(B));
    M: constant Mat3 := B;
  begin
    B.U.A :=  IDet*M.V.B;
    B.U.B := -IDet*M.U.B;
    B.V.A := -IDet*M.V.A;
    B.V.B :=  IDet*M.U.A;
    B.U.C := -(B.U.A*M.U.C+B.U.B*M.V.C);
    B.V.C := -(B.V.A*M.U.C+B.V.B*M.V.C);
  end Inv;

  procedure Scale(S: in Scalar; B: in out Bas3) is
  begin
    B.U.A := Sqr(S)*B.U.A;
    B.U.B := S*B.U.B;
    B.V.A := Sqr(S)*B.V.A;
    B.V.B := S*B.V.B;
  end Scale;

  procedure Center(B: in out Bas3; R0,S0: out Scalar) is
--- eliminate constants
  begin
    if B.U.C=Scal_Zero and then B.V.C=Scal_Zero then
      R0 := Scal_Zero;
      S0 := Scal_Zero;
    else
      S0 := Inv(Det(B));
      R0 := (B.U.B*B.V.C-B.V.B*B.U.C)*S0;
      S0 := (B.V.A*B.U.C-B.U.A*B.V.C)*S0;
      B.U.C := Scal_Zero;
      B.V.C := Scal_Zero;
    end if;
  end Center;

  function Dom(U,V: in Vec3; R: in RadPair) return RadPair is
    --- functions A*U+B*V+C
    UA:  constant Scalar := Scal(MaxNorm(U.A));
    UB:  constant Scalar := Scal(MaxNorm(U.B));
    VA:  constant Scalar := Scal(MaxNorm(V.A));
    VB:  constant Scalar := Scal(MaxNorm(V.B));
    Det: constant Scalar := UA*VB-UB*VA;
    Err: Boolean;
    D1,D2: Rep;
    N1,N2: Scalar;
  begin
    D1 := Inf(Scal(R(1))-Norm(U.C));
    D2 := Inf(Scal(R(2))-Norm(V.C));
    Err := (D1<Zero) or (D2<Zero);
    N1 := (D1*VB-D2*UB)/Det;
    N2 := (D2*UA-D1*VA)/Det;
    D1 := Inf(N1);
    D2 := Inf(N2);
    Err := Err or (D1<Zero) or (D2<Zero);
    if Err then
      Error("MiniFuns.Ops.Dom error: domain violation",T3);
      if M_Trunc then return R; end if;
    end if;
    return (D1,D2);
  end Dom;

  procedure Show1(N: in String; W: in Vec3) is
  begin
    Show1(N & "A =",W.A);
    Show1(N & "B =",W.B);
    Show1(N & "C =",W.C);
    New_Line;
  end Show1;

  procedure Show1(N: in String; B: in Mat3) is
  begin
    Show1(N & "U.",B.U);
    Show1(N & "V.",B.V);
  end Show1;

----------------------------------------

  function InvDet(A: Args) return Scalar is
  begin
    return Inv(Scal_One-A.A*Scal(A.B));
  end InvDet;

  procedure Approx(B: in Bas3; A: out Args) is
    R0,S0: Scalar;
    Bt: Bas3 := B;
  begin
    Normalize(Bt);
    Center(Bt,R0,S0);
    A.S0 := Approx(S0);
    A.T0 := Approx(Sqrt(R0));
    A.B := Approx(Bt.U.B);
    A.A := Approx(Bt.V.A);
  end Approx;

  procedure Convert(A: in Args; B: out Bas3) is
    R0: constant Scalar := A.T0*Scal(A.T0);
    S0: constant Scalar := Scal(A.S0);
  begin
    B.U.A := Scal_One;
    B.U.B := Scal(A.B);
    B.U.C := Scal_Zero;
    B.V.A := Scal(A.A);
    B.V.B := Scal_One;
    B.V.C := Scal_Zero;
    Translate(R0,S0,B.U);
    Translate(R0,S0,B.V);
  end Convert;

  procedure Coord(A1,A2: in Args; B: out Mat3) is
--- coordinates of A1 in "basis" A2
    B1,B2: Mat3;
  begin
    Convert(A1,B1);
    Convert(A2,B2);
    Inv(B2);
    Prod(B2,B1,B);
  end Coord;

  function CoordT2(A: Args) return Vec3 is
    W: Vec3;
  begin
    W.A := InvDet(A);
    W.B := (-A.B)*W.A;
    W.C := A.T0*Scal(A.T0);
    return W;
  end CoordT2;

  function CoordS1(A: Args) return Vec3 is
    W: Vec3;
  begin
    W.B := InvDet(A);
    W.A := (-A.A)*W.B;
    W.C := Scal(A.S0);
    return W;
  end CoordS1;

  procedure Scale(S: in Scalar; A1,A2: in Args; B: out Mat3) is
--- coordinates of scaled A1 in basis A2
    B1,B2: Mat3;
  begin
    Convert(A1,B1);
    Scale(S,B1);
    Convert(A2,B2);
    Inv(B2);
    Prod(B2,B1,B);
  end Scale;

------------------------------------

  procedure ValTS(T,S: in Scalar; A: in Args; U,V: out Scalar) is
    T0: constant Scalar := Scal(A.T0);
    S0: constant Scalar := Scal(A.S0);
  begin
    V := A.A*(T+T0)*(T-T0)+(S-S0);
    U := (T+T0)*(T-T0)+A.B*(S-S0);
  end ValTS;

  procedure Val(X,Y: in Scalar; A: in Args; U,V: out Scalar) is
  begin
    ValTS(X+Y,X-Y,A,U,V);
  end Val;

  function ValTSRadii(T,S: Scalar; A: Args) return RadPair is
    U,V: Scalar;
  begin
    ValTS(T,S,A,U,V);
    return (MaxNorm(U),MaxNorm(V));
  end ValTSRadii;

  function ValRadii(X,Y: Scalar; A: Args) return RadPair is
  begin
    return ValTSRadii(X+Y,X-Y,A);
  end ValRadii;

  function ValS1(U,V: Scalar; A: Args) return Scalar is
  begin
    return Scal(A.S0)+(V-A.A*U)*InvDet(A);
  end ValS1;

  function ValT2(U,V: Scalar; A: Args) return Scalar is
  begin
    return A.T0*Scal(A.T0)+(U-A.B*V)*InvDet(A);
  end ValT2;

  function TNorm(A: Args; R: RadPair) return Scalar is
    U: constant Scalar := Scal(R(1));
    V: constant Scalar := Scal(R(2));
  begin
    return Sqrt(A.T0*Scal(A.T0)+(U+abs(A.B)*V)*InvDet(A));
  end TNorm;

  function TNorm(D: Domain) return Scalar is
  begin
    return TNorm(D.A,D.R);
  end TNorm;

  procedure CheckDomainTS(E: in String; D: in Domain; T,S: in Scalar) is
    U,V: Scalar;
  begin
    ValTS(T,S,D.A,U,V);
    if MaxNorm(U)>D.R(1) then
      Error(E & " error: |U|>Rho1",MaxNorm(U),D.R(1),T3);
    end if;
    if MaxNorm(V)>D.R(2) then
      Error(E & " error: |V|>Rho2",MaxNorm(V),D.R(2),T3);
    end if;
  end CheckDomainTS;

  procedure CheckDomain(E: in String; D: in Domain; X,Y: in Scalar) is
  begin
    CheckDomainTS(E,D,X+Y,X-Y);
  end CheckDomain;

end MiniFuns.Ops;