File : lattice.adb

with Messages, Reps.Ops;
use Messages, Reps.Ops;

package body Lattice is

  function ToInt(I: Integer) return Int is
--- such sillyness should not be necessary
  begin
    if I >= 0 then
      return Int(I);
    else
      return Int(I-Integer'First)-(Int(-1-Integer'First)+1);
    end if;
  end ToInt;

  function Pt(I1,I2,I3: Integer) return IPoint is
  begin
    return (I1,I2,I3);
  end Pt;

  function Coord(L: IPoint) return RPoint is
    E: RPoint;
  begin
    E(1) := Rep(L(1)+L(2))*Half;
    E(2) := Rep(L(1)-L(2))*S3/Six;
    E(3) := Rep(L(3))*S6/Three;
    return E;
  end Coord;

  function RDist(L1,L2: IPoint) return Rep is
    R: constant RPoint := Coord(L2)-Coord(L1);
  begin
    return Sqrt(R*R);
  end RDist;

  function "*"(P,Q: IPoint) return Rep is
  begin
    return (P(1)*Q(2)+P(2)*Q(1)+2*(P(1)*Q(1)+P(2)*Q(2)+2*P(3)*Q(3)))/Six;
  end "*";

  procedure Translate(L: in IPoint; K: in out Int) is
  begin
    K := K+ToInt(L(1))+Fac1*(ToInt(L(2))+Fac1*ToInt(L(3)));
  end Translate;

  procedure Translate(L: in IPoint; B: in out IntVec) is
    D: constant Int := ToInt(L(1))+Fac1*(ToInt(L(2))+Fac1*ToInt(L(3)));
  begin
    for K in B'Range loop
      B(K) := B(K)+D;
    end loop;
  end Translate;

  procedure CountNbors(B: in BallPtr; N,Nfix: out Integer) is
  begin
    N := 0;
    Nfix := 0;
    if B=null then return; end if;
    for K in 1 .. 12 loop
      if B.Nb(K) /= null then
        N := N+1;
        if B.Nb(K).ID=0 then Nfix := Nfix+1; end if;
      end if;
    end loop;
  end CountNbors;

  procedure FreeAll(B: in out PolysPtr) is
  begin
    if B=null then return; end if;
    for K in B'Range loop
      if B(K) /= null then Free(B(K)); end if;
    end loop;
    Free(B);
  end FreeAll;

  procedure ResetParam(B: in PolyPtr) is
  begin
    B.L := 0; B.M := 1; B.R := Zero;
  end ResetParam;

  procedure CopyParam(B: in Poly; C: in out Poly) is
  begin
    C.L := B.L; C.M := B.M; C.R := B.R;
  end CopyParam;

  procedure CopyParam(B,C: in PolyPtr) is
  begin
    C.L := B.L; C.M := B.M; C.R := B.R;
  end CopyParam;

  function NumMovable(B: PolyPtr) return Natural is
  begin
    for K in 1 .. B.L loop
      if B.V(K) < Fac3 then return K-1; end if;
    end loop;
    return B.L;
  end NumMovable;

  function RDist(K: Int; B: PolyPtr) return Rep is
    P: constant IPoint := UnPack(K);
    D: Rep := Rep'Last;
  begin
    for N in 1 .. B.L loop
      D := Rep'Min(RDist(P,UnPack(B.V(N))),D);
    end loop;
    return D;
  end RDist;

  procedure AddSpace(N: in Positive; B: in out PolyPtr) is
  begin
    if B.L+N <= B.V'Last then return; end if;
    declare
      C: PolyPtr := new Poly(B.L+N);
    begin
      CopyParam(B,C);
      C.V(1 .. C.L) := B.V(1 .. C.L);
      Free(B);
      B := C;
    end;
  end AddSpace;

  procedure Translate(L: in IPoint; B: in Poly; C: out Poly) is
    D: constant Int := ToInt(L(1))+Fac1*(ToInt(L(2))+Fac1*ToInt(L(3)));
  begin
    CopyParam(B,C);
    for I in 1 .. B.L loop
      C.V(I) := B.V(I)+D;
    end loop;
  end Translate;

  function ">"(A,B: PolyPtr) return Boolean is
  begin
    for N in 1 .. Integer'Min(A.L,B.L) loop
      if A.V(N) /= B.V(N) then return A.V(N)>B.V(N); end if;
    end loop;
    return A.L>B.L;
  end ">";

  procedure Trim(B: in out PolysPtr; Size: in Natural) is
