INTERFACE SPPWFunction;

(* Piecewise functions on the sphere *)

IMPORT LR4, LR3, Wr, Rd, SPTriang, SPFunction, SPMatrix;
FROM SPTriang IMPORT TriangleNumber, Arc;

TYPE
  LONG = LONGREAL;

TYPE
  T <: PublicT;

  PublicT = SPFunction.T OBJECT
      tri: REF SPTriang.T;       (* The domain partition. *)
      supp: LR4.T;               (* A plane such that "f" is zero on the neg. side *)
      data: REF Triangles;       (* Data for triangles where "f" is non-NULL. *) 
    METHODS
      write(wr: Wr.T); 
        (* 
          Writes the function (but not the triangulation) to "wr". *)
          
      read(rd: Rd.T);  
        (*
          Reads the function from "rd", in the format
          created by "write". Assumes the triangulation "tri"
          is the correct one. *)
      
      locate(READONLY p: Point): TriangleNumber;
        (*
          Finds the triangle of "f.tri" that contains "p".
          Most efficient if "p" is in the same triangle as
          the previous call. *)
    END;

  Point = SPFunction.Point;       (* Cartesian coordinates of point *)
  Barycentric = LR3.T;            (* Triangle-relative coordinates of point *)
  Gradient = SPFunction.Gradient;
  
  Triangles = ARRAY OF TriangleData;

  TriangleData = OBJECT
      face: TriangleNumber; (* Number of triangle *)
    METHODS 
      eval(READONLY p: Point; READONLY b: Barycentric): LONGREAL;
      grad(READONLY p: Point; READONLY b: Barycentric): Gradient;
      type(): CHAR;          (* The actual type of TriangleData *)
      write(wr: Wr.T);       (* Writes the triangle's function to "wr" *)
      read(rd: Rd.T);        (* Reads the triangle's function from "rd" *)
      toMaple(wr: Wr.T);     (* Writes the triangle's function in Maple format *)

      copy(): TriangleData;  
        (* 
          Makes a clone of this TriangleData, and initializes it 
          with zero function. *)

      add(a: LONG; h: TriangleData);  
        (* 
          Adds the function of triangle "h", times the factor "a",
          to the function of this triangle (which musthave the
          same actual type as "h"). *)
    END;
    
  Arcs = ARRAY OF Arc;
  
  Basis = ARRAY OF T; 
    (* 
      A basis for a space of piecewise spherical functions. *)
      
  BasisCoeffs = ARRAY OF LONGREAL;
  
PROCEDURE WriteBasis(wr: Wr.T; READONLY b: Basis);
  (*
    Writes a basis "b" to "wr", in a format that can be 
    read back by "Read".  Does NOT write the triangulation. *)
   
PROCEDURE ReadBasis(rd: Rd.T; tri: REF SPTriang.T): REF Basis;
  (*
    Reads from "rd" a basis for a space of piecewise
    spherical functions defined on the triangulation "tri".
    Assumes the contents of "rd" was created by "write" above,
    and that the triangulation was exactly the same 
    (including vertex, face,and edge numbers). *)
   
PROCEDURE ComputeSupportingPlane(
    READONLY d: Triangles; (* A list of triangles. *)
    READONLY side: Arcs;   (* "side" table of trangulation. *)
  ): LR4.T;
  (*
    Returns a plane "P" such that all spherical triangles "d[i]"
    are contained in the positive half-space of "P".
    The second parameter should be the "side" table
    of the triangulation. *)
    
PROCEDURE DotProduct(
    f, g: T; 
    order: CARDINAL;
    F, G: Func;
    tri: REF SPTriang.T;
  ): LONG;
  (* 
    Computes the dot product of "F(f)" and "G(g)", 
    that is the integral of "F(f(p),p)*G(g(p),p)" over the whole sphere. 
    The dot product is approximated by a weigted summation over a finite 
    set of sample points.
    
    If "tri # NIL", the procedure uses "M = (order+1)(order+2)/2"
    sample points in each triangle of "tri".  If "tri = NIL", 
    the procedure requires that "f.tri = g.tri", and uses that
    triangulation. *)
    
PROCEDURE DotProductGrad(
    f, g: T; 
    order: CARDINAL;
    tri: REF SPTriang.T;
  ): LONG;
  (*
    Computes the dot product of "Grad(f)" and "Grad(g)", that is the
    integral of "LR3.Dot(Grad(f),Grad(g))" over the whole sphere.
    The integral is approximated as described for "DotProduct". *)
    
TYPE Func = PROCEDURE (u: LONG; READONLY p: Point): LONG;
    
PROCEDURE Id(u: LONG; READONLY p: Point): LONG;
  (*
    The identity function: returns "u", ignores "p". *)
    
PROCEDURE LinComb(
    READONLY a: ARRAY OF LONG; 
    READONLY f: ARRAY OF T;
  ): T;
  (*
    Returns the linear combination "SUM(a[i]*f[i])",
    a new spherical function. All functions "f[i]"
    must have the same triangulation, which will 
    also be used for the result. *)
    
PROCEDURE BuildMatrix(
    READONLY F, G: Basis; 
    order: CARDINAL;
    tri: REF SPTriang.T;
  ): SPMatrix.T; 
  (*
    Builds the rigidity matrix "m" for bases "F" and "G",
    i.e. "m[i,j] = Dot(F[i], G[j])".  The parameters
    "order" and "tri" define the dot product
    used; see the documentation of "DotProd". *)  
    
PROCEDURE BuildMatrixGrad(
    READONLY F, G: Basis; 
    order: CARDINAL;
    tri: REF SPTriang.T;
  ): SPMatrix.T;
  (*
    Builds the rigidity matrix "m" for bases "F" and "G", i.e. "m[i,j]
    = Dot(Grad(F[i]), Grad(G[j]))".  The parameters "order" and "tri"
    define the dot product used; see the documentation of
    "DotProdGrad". *)
    
PROCEDURE ChangeCoordsBasis( 
    READONLY a: SPMatrix.Vector;
    READONLY FG: SPMatrix.T;
    READONLY GG: SPMatrix.T;
    VAR b : SPMatrix.Vector;
  );  
  (* 
    Given the coefficients "a[0..m-1]" of a function "f(p)" in terms of
    a functional basis "F[0..m-1]", computes the coefficients
    "b[0..n-1]" of the least squares approximation "g(p)" to "f(p)" in
    some other basis "G[0..n-1]".  Requires the rigidity matrices
    "FG[i,j] = Dot(F[i].G[j]) and "GG[i,j] = Dot(G[i],G[j])". *)
       
END SPPWFunction.