MODULE PZTestGaussian EXPORTS Main;

(* Basic test of the parametric filter routine. *)
(* Last edited on 1999-08-06 21:26:23 by hcgl *)

IMPORT  ParseParams, Process, Wr;
IMPORT  Thread, Stdio, Fmt;
IMPORT  LR3;
IMPORT PZTypes, PZSmooth, PZLR3Chain, PZLRChain;

FROM Math IMPORT exp, sqrt;
FROM Stdio IMPORT stderr, stdout;
FROM PZTypes IMPORT LONG, LONGS, NAT, INT, BOOL;

(*
  Builds a Dirac-like curve with the pulse centered at t=0.
  Convolves it with a Gaussian of specified duration. 
  Writes to standard output a table with 7 columns:
    time
    computed convolution samples:
      value u
      derivative
      second derivative
    theoretically expected values:
      value 
      derivative
      second derivative
*)
 
<* FATAL Wr.Failure, Thread.Alerted *>

TYPE
  LONG = LONGREAL;
  
  MismatchOptions = PZMismatch.Options;
  Options = RECORD
      sigma: LONG;  (* The characteristic duration of the gaussian. *)
    END;
    
CONST
  Pi = 3.14159265358979d0;
    
PROCEDURE Main() =
  BEGIN
    WITH
      o = GetOptions(),
      (* A dirac pulse centered at t=0: *)
      t = PZLRChain.T{ 
        -10.0000d0, 
        -00.0001d0, 
        000.0000d0,
        +00.0001d0
      },
      tPeriod = 20.0d0,
      tStart = t[0],
      p = PZLR3Chain.T{
        LR3.T{     0.0d0, 0.0d0, 0.0d0},
        LR3.T{     0.0d0, 0.0d0, 0.0d0},
        LR3.T{ 10000.0d0, 0.0d0, 0.0d0},
        LR3.T{     0.0d0, 0.0d0, 0.0d0}
      },
      (* Filtered curve: *)
      NSamples = 200,
      u = NEW(REF PZLR3Chain.T, NSamples)^,
      v = NEW(REF PZLR3Chain.T, NSamples)^,
      w = NEW(REF PZLR3Chain.T, NSamples)^
    DO 
      PZSmooth.Polygonal(o.sigma, p, t, tPeriod, tStart, u, v, w);
      FOR i := 0 TO LAST(u) DO
        WITH
          sigma2 = o.sigma*o.sigma,
          t = tStart + tPeriod * FLOAT(i,LONG)/FLOAT(NSamples, LONG),
          t2 = t*t,
          ui = u[i][0],
          vi = v[i][0],
          wi = w[i][0],
          A = 1.0d0/(o.sigma * sqrt(Pi)),
          ux = + A * exp(-t2/sigma2),
          vx = - ux * 2.0d0 * t / sigma2,
          wx = + ux * (4.0d0 * t2/sigma2 - 2.0d0) / sigma2
        DO
          Wr.PutText(stdout, " ");
          Wr.PutText(stdout, FLR(t, 11, 6));
          Wr.PutText(stdout, " ");
          Wr.PutText(stdout, FLR(ui, 11, 6));
          Wr.PutText(stdout, " ");
          Wr.PutText(stdout, FLR(vi, 11, 6));
          Wr.PutText(stdout, " ");
          Wr.PutText(stdout, FLR(wi, 11, 6));
          Wr.PutText(stdout, " ");
          Wr.PutText(stdout, FLR(ux, 11, 6));
          Wr.PutText(stdout, " ");
          Wr.PutText(stdout, FLR(vx, 11, 6));
          Wr.PutText(stdout, " ");
          Wr.PutText(stdout, FLR(wx, 11, 6));
          Wr.PutText(stdout, "\n");
        END
      END;
      Wr.Close(stdout)
    END
  END Main;
  
<*UNUSED*>
PROCEDURE Debug(msg: TEXT; READONLY x: LONGS) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(stderr, "\n");
    Wr.PutText(stderr, msg);
    FOR i := 0 TO LAST(x) DO
      IF i > 0 AND i MOD 5 = 0 THEN 
        Wr.PutText(stderr, "\n")
      ELSE
        Wr.PutText(stderr, " ")
      END;
      Wr.PutText(stderr, 
        Fmt.Pad(Fmt.LongReal(x[i], prec := 8, style := Fmt.Style.Fix), 12)
      )
    END;
    Wr.PutText(stderr, "\n");
  END Debug;
  
PROCEDURE FLR(x: LONG; w: CARDINAL; prec: CARDINAL): TEXT = 
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, prec := prec, style := Fmt.Style.Fix), w)
  END FLR;
  
PROCEDURE GetOptions(): Options =
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY

        pp.getKeyword("-sigma");
        o.sigma := pp.getNextLongReal();
 
        pp.finish();                                       
      
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PZTestGaussian \\\n");
          Wr.PutText(stderr, "  -sigma \n");
        Process.Exit(1);
      END;
    END;
    RETURN o
    
  END GetOptions;


BEGIN
  Main()
END PZTestGaussian.