MODULE SPGGMatrix EXPORTS Main;

(* Computes the scalar product matrices "G" and "S" for a PW spherical basis *)

IMPORT Wr, Fmt, Thread, ParseParams, Process;
IMPORT SPTriang, SPMatrix, SPPWFunction; 
FROM SPFile IMPORT OpenRead, OpenWrite;

FROM Math IMPORT sqrt;

FROM Stdio IMPORT stderr;

TYPE 
  Options = RECORD
      triName: TEXT;       (* Triangulation file for basis. *)
      baseName: TEXT;      (* Function file of basis. *)
      outName: TEXT;       (* Prefix for matrix files. *)
      order: CARDINAL;     (* Order of sampling grid in each triangle. *)
      writeG: BOOL;        (* TRUE to compute and write the "G" matrix. *)
      writeS: BOOL;        (* TRUE to compute and write the "S" matrix. *)
    END;
  
TYPE 
  LONG = LONGREAL;
  BOOL = BOOLEAN;

<* FATAL Wr.Failure, Thread.Alerted *>

PROCEDURE Main() =
  VAR G, S: SPMatrix.T;
  BEGIN
    WITH 
      o = GetOptions(),
      rdtri = OpenRead(o.triName & ".tri"),
      tri = MakeTriRef(SPTriang.Read(rdtri)),
      
      rdbas = OpenRead(o.baseName & ".bas"),
      basis = SPPWFunction.ReadBasis(rdbas, tri)
    DO 
      IF o.writeG THEN
        G := SPPWFunction.BuildMatrix(basis^, basis^, o.order, tri);
        WITH wr = OpenWrite(o.outName & "-G.mat") DO
          SPMatrix.Write(wr, G)
        END;
        CheckAngles(G);
      END;
      IF o.writeS THEN
        S := SPPWFunction.BuildMatrixGrad(basis^, basis^, o.order, tri);
        WITH wr = OpenWrite(o.outName & "-S.mat") DO
          SPMatrix.Write(wr, S)
        END;
      END;
    END;  
  END Main;

PROCEDURE MakeTriRef(READONLY tri: SPTriang.T): REF SPTriang.T = 
  (*
    Droga de Modula-3... *)
  BEGIN
    WITH t = NEW(REF SPTriang.T) DO t^ := tri; RETURN t END
  END MakeTriRef;
  
PROCEDURE CheckAngles(M: SPMatrix.T) =
  BEGIN
    <* ASSERT M.rows = M.cols *>
    WITH
      N = M.rows,
      elem = M.r^,
      d = NEW(REF ARRAY OF LONG, N)^
    DO
      FOR k := 0 TO LAST(elem) DO
        WITH e = elem[k] DO 
          IF e.row = e.col THEN d[e.row] := sqrt(e.va) END 
        END
      END;
      FOR k := 0 TO LAST(elem) DO
        WITH 
          e = elem[k],
          c = e.va/d[e.row]/d[e.col]
        DO 
          IF (ABS(c) > 0.9999d0) AND (e.row # e.col) THEN 
            Wr.PutText(stderr, "basis elements " & FI(e.row,4) & " and " & FI(e.col,4));
            Wr.PutText(stderr, " cosine = " & FLR(c, 15,9) & "\n");
          END
        END
      END;
    END
  END CheckAngles;

PROCEDURE FLR(x: LONGREAL; wid, prec: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, style:= Fmt.Style.Fix, prec:= prec), wid)
  END FLR;

PROCEDURE FI(x: CARDINAL; wid: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Int(x), wid)
  END FI;

PROCEDURE GetOptions (): Options =
  VAR o: Options;
  BEGIN
    TRY
      WITH  pp = NEW(ParseParams.T).init(stderr) DO 
          
        pp.getKeyword("-triName");                               
        o.triName := pp.getNext();  

        pp.getKeyword("-baseName");                               
        o.baseName := pp.getNext();  

        pp.getKeyword("-outName");                               
        o.outName := pp.getNext();  
                  
        pp.getKeyword("-order");
        o.order := pp.getNextInt(0,LAST(CARDINAL));
        
        o.writeG := pp.keywordPresent("-writeG");
        o.writeS := pp.keywordPresent("-writeS");
      
        pp.finish();
      END
    EXCEPT                                                            
    | ParseParams.Error =>                                              
        Wr.PutText(stderr, "Usage: SPMatrix \\\n");
        Wr.PutText(stderr, "  -triName <triName> \\\n");
        Wr.PutText(stderr, "  -baseName <baseName> \\\n");
        Wr.PutText(stderr, "  -outName <name> \\\n");
        Wr.PutText(stderr, "  -order <num> \\\n");
        Wr.PutText(stderr, "  [ -writeG ] [ -writeS ]\n");
        Process.Exit (1);
    END;
    RETURN o
  END GetOptions;
  
BEGIN
  Main();
END SPGGMatrix.