MODULE SPIntTest EXPORTS Main;

(* Tests the integration of arbitrary funtions on the sphere *)

IMPORT Wr, Math, Fmt, Thread, ParseParams;
IMPORT LR3, Process, SPTriang;

FROM SPIntegral IMPORT OnTriangle; 
FROM Stdio IMPORT stderr;

TYPE 
  Options = RECORD
      order: CARDINAL;     (* Order of sampling grid in each triangle *)
    END;

TYPE 
  LONG = LONGREAL;
  Point = SPTriang.Coords;
  
<* FATAL Wr.Failure, Thread.Alerted *>

PROCEDURE Main() =
  VAR 
    s: LONGREAL := 0.0d0;
  BEGIN
    WITH 
      o = GetOptions()
    DO  

      PROCEDURE Func( READONLY p: Point):LONGREAL =
        BEGIN 
          (* RETURN 1.0d0 *)
          (* RETURN ABS(p[0]) *)
          (* RETURN p[0]*p[0] *)
             RETURN p[0]*p[0]*p[0]*p[0] 
          (* WITH u = (p[0] + p[1] + p[2]) 
           DO 
             RETURN (u*u*u*u)/9.0d0
           END; *)
          (*  WITH u = (p[0] + p[1]) 
           DO 
             RETURN (u*u*u*u)/4.0d0
           END; *) 
          (* RETURN p[0]*p[0] + p[1]*p[1] + p[2]*p[2] *) 
          (* WITH u = (p[0] + p[1] + p[2]) 
           DO 
             RETURN (u*u)/3.0d0
           END; *)
          (* RETURN ((p[0] + 1.0d0)/2.0d0)^8 *)
        END Func;

      BEGIN
        FOR x := -1 TO +1 BY 2 DO
          FOR y := -1 TO +1 BY 2 DO
            FOR z := -1 TO +1 BY 2 DO
              WITH
                p = LR3.T{FLOAT(x, LONG), 0.0d0, 0.0d0},
                q = LR3.T{0.0d0, FLOAT(y, LONG), 0.0d0},
                r = LR3.T{0.0d0, 0.0d0, FLOAT(z, LONG)}
              DO
                s := s + OnTriangle(p, q, r, Func, o.order);
                Wr.PutText(stderr, "Partial integral = ");
                Wr.PutText(stderr, FLR(s, 15, 8));
                Wr.PutText(stderr, "\n");
              END
            END
          END
        END
      END
    END;  
  END Main;
  
PROCEDURE FLR(x: LONGREAL; wid, prec: CARDINAL): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, style:= Fmt.Style.Sci, prec:= prec), wid)
  END FLR;

PROCEDURE GetOptions (): Options =
  VAR o: Options;
  (*  obs, upp: HLR3.Point;*)
  BEGIN
    TRY
      WITH  pp = NEW(ParseParams.T).init(stderr) DO 
          
        pp.getKeyword("-order");
        o.order := pp.getNextInt(0,LAST(CARDINAL));
        
        pp.finish();
      END
    EXCEPT                                                            
    | ParseParams.Error =>                                              
        Wr.PutText(stderr, "Usage: SPIntTest \\\n");
        Wr.PutText(stderr, "  -order <num>\n");
        Process.Exit (1);
    END;
    RETURN o
  END GetOptions;
  
BEGIN
  Main();
END SPIntTest.