PROCEDURE InterpolateTriangD3C0
    (* Compute the vertex-associated coefficients: *)
      
    FOR i := 0 TO 2 DO 
      c0[CoeffIndex(6,0,0, ie)] := 0.25d0 * f[ie];
      c1[CoeffIndex(5,0,0, ie)] := 0.75d0 * f[ie];
    END;

    FOR i := 0 TO 2 DO 
      c0[CoeffIndex(3,0,0, i)] := lf[i] - f[i];
      c1[CoeffIndex(2,0,0, i)] := 2.0d0 * f[i] - lf[i];
    END;

    c0[0] := -0.50d0 * f[0]; c1[0] := +1.50d0 * f[0];
    c0[6] := -0.50d0 * f[1]; c1[3] := +1.50d0 * f[1];
    c0[9] := -0.50d0 * f[2]; c1[5] := +1.50d0 * f[2];



    c0[0] := 0.25d0 * f[0]; c1[0] := 0.75d0 * f[0];
    c0[6] := 0.25d0 * f[1]; c1[3] := 0.75d0 * f[1];
    c0[9] := 0.25d0 * f[2]; c1[5] := 0.75d0 * f[2];


  CONST
    alpha = 0.1464466095d0;
    beta =  0.8535533905d0;



(*PROCEDURE Area(READONLY p, q, r: Point): LONGREAL;*)
  (*
    Area of the spherical triangle with vertices "p,q,r". *)

PROCEDURE Area(READONLY p, q, r: Point): LONGREAL =
  (*
    Area of the spatial triangle with vertices "p,q,r" *)
  CONST Pi = 3.14159265358979d0;
  BEGIN
    WITH
      A = LR3Extras.Cross(q,r),
      B = LR3Extras.Cross(r,p),
      C = LR3Extras.Cross(p,q),
      AB = acos(-LR3.Cos(A,B)),
      BC = acos(-LR3.Cos(B,C)),
      AC = acos(-LR3.Cos(A,C)),
      area = AB + BC + AC - Pi 
    DO
      RETURN area
    END
  END Area;

(* From SPPWFunction.Dotproduct: *)

 TYPE
  NAT = CARDINAL; 
    
   fO, gO: REF SPOverlapTable.T;

    In that case, "fO[i]"
    must be the list of numbers of triangles in "f.tri" that overlap
    triangle "i" of "tri". If "fO = NIL", the procedure requires that
    "f.tri = tri".  The analogous requirements hold for "gO" and
    "g.tri".
    
    IF fO = NIL THEN 
      <* ASSERT tri = f.tri *>
    ELSE
      <* ASSERT NUMBER(fO^) = NUMBER(tri.side^) *>
    END;
    IF gO = NIL THEN 
      <* ASSERT tri = g.tri *>
    ELSE
      <* ASSERT NUMBER(gO^) = NUMBER(tri.side^) *>
    END;
    
    WITH
      trivList = NEW(REF ARRAY OF NAT, 1)
    DO
      VAR fTs, gTs: REF ARRAY OF NAT;
      BEGIN

        trivList[0] := i;
        (* Find the triangles "fTs" of "f.tri" that overlap triangle "i" of "tri" *)
        IF fO = NIL THEN fTs := trivList ELSE fTs := fO[i] END;

        (* Find the triangles "gTs" of "g.tri" that overlap triangle "i" of "tri" *)
        IF gO = NIL THEN gTs := trivList ELSE gTs := gO[i] END;
      END;  



    END
    
    
    If "tri # f.tri", "fO[i]" must be the list of numbers of triangles
    in "f.tri" that overlap triangle "i" of "tri".  The analogous
    requirements hold for "gO" and "g.tri".

    RETURN DiffTriangIntegral(f.tri, g.tri, f, g, H, order, fZero, gZero)
       RETURN DiffTriangIntegral(f.tri, g.tri, f, g, H, order, TRUE, TRUE)
   
(* *)

    PROCEDURE Small(READONLY p, q, r: Point): BOOLEAN =
      BEGIN 
        WITH 
          dpq2 = LR3.DistSqr(p, q),
          dpr2 = LR3.DistSqr(p, r),
          dqr2 = LR3.DistSqr(q, r)
        DO 
          RETURN (dpq2 <= tol) AND (dpr2 <= tol) AND (dqr2 <= tol)
        END;
      END Small;
    