--- has undesired effects if B contains duplicate elements
    C: PolysPtr;
  begin
    if Size >= B'Last then return; end if;
    C := new Polys(1 .. Size);
    C.all(1 .. Size) := B.all(1 .. Size);
    for K in Size+1 .. B'Last loop
       Free(B(K));
    end loop;
    Free(B);
    B := C;
  end Trim;

  procedure SetDim(N1,N2,N3: in Positive) is
  begin
    Dim1 := N1; Dim2 := N2; Dim3 := N3;
    IZero := Pack((0,0,0));
  end SetDim;

  function InPatch(L: IPoint) return Boolean is
  begin
    return (abs(L(1)) <= Dim1) and (abs(L(2)) <= Dim2) and (abs(L(3)) <= Dim3);
  end InPatch;

  function FitsType(L: IPoint) return Boolean is
    Type12: constant Integer := (L(1)-L(2)) mod 3;
    Type3:  constant Integer := L(3) mod LayerPeriod;
  begin
    return Type3=Type12;
  end FitsType;

  function Check(I: Int) return Boolean is
    P: constant IPoint := UnPack(I);
  begin
    return InPatch(P) and then FitsType(P);
  end Check;

  function Pack(L: IPoint) return Int is
  begin
    return Int(Dim1+L(1))+Fac1*(Int(Dim2+L(2))+Fac1*Int(Dim3+L(3)));
  end Pack;

  procedure Pack(L: in IPoint; L4: in Natural; K: out Int) is
  begin
    K := Int(Dim1+L(1))+Fac1*(Int(Dim2+L(2))+Fac1*Int(Dim3+L(3)));
    if L4>0 then K := K+Fac3; end if;
  end Pack;

  function UnPack(K: Int) return IPoint is
    X: constant Integer := Integer( K       and Mask1)-Dim1;
    Y: constant Integer := Integer((K/Fac1) and Mask1)-Dim2;
    Z: constant Integer := Integer((K/Fac2) and Mask1)-Dim3;
  begin
    return (X,Y,Z);
  end UnPack;

  procedure UnPack(K: in Int; L: out IPoint; L4: out Natural) is
  begin
    L(1) := Integer( K       and Mask1)-Dim1;
    L(2) := Integer((K/Fac1) and Mask1)-Dim2;
    L(3) := Integer((K/Fac2) and Mask1)-Dim3;
    if K<Fac3 then L4 := 0; else L4 := 1; end if;
  end UnPack;

  procedure Copy(B: in Blob; P: in out PolyPtr) is
  begin
    if P=null then
      P := new Poly(B.L);
    elsif P.V'Last<B.L then
      Free(P);
      P := new Poly(B.L);
    else
      ResetParam(P);
    end if;
    P.L := B.L;
    for N in 1 .. P.L loop
      Pack(B.V(N).IP,B.V(N).ID,P.V(N));
    end loop;
  end Copy;

  procedure AddBallPlus(A: in Ball; B: in out PolyPtr) is
--- assumes ball has correct ID
    K: constant Integer := B.L+1;
  begin
    AddSpace(13,B);
    Pack(A.IP,A.ID,B.V(K));
    for N in 1 .. 12 loop
      Pack(A.Nb(N).IP,A.Nb(N).ID,B.V(K+N));
    end loop;
    B.L := K+12;
  end AddBallPlus;

  procedure FreeAll(E: in out PExpPtr) is
  begin
    for I in E'Range loop FreeAll(E(I)); end loop;
    Free(E);
  end FreeAll;

  function IDiff(A,B: PolyPtr) return Integer is
  begin
    for N in 1 .. Integer'Min(A.L,B.L) loop
      if A.V(N) /= B.V(N) then
        if A.V(N)>B.V(N) then return 1; end if;
        return -1;
      end if;
    end loop;
    if A.L>B.L then return 1; end if;
    if A.L<B.L then return -1; end if;
    return 0;
  end IDiff;

  procedure Find(P: in PolysPtr; B: in PolyPtr; Pos: out Integer; OK: out Boolean) is
--- assuming P is full and sorted
    L: Integer := P'First;
    R: Integer := P'Last;
    D: Integer;
  begin
    OK := True;
    loop
      Pos := (L+R)/2;
      D := IDiff(B,P(Pos));
      OK := D=0;
      if D>0 then
        if Pos=L then return; end if;
        R := Pos-1;
      elsif D<0 then
        if Pos=R then return; end if;
        L := Pos+1;
      else return; end if;
    end loop;
  end Find;

  procedure GetR(B: in PolyPtr; P: in PolysPtr; R: out Rep) is
    OK: Boolean;
    N: Integer;
  begin
    Find(P,B,N,OK);
    if not OK then Message("Lattice.GetR",Data_Problem); end if;
    R := P(N).R;
  end GetR;

  procedure PutR(B: in PolyPtr; P: in PolysPtr; R: in Rep) is
    OK: Boolean;
    N: Integer;
  begin
    Find(P,B,N,OK);
    if not OK then Message("Lattice.GetR",Data_Problem); end if;
    P(N).R := R;
  end PutR;

end Lattice;