Code for Computing Dual Graphs of Special Fibers of Shimura Curves

The following is MAGMA code for computing the dual graph of the special fiber of the Shimura curve over QQ of discriminant Dp and level N (here DNp is squarefree and D is divisible by an odd number of primes). Of course this is heavily based on the Echidna/Magma code of David Kohel.
As with any code, buyer beware.
function LocalDualGram(L,q)
  Mat := MatrixAlgebra(ResidueClassRing(q), Rank(L));
  U := Kernel(Mat!GramMatrix(L));
  B := [ L!v : v in Basis(U) ] cat [ q*v : v in Basis(L) ];
  return MinkowskiGramReduction((1/q)*GramMatrix(sub< L | B >):
  Canonical := true);
end function;

function fiber(D,N,p)
  Q := QuaternionAlgebra(D);
  O := QuaternionOrder(Q,N);
  L := LeftIdealClasses(O);
  OO := QuaternionOrder(Q,N*p);
  LL := LeftIdealClasses(OO);
  OriginList := [ Explode([Index(L,I): I in L | IsIsomorphic(I,lideal)]) : II in LL];
  ALList := [0 : II in LL];
  OrderGrams := [ ReducedGramMatrix(RightOrder(II)) : II in LL];
  for A in SequenceToSet(OrderGrams) do
   Lat := LatticeWithGram(A);
   K := [ j : j in [1..#OrderGrams] | OrderGrams[j] eq A ];
   B := LocalDualGram(Lat,p);
   for i in K do

   ALList[i] := Explode([ j : j in K | IsIsometric(LatticeWithGram(1/(Norm(I)*Norm(J)) * GramMatrix(Conjugate(J)*I)),LatticeWithGram(B)) where I := LL[i] where J := LL[j] ]);
   end for;
  end for;
  W := ZeroMatrix(Integers(),#LL,#LL);
  for i in [1..#LL] do
   W[i,ALList[i]] := 1;
  end for;
  TerminusList := RowSequence(Transpose(W*Matrix(Integers(),   #LL,1,OriginList)))[1];
  ALOperators := AssociativeArray();
  for l in PrimeDivisors(D*N) do
   ALList := [0 : II in LL];
   for A in SequenceToSet(OrderGrams) do
    Lat := LatticeWithGram(A);
    K := [ j : j in [1..#OrderGrams] | OrderGrams[j] eq A ];
    B := LocalDualGram(Lat,l);
    for i in K do
    ALList[i] := Explode([ j : j in K | IsIsometric(LatticeWithGram(1/(Norm(I)*Norm(J)) * GramMatrix(Conjugate(J)*I)),LatticeWithGram(B)) where I := LL[i] where J := LL[j] ]);
    end for;
   end for;
   ALOperators[l] := ALList;
  end for;
  return OriginList , TerminusList, [ExactQuotient(#Units(RightOrder(II)),2) : II in LL], ALOperators;
end function;