(* Old Area *)    
  
  
    WITH
      pq = LR3.Sub(q,p),
      pr = LR3.Sub(r,p),
      c = LR3Extras.Cross(pq, pr)
    DO
      RETURN 0.5d0 * LR3.Norm(c)
    END




PROCEDURE OnTriangle(
    READONLY p, q, r: SPTriang.Coords;
    func: Function;
    small: TrianglePredicate;
    depth: CARDINAL := LAST(CARDINAL);
  ): LONGREAL =
  CONST mu = 0.8d0; nu = 0.2d0;
  BEGIN
    (* PrintTriangle(depth, p, q, r); *)
    IF depth = 0 OR small(p,q,r) THEN
      WITH
        area = Area(p,q,r),
        pf = LR3.Dir(LR3.Mix(mu, p, nu, LR3.Mix(0.5d0, q, 0.5d0, r))),
        qf = LR3.Dir(LR3.Mix(mu, q, nu, LR3.Mix(0.5d0, r, 0.5d0, p))),
        rf = LR3.Dir(LR3.Mix(mu, r, nu, LR3.Mix(0.5d0, p, 0.5d0, q))),
        avg = (func(pf) + func(qf) + func(rf))/3.0d0
      DO
        RETURN avg * area
      END
    ELSE
      WITH
        pq = LR3.Dir(LR3.T{(p[0]+q[0])/2.0d0, (p[1]+q[1])/2.0d0, (p[2]+q[2])/2.0d0}),
        rq = LR3.Dir(LR3.T{(r[0]+q[0])/2.0d0, (r[1]+q[1])/2.0d0, (r[2]+q[2])/2.0d0}),
        pr = LR3.Dir(LR3.T{(p[0]+r[0])/2.0d0, (p[1]+r[1])/2.0d0, (p[2]+r[2])/2.0d0}),
        I1 = OnTriangle(q, pq, rq, func, small, depth-1),
        I2 = OnTriangle(p, pq, pr, func, small, depth-1),
        I3 = OnTriangle(r, pr, rq, func, small, depth-1),
        I4 = OnTriangle(pr, pq, rq, func, small, depth-1)
      DO
        RETURN I1+I2+I3+I4
      END
    END
  END OnTriangle;

PROCEDURE OnTriangle(
    READONLY p, q, r: Point;
    func: Function;
    order: CARDINAL;
  ): LONG =
  VAR 
    sf, cf, 
    sw, cw, 
    xf, xw, 
    cor: LONGREAL:= 0.0d0;
    CONST Pi = 3.14159265358979d0;  
  BEGIN
    WITH
      mm = order,
      fmm = FLOAT(mm, LONG),
      n = LR3.Dir(LR3Extras.Cross(LR3.Sub(q, p), LR3.Sub(r, p))),
      A = LR3Extras.Cross(q,r),
      B = LR3Extras.Cross(r,p),
      C = LR3Extras.Cross(p,q),
      AB = acos(-LR3.Cos(A,B)),
      BC = acos(-LR3.Cos(B,C)),
      AC = acos(-LR3.Cos(A,C)),
      area = AB + BC + AC - Pi
    DO
    FOR k:= 0 TO mm DO 
        FOR j:= 0 TO mm-k DO 
          WITH
            i = mm-(k+j),
            fi = FLOAT(i,LONG), fj = FLOAT(j,LONG), fk = FLOAT(k,LONG),
            a = LR3.T{
              (p[0]*fi + q[0]*fj + r[0]*fk)/fmm,
              (p[1]*fi + q[1]*fj + r[1]*fk)/fmm,
              (p[2]*fi + q[2]*fj + r[2]*fk)/fmm
            },
            t = LR3.Dir(a),
            sqra = LR3.NormSqr(a),
            na = LR3.Dot(n, t)
         DO
            xw := na/sqra;
            IF i # 0 AND j # 0 AND k # 0 THEN 
              cor := 1.0d0
            ELSIF (i = 0 AND j = 0) THEN cor := AB/(2.0d0*Pi)
            ELSIF (i = 0 AND k = 0) THEN cor := AC/(2.0d0*Pi)
            ELSIF (j = 0 AND k = 0) THEN cor := BC/(2.0d0*Pi)
            ELSIF i = 0 OR j = 0 OR k = 0 THEN
              cor:= 1.0d0/2.0d0         
            END;
            xw := cor*xw;
            xf := xw*func(t);
            WITH 
              (* Kahan's summation formula: *)
              yf = xf - cf,
              sNewf = sf + yf,
              yw = xw - cw,
              sNeww = sw + yw
            DO  
              cf := (sNewf - sf) - yf;
              sf := sNewf;
              cw := (sNeww - sw) - yw;
              sw := sNeww;
            END;
          END;   
        END;
      END;
      RETURN (area/sw)*sf;  
    END;
  END OnTriangle;
  
PROCEDURE OnTriangle(
    READONLY p, q, r: Point;
    func: Function;
    order: CARDINAL;
  ): LONG =
  VAR 
    s, c: LONGREAL:= 0.0d0;
    np: CARDINAL := 0;
  BEGIN
    WITH
      mm = 2*order + 1,
      fmm = FLOAT(mm, LONG)
    DO
      FOR k:= 1 TO mm-2 BY 2 DO 
        FOR j:= 1 TO mm-(k+1) BY 2 DO 
          WITH
            i = mm-(k+j),
            fi = FLOAT(i,LONG), fj = FLOAT(j,LONG), fk = FLOAT(k,LONG),
            a = LR3.T{
              (p[0]*fi + q[0]*fj + r[0]*fk)/fmm,
              (p[1]*fi + q[1]*fj + r[1]*fk)/fmm,
              (p[2]*fi + q[2]*fj + r[2]*fk)/fmm
            },
            t = LR3.Dir(a),
            x = func(t),
            (* Kahan's summation formula: *)
            y = x - c,
            sNew = s + y
          DO
            c := (sNew - s) - y;
            s := sNew
          END;
          INC(np)
        END;
      END
      END;
    RETURN (Area(p,q,r)/FLOAT(np,LONG))*s;
  END OnTriangle;

PROCEDURE DotProductGrad(
    f, g: T; 
    tol: LONG; 
    depth: CARDINAL := LAST(CARDINAL);
  ): LONG = 
    
    VAR dotProd: LONG := 0.0d0;
    
    PROCEDURE Func(READONLY p: Point): LONG = 
      BEGIN 
        RETURN LR3.Dot(f.grad(p), g.grad(p))
      END Func;
    
    PROCEDURE Small(READONLY p, q, r: Point): BOOLEAN =
      BEGIN 
        WITH 
          dpq2 = LR3.DistSqr(p, q),
          dpr2 = LR3.DistSqr(p, r),
          dqr2 = LR3.DistSqr(q, r)
        DO 
          RETURN (dpq2 <= tol) AND (dpr2 <= tol) AND (dqr2 <= tol)
        END;
      END Small;
     
  BEGIN
    IF Disjoint(f.supp, g.supp) THEN RETURN 0.0d0 END;
    FOR i:= 0 TO LAST(f.data^) DO 
      WITH 
        nf = f.data[i].face,
        e1 = f.tri.side[nf],
        T10 = Dest(Lnext(e1)).c,
        T11 = Org(e1).c,
        T12 = Dest(e1).c,
        T1 = ARRAY[0..2] OF Point{T10, T11, T12}
      DO   
        FOR j:= 0 TO LAST(g.data^) DO 
          WITH  
            ng = g.data[j].face,
            e2 = g.tri.side[ng],
            T20 = Dest(Lnext(e2)).c,
            T21 = Org(e2).c,
            T22 = Dest(e2).c,
            T2 = ARRAY[0..2] OF Point{T20, T21, T22}
          DO 
            dotProd := dotProd + OnTwoTriangle(T1, T2, Func, Small, depth);
          END;
        END;
      END;  
    END;
    RETURN dotProd;
  END DotProductGrad;
    PROCEDURE DiscreteDotProduct(f, g: T; m: CARDINAL; F, G: Func := Id): LONG = 
  VAR fi, gi: CARDINAL;
      s: LONG := 0.0d0;
  BEGIN
    <* ASSERT f.tri = g.tri *>
    WITH
      tri = f.tri,
      c2b = tri.c2b^,
      b2c = tri.b2c^,
      fd = f.data^, fN = NUMBER(fd),
      gd = g.data^, gN = NUMBER(gd)
    DO
      IF Disjoint(f.supp, g.supp) AND F = Id AND G = Id THEN RETURN 0.0d0 END;
      SortTriangles(fd); fi := 0;
      SortTriangles(gd); gi := 0;
      FOR i := 0 TO LAST(tri.side^) DO
        VAR ft, gt: TriangleData;
        BEGIN
          IF fi < fN AND fd[fi].face = i THEN ft := fd[fi]; INC(fi) ELSE ft := NIL END;
          IF gi < gN AND gd[gi].face = i THEN gt := gd[gi]; INC(gi) ELSE gt := NIL END;

          IF (F = Id AND ft = NIL) OR (G = Id AND gt = NIL)
          THEN
            (* We can ignore this triangle, "F(f(p),p)*G(g(p),p) = 0" in it. *)
          ELSE
            s := s + DiscreteDotInTriangle(c2b[i], b2c[i], ft, F, gt, G, m);
          END;  
        END
      END;
    END;
    RETURN s
  END DiscreteDotProduct;




    PROCEDURE Small(READONLY p, q, r: Point): BOOLEAN =
      BEGIN 
        WITH 
          dpq2 = LR3.DistSqr(p, q),
          dpr2 = LR3.DistSqr(p, r),
          dqr2 = LR3.DistSqr(q, r)
        DO 
          RETURN (dpq2 <= tol) AND (dpr2 <= tol) AND (dqr2 <= tol)
        END;
      END Small;
    


PROCEDURE DiscreteDotProduct(f, g: T; m: CARDINAL; F, G: Func := Id): LONG;
  (* 
    Functions "f" and "g" must have the same triangulation.
    Computes the discrete dot product of "F(f)" and "G(g)", 
    defined on each triangle as "(A/N)*SUM(ff(p[k])*gg(p[k]))", 
    where "ff(p) = F(f(p),p)", "gg(p) = G(g(p),p)",
    "N = m*(m+1)/2", "A" is the area of the triangle, 
    and the "p[k]" are a triangular arrangement of 
    "N" points in that triangle. *)
    



  PROCEDURE ComputeVertexSupportingPlane(e: Arc): LR4.T =
    BEGIN
      WITH
        m = Degree(e),
        s = NEW(REF ARRAY OF Arc, m)^
      DO
        FOR i := 0 TO m-1 DO
          s[i] := e; e := Onext(e)
        END;
        RETURN ComputeSupportingPlane(s)
      END
    END ComputeVertexSupportingPlane;
    


  <*UNUSED*>
  PROCEDURE ShowTriangle(READONLY P, Q, R: Coords) =
    VAR pa, pb, qa, qb, ra, rb: LONGREAL;
    BEGIN
      Homog(P, pa, pb);
      Homog(Q, qa, qb);
      Homog(R, ra, rb);
      f.setFillColor(PSPlot.Gray(0.90));
      f.triangle(pa, pb, qa, qb, ra, rb);
      f.segment(pa, pb, qa, qb);
      f.segment(qa, qb, ra, rb);
      f.segment(ra, rb, pa, pb);
  END ShowTriangle;



PROCEDURE AddCoords(READONLY a, b: Coords): Coords =
  VAR r: Coords;
  BEGIN
    FOR i := 0 TO 2 DO r[i] := a[i] + b[i] END;
    RETURN r
  END AddCoords;
  
PROCEDURE ScaleCoords(READONLY a: Coords; s: LONGREAL): Coords =
  VAR r: Coords;
  BEGIN
    FOR i := 0 TO 2 DO r[i] := a[i] * s END;
    RETURN r
  END ScaleCoords;
                        


-- grau 5

      c102 = a0 * c202 + a1 * c112 + a2 * c103,
      d102d0 = c202 + a0*d202d0 + a1*d112d0 + a2*d103d0,
      d102d1 = a0*d202d1 + c112 + a1*d112d1 + a2*d103d1,
      d102d2 = a0*d202d2 + a1*d112d2 + c103 + a2*d103d2,
      
      c030 = a0 * c130 + a1 * c040 + a2 * c031,
      d030d0 = c130 + a0*d130d0 + a1*d040d0 + a2*d031d0,
      d030d1 = a0*d130d1 + c040 + a1*d040d1 + a2*d031d1,
      d030d2 = a0*d130d2 + a1*d040d2 + c031 + a2*d031d2,
      
      c021 = a0 * c121 + a1 * c031 + a2 * c022,
      d021d0 = c121 + a0*d121d0 + a1*d031d0 + a2*d022d0,
      d021d1 = a0*d121d1 + c031 + a1*d031d1 + a2*d022d1,
      d021d2 = a0*d121d2 + a1*d031d2 + c022 + a2*d022d2,
      
      c012 = a0 * c112 + a1 * c022 + a2 * c013,
      d012d0 = c112 + a0*d112d0 + a1*d022d0 + a2*d013d0,
      d012d1 = a0*d112d1 + c022 + a1*d022d1 + a2*d013d1,
      d012d2 = a0*d112d2 + a1*d022d2 + c013 + a2*d013d2, 
      
      c003 = a0 * c103 + a1 * c013 + a2 * c004,
      d003d0 = c103 + a0*d103d0 + a1*d013d0 + a2*d004d0,
      d003d1 = a0*d103d1 + c013 + a1*d013d1 + a2*d004d1,
      d003d2 = a0*d103d2 + a1*d013d2 + c004 + a2*d004d2,
      

--- grau 6
      c102 = a0 * c202 + a1 * c112 + a2 * c103,
      d102d0 = c202 + a0*d202d0 + a1*d112d0 + a2*d103d0,
      d102d1 = a0*d202d1 + c112 + a1*d112d1 + a2*d103d1,
      d102d2 = a0*d202d2 + a1*d112d2 + c103 + a2*d103d2,
      
      c030 = a0 * c130 + a1 * c040 + a2 * c031,
      d030d0 = c130 + a0*d130d0 + a1*d040d0 + a2*d031d0,
      d030d1 = a0*d130d1 + c040 + a1*d040d1 + a2*d031d1,
      d030d2 = a0*d130d2 + a1*d040d2 + c031 + a2*d031d2,
      
      c021 = a0 * c121 + a1 * c031 + a2 * c022,
      d021d0 = c121 + a0*d121d0 + a1*d031d0 + a2*d022d0,
      d021d1 = a0*d121d1 + c031 + a1*d031d1 + a2*d022d1,
      d021d2 = a0*d121d2 + a1*d031d2 + c022 + a2*d022d2,
      
      c012 = a0 * c112 + a1 * c022 + a2 * c013,
      d012d0 = c112 + a0*d112d0 + a1*d022d0 + a2*d013d0,
      d012d1 = a0*d112d1 + c022 + a1*d022d1 + a2*d013d1,
      d012d2 = a0*d112d2 + a1*d022d2 + c013 + a2*d013d2, 

      c003 = a0 * c103 + a1 * c013 + a2 * c004,
      d003d0 = c103 + a0*d103d0 + a1*d013d0 + a2*d004d0,
      d003d1 = a0*d103d1 + c013 + a1*d013d1 + a2*d004d1,
      d003d2 = a0*d103d2 + a1*d013d2 + c004 + a2*d004d2,
      
---
      


(* From SPPWFunction.m3 98-01-09 *)

PROCEDURE TriLocate(READONLY p: Point; e: Arc): Arc;
  (*
    Given a point "p" on the sphere, and an arc "e" 
    of a spherical triangualtion, returns an arc
    "a" such that Left9a) is the triangle that 
    contains "p". Uses a ``walking'' algorithm. *)
 
PROCEDURE TriLocate(READONLY p: Point; e: Arc): Arc =
  VAR count: CARDINAL := 0;
  BEGIN
    WITH q = LR3.T{0.01676565151d0, 0.02173421513d0, 0.03145257313d0} DO 
    LOOP 
      WITH 
        a = Org(e).c,
        b = Dest(e).c,
        c = Dest(Onext(e)).c
      DO
        IF count > 100 THEN 
          Wr.PutText(stderr, "SPPWFunction.TriLocate: Loop?\n");
          Wr.PutText(stderr, "a = "); PrintPt(a); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "b = "); PrintPt(b); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "c = "); PrintPt(c); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "\n");
          IF count > 110 THEN Process.Exit(1) END;
        END;
        IF SamePoint(p, a) THEN RETURN e
        ELSIF SamePoint(p, b) THEN RETURN e
        ELSIF SamePoint(p, c) THEN RETURN e
        ELSIF PositiveSide(q, b, a, p) THEN e := Sym(e)
        ELSIF NOT PositiveSide(q, c, a, p) THEN e := Onext(e)
        ELSIF NOT PositiveSide(q, b, c, p) THEN e := Dprev(e) 
        ELSE RETURN e 
        END
      END;
      INC(count)
    END
   END
  END TriLocate;

(* *)

    PROCEDURE Det(u, v, w: Coords): LONGREAL =
      BEGIN
        WITH
          u0 = FLOAT(u[0], LONGREAL),
          u1 = FLOAT(u[1], LONGREAL),
          u2 = FLOAT(u[2], LONGREAL),
          v0 = FLOAT(v[0], LONGREAL),
          v1 = FLOAT(v[1], LONGREAL),
          v2 = FLOAT(v[2], LONGREAL),
          w0 = FLOAT(w[0], LONGREAL),
          w1 = FLOAT(w[1], LONGREAL),
          w2 = FLOAT(w[2], LONGREAL),
          det = 
            + u0*(v1*w2 - v2*w1) 
            - u1*(v0*w2 - v2*w0) 
            + u2*(v0*w1 - v1*w0)
        DO
          (* Wr.PutText(stderr, "Det = " & Fmt.LongReal(det) & "\n"); *)
          RETURN det
       END;
    END Det;
  
  
    PROCEDURE Alt(
        p,q,r:Coords;
        d:LONGREAL;
      ): LONGREAL =

      VAR p1,p2,p3: Coords;
          a,b,c:LONGREAL;
      BEGIN 

        p1[0]:= p[1];
        p1[1]:= p[2];
        p1[2]:= 1.0d0;

        p2[0]:= q[1];
        p2[1]:= q[2];
        p2[2]:= 1.0d0;

        p3[0]:= r[1];
        p3[1]:= r[2];
        p3[2]:= 1.0d0;

        a:= ABS(LR3Extras.Det(p1,p2,p3));

        p1[0]:= p[0];
        p1[1]:= p[2];
        p1[2]:= 1.0d0;

        p2[0]:= q[0];
        p2[1]:= q[2];
        p2[2]:= 1.0d0;

        p3[0]:= r[0];
        p3[1]:= r[2];
        p3[2]:= 1.0d0;


        b:= ABS(Det(p1,p2,p3));

        p1[0]:= p[0];
        p1[1]:= p[1];
        p1[2]:= 1.0d0;

        p2[0]:= q[0];
        p2[1]:= q[1];
        p2[2]:= 1.0d0;

        p3[0]:= r[0];
        p3[1]:= r[1];
        p3[2]:= 1.0d0;

        c:= ABS(LR3Extras.Det(p1,p2,p3));
        WITH
          alt = d/Math.sqrt(a*a + b*b + c*c)
        DO
          (* Wr.PutText(stderr, "Alt = " & Fmt.LongReal(alt) & "\n"); *)
          RETURN alt
        END;
      END Alt;  

    VAR d, volume,altura,area : LONGREAL;
      d := ABS(Det(p,q,r));
      volume := d/6.0d0;
      altura := Alt(p,q,r,d);
      area := (3.0d0*volume)/altura;
    
  PROCEDURE NormalizeSiteCoords(VAR c: Coords) =
    BEGIN
      WITH 
        x = c[0],
        y = c[1],
        z = c[2],
        xx = FLOAT(x, LONGREAL),
        yy = FLOAT(y, LONGREAL),
        zz = FLOAT(z, LONGREAL),
        dd = Math.sqrt(xx*xx + yy*yy + zz*zz)
      DO
        x := FLOAT(xx/dd, LONGREAL);
        y := FLOAT(yy/dd, LONGREAL);
        z := FLOAT(zz/dd, LONGREAL)
      END
    END NormalizeSiteCoords;
 
  PROCEDURE  Integral(READONLY p, q, r: SPTriang.Coords): LONGREAL =
    BEGIN 
      WITH
        area = Area(p, q, r),
        intg = ((func(p) + func(q) + func(r))*area)/3.0d0
      DO
        (* Wr.PutText(stderr, "Integral = " & Fmt.LongReal(intg) & "\n"); *)
        RETURN intg
      END
    END Integral;


          Math.sqrt((p[0] - q[0])*(p[0] - q[0]) + (p[1] - q[1])*(p[1] - q[1]) + 
            (p[2] - q[2])*(p[2] - q[2])),
          dpr = Math.sqrt((p[0] - r[0])*(p[0] - r[0]) + (p[1] - r[1])*(p[1] - r[1]) + 
            (p[2] - r[2])*(p[2] - r[2])),   
          dqr = Math.sqrt((q[0] - r[0])*(q[0] - r[0]) + (q[1] - r[1])*(q[1] - r[1]) + 
            (q[2] - r[2])*(q[2] - r[2]))  



  PROCEDURE B(
      e: Arc;
      VAR b: Basis;
      tri: REF SPTriang.T;
      d: CARDINAL;
      VAR nb: CARDINAL; 
    ) = 
    
      VAR a: Arc;
        bar: ARRAY[0..2] OF LONGREAL;
        vc, w : ARRAY[0..9] OF LONGREAL;
        
    BEGIN
      WITH 
        elem = NEW(T),
        order = SPTriang.Degree(e)
      DO 
        b[nb] := elem;
        elem.tri := tri;
        elem.supp := ComputeVertexSupportingPlane(e);
        elem.data:= MakeTriangles(order, ncv);
        FOR i := 0 TO 9 DO vc[i] := 0.0d0 END;
        WITH 
            t0 = NARROW(elem.data[0], TriangleData)
          DO 
            t0.face := Left(e).num;
            t0.degree := d;
            ComputeVC7AndVC8(e, vc);
            SetCoeffsVertex(t0.coeffs^, vc, e);
        END;
      
        a := Onext(e); 
        FOR i := 0 TO 9 DO vc[i] := 0.0d0 END;
        vc[5] := 1.0d0;
        
        FOR k:= 1 TO order-3 DO 
          WITH 
            t0 = NARROW(elem.data[k], TriangleData)
          DO 
            t0.face := Left(a).num;
            t0.degree := d;
            ComputeVC7AndVC8(a, vc);
            SetCoeffsVertex(t0.coeffs^, vc, a);
          
          END;
          bar := ComputeEdgeContinuityCoeffs(Sym(Onext(e)), tri);
          vc[4] := vc[4]*bar[0] + bar[1]*vc[5] + bar[2]*vc[2];
          vc[2] := vc[1]*bar[0] + bar[1]*vc[2] + bar[2]*vc[0];
          vc[1] := vc[2];
          vc[3] := vc[5];
          vc[5] := 0.0d0;
          a := Onext(a);
        END;
         
        vc[5] := ( - bar[0]*vc[4] - bar[2]*vc[2])/bar[1];
        WITH 
          t0 = NARROW(elem.data[order-2], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.degree := d;
          ComputeVC7AndVC8(a, vc);
          SetCoeffsVertex(t0.coeffs^, vc, a);
        END;
        
        a := Onext(a);
        FOR i := 0 TO 9 DO w[i] := 0.0d0 END;
        w[3] := vc[5];
        w[1] := vc[2];
        
       WITH 
          t0 = NARROW(elem.data[order-1], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.degree := d;
          ComputeVC7AndVC8(a, w);
          SetCoeffsVertex(t0.coeffs^, w, a);
        END;
      END;
      Print(b[nb]);
      INC(nb);  
    END B;
    
   PROCEDURE C(
      e:Arc;
      VAR b:Basis;
      tri: REF SPTriang.T;
      d:CARDINAL;
      VAR nb:CARDINAL; 
    )=
      
      VAR a: Arc;
        bar: ARRAY[0..2] OF LONGREAL;
        vc, w :  ARRAY[0..9] OF LONGREAL;
            
    BEGIN
      WITH 
        elem = NEW(T),
        order = SPTriang.Degree(e)
      DO 
        b[nb] := elem;
        elem.tri := tri;
        elem.supp := ComputeVertexSupportingPlane(e);
        elem.data:= MakeTriangles(order, ncv);
        a := e;
        FOR k := 0 TO 1 DO 
          FOR i := 0 TO 9 DO vc[i] := 0.0d0 END;
          WITH 
            t0 = NARROW(elem.data[k], TriangleData)
          DO 
            t0.face := Left(a).num;
            t0.degree := d;
            ComputeVC7AndVC8(a, vc);
            SetCoeffsVertex(t0.coeffs^, vc, a);
          END;
          a := Onext(a); 
        END;
      
        FOR i := 0 TO 9 DO vc[i] := 0.0d0 END;
        vc[5] := 1.0d0;
        
        WITH 
          t0 = NARROW(elem.data[order-3], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.degree := d;
          ComputeVC7AndVC8(a, vc);
          SetCoeffsVertex(t0.coeffs^, vc, a);
        END;
        
        bar := ComputeEdgeContinuityCoeffs(Sym(Onext(a)), tri);  
        vc[4] := vc[4]*bar[0] + bar[1]*vc[5] + bar[2]*vc[2];
        vc[2] := vc[1]*bar[0] + bar[1]*vc[2] + bar[2]*vc[0];
        vc[1] := vc[2];
        vc[3] := vc[5];
        vc[5] := (- bar[0]*vc[4] - bar[2]*vc[2])/bar[1];
        a:= Onext(a);
        WITH 
          t0 = NARROW(elem.data[order-2], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.degree := d;
          ComputeVC7AndVC8(a, vc);
          SetCoeffsVertex(t0.coeffs^, vc, a);
        END;
       
        FOR i := 0 TO 9 DO w[i] := 0.0d0 END;
        w[3] := vc[5];
        w[1] := vc[2];
        a:= Onext(a);     
      
        WITH 
          t0 = NARROW(elem.data[order-1], TriangleData)
        DO 
          t0.face := Left(a).num;
          t0.degree := d;
          ComputeVC7AndVC8(a, w);
          SetCoeffsVertex(t0.coeffs^, w, a);
        END;
      END;
      Print(b[nb]);
      INC(nb);  
    END C;
    
  


        PROCEDURE Changeijk(VAR ijk: SPHomoLabel.T; i,j,k : CARDINAL)=
          VAR aux: SPHomoLabel.T;
            
          BEGIN
            aux[0]:= i;
            aux[1]:= j;
            aux[2]:= k;
            ijk:= aux
          END Changeijk;
      

      PROCEDURE Swap(VAR s, t : LONGREAL)=
        VAR aux: LONGREAL;
        BEGIN
          aux:= s;
          s:= t;
          t:= s
        END Swap;


PROCEDURE CollectOut(e: Arc; nSites: CARDINAL): REF ARRAY OF Arc =
  BEGIN
    WITH ro = NEW(REF ARRAY OF Arc, nSites) DO
      CollectVertices(e, ro^);
      RETURN ro
    END
  END CollectOut;



          IF d = 5 THEN 
          ELSE 
            IF a = e THEN
              FOR i:= 0 TO 5 DO
                c[i]:= vc[i]
              END;
            ELSIF Sym(Onext(e)) = a  THEN 
              FOR i:= 0 TO 5 DO 
                ijk := SPHomoLabel.IndexToLabel(i,6);
                ijk := SPHomoLabel.T{ijk[1], ijk[2], ijk[0]};
                indr:= SPHomoLabel.LabelToIndex(ijk,6);
                c[indr]:= vc[i]
              END;
            ELSIF  a = Lnext(e) THEN 
              FOR i:= 0 TO 5 DO 
                ijk:= SPHomoLabel.IndexToLabel(i,6);
                ijk := SPHomoLabel.T{ijk[2], ijk[0], ijk[1]};
                indr:= SPHomoLabel.LabelToIndex(ijk,6);
                c[indr]:= vc[i]
              END;
            END;  
          END;