MODULE BITSpectrum EXPORTS Main;

(* Writes the Fourier power spectrum of a uniformly sampled curve. *)

(* Last edited on 1998-12-09 03:15:14 by hcgl *)

IMPORT  ParseParams, Process, Wr, FileWr, Rd, Lex;
IMPORT  OSError, Thread, Stdio, Fmt, FPut, FileFmt, FGet;
FROM Math IMPORT sqrt;
FROM Stdio IMPORT stderr;
FROM PZTypes IMPORT LONG, LONGS, NAT, INT, BOOL;
FROM Stdio IMPORT stdin;
FROM PZTypes IMPORT LONG, LONGS, NAT, INT, BOOL;

<* FATAL Wr.Failure, Thread.Alerted, OSError.E,  Rd.Failure  *>

TYPE
  LONG = LONGREAL;
  
  FreqData = RECORD
    freq: LONG;
    lambda: LONG;
    rms: LONG
  END;
  
  FreqPower = ARRAY OF FreqData;
  
  MismatchOptions = PZMismatch.Options;
  Options = RECORD
      outName: TEXT;     (* Output file name (without extension).  *)
      nFreqs: CARDINAL;  (* Number of frequences in output file     *)   
  END;
  (*
  
  
  *)
  
PROCEDURE Main() =
  VAR 
     k : CARDINAL;
  BEGIN
    WITH
      o = GetOptions(),
      nFreqs = o.nFreqs,
      nSpectra = NEW(REF  LONGS, nFreqs)^,
      accPower = NEW(REF LONGS, nFreqs)^,
      freq    = NEW(REF LONGS, nFreqs)^,
      lambda    = NEW(REF LONGS, nFreqs)^,
      accMean = NEW(REF LONGS,nFreqs)^
    DO 
      FOR k := 0 TO nFreqs-1 DO 
        nSpectra[k] := 0.0d0; accPower[k] := 0.0d0;
        freq[k] := 0.0d0; lambda[k] := 0.0d0 
      END;
      LOOP
        Lex.Skip(stdin);
        IF Rd.EOF(stdin) THEN EXIT END;
        ReadAndAccumulateFreqSpectrum
          (stdin, nSpectra, accPower, freq, lambda, nFreqs)
      END;
      k:= 0; 
      WHILE   (k < nFreqs) DO 
        IF nSpectra[k] > 0.0d0 THEN 
           accMean[k] := accPower[k]/nSpectra[k];
           k := k+1
        END   
      END;
      WriteMeanFreqPower(o.outName, accMean, freq, lambda);
    END 
  END Main;

 
PROCEDURE ReadAndAccumulateFreqSpectrum(
    rd: Rd.T;
    VAR nSpectra: LONGS;
    VAR accPower: LONGS;
    VAR freq: LONGS;
    VAR lambda: LONGS;
    nFreqs: CARDINAL
  ) =
  BEGIN  
    WITH 
      s = ReadFreqPower(rd, nFreqs)^,
      n = NUMBER(s)  
    DO  
     FOR k := 0 TO n-1 DO
       nSpectra[k] := nSpectra[k] + 1.0d0;
       accPower[k]:= accPower[k] + s[k].rms*s[k].rms;
       lambda[k] := s[k].lambda;
       freq[k] := s[k].freq;
     END;  
    END; 
  END ReadAndAccumulateFreqSpectrum;


PROCEDURE  ReadFreqPower(rd: Rd.T; nFreqs: CARDINAL): REF FreqPower = 
  VAR k: CARDINAL;
  BEGIN  
    WITH 
      (* rd = FileRd.Open(fileName), *)
      su = NEW(REF FreqPower, nFreqs),
      s = su^
    DO 
      EVAL FileFmt.ReadComment(rd, '#');
      (* Wr.PutText(stderr, "reading file  ...\n");  *)
      (* Wr.PutText(stderr, cmt); Wr.PutText(stderr, "\n");  *)
      k:= 0;
      FGet.Skip(rd);
      FOR k := 0 TO nFreqs-1 DO
        s[k].freq :=   FGet.LongReal(rd); 
        (* Wr.PutText(stderr, Fmt.LongReal(s[k].freq) & "  " & "\n");  *)
        s[k].lambda :=   FGet.LongReal(rd);
        (* Wr.PutText(stderr, Fmt.LongReal(s[k].lambda) & "  "); *)
        s[k].rms :=   FGet.LongReal(rd);
        (* Wr.PutText(stderr, Fmt.LongReal(s[k].rms) & "\n"); *)
        FGet.Skip(rd);
     END;
     RETURN su
   END;  
 END ReadFreqPower;   


PROCEDURE WriteMeanFreqPower(
  outName: TEXT;
  READONLY accMean: LONGS;
  READONLY freq: LONGS;
  READONLY lambda: LONGS;
 )=
  BEGIN 
    WITH 
      outFile = outName & ".fpw",
      wr = FileWr.Open(outFile)
    DO 
      Wr.PutText(stderr, "writing file " & outFile & " ...\n");
      Wr.PutText(wr, "# freq  lambda    mean \n");
      FOR k:= 0 TO LAST(accMean) DO
        IF lambda[k] >= 0.0d0 THEN
          FPut.LongReal(wr, freq[k], 2, Fmt.Style.Fix); 
          Wr.PutText(wr, "    ");
          FPut.LongReal(wr, lambda[k], 2, Fmt.Style.Fix); 
          Wr.PutText(wr, "    ");
          FPut.LongReal(wr, sqrt(accMean[k]), 4, Fmt.Style.Fix); 
          Wr.PutText(wr, "\n");
        END
      END;
      Wr.Flush(wr);
    END;
  END WriteMeanFreqPower;


PROCEDURE GetOptions(): Options =
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        pp.getKeyword("-outName");
        o.outName := pp.getNext();
        
        pp.getKeyword("-nFreqs");
        o.nFreqs := pp.getNextInt(1);
        
        pp.finish();      
        
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: BITSpectrum \\\n");
          Wr.PutText(stderr, "  -nFreqs CARDINAL  \\\n");
          Wr.PutText(stderr, "  -outName NAME\n");
          Process.Exit(1);
      END;
    END;
    RETURN o
    
  END GetOptions;


BEGIN
  Main()
END BITSpectrum.