MODULE BITSpectrum EXPORTS Main;

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

(* Last edited on 1998-11-09 19:16:43 by hcgl *)

IMPORT  ParseParams, Process, FileRd, Wr, FileWr, Rd;
IMPORT  OSError, Thread, Stdio, Fmt, FPut, FileFmt, FGet;
FROM Math IMPORT sqrt;
FROM Stdio IMPORT stderr;
FROM PZTypes IMPORT LONG, LONGS, NAT, INT, BOOL;
 
<* FATAL Wr.Failure, Thread.Alerted, OSError.E,  Rd.Failure  *>

CONST MaxBands = 256;

TYPE
  LONG = LONGREAL;
  
  BandData = RECORD
    lo: LONG;
    hi: LONG;
    rms: LONG
  END;
  
  Spectrum = ARRAY OF BandData;
  
  MismatchOptions = PZMismatch.Options;
  Options = RECORD
      inName: TEXT;      (* Prefix for input file names. *)
      nFiles: CARDINAL;  (* Number of files to read. *)
      outName: TEXT;     (* Prefix for output file names.  *)
      tag: TEXT;         (* Suffix for input and output file names. *)
  END;
  (*
    This program reads one or more band-power spectrum files with names
    of the form `xxx-iiiiii-yyy.bpw' where `xxx' is "o.inName", 
    `yyy' is "o.tag", and `iiiiii' is a six-digit serial number.
    
    Each file must contain sero or more records of the form "LAMBDALO
    LAMBDAHI RMS" where "RMS" is the square root of the total power in
    the wavelength range "[LAMBDALO _ LAMBDAHI)". All spectra must
    contain the same bands in the same order, except that some bands
    at the end may be missing in some files.
    
    The program computes the power in each band, averaged over all
    spectra that include that band, and writes the square root of the
    result as a single file `zzz-mean-yyy.bpw' where `zzz' is
    "o.outName". *)
  
PROCEDURE Main() =
  VAR 
     k : CARDINAL;
  BEGIN
    WITH
      o = GetOptions(),
      nFiles = o.nFiles,
      tag = o.tag,
      nSpectra = NEW(REF  LONGS, MaxBands)^,
      accPower = NEW(REF LONGS, MaxBands)^,
      lambdaLo = NEW(REF LONGS, MaxBands)^,
      lambdaHi = NEW(REF LONGS, MaxBands)^,
      accMean = NEW(REF LONGS,MaxBands)^
    DO 
      FOR k := 0 TO MaxBands-1 DO 
        nSpectra[k] := 0.0d0; accPower[k] := 0.0d0;
        lambdaLo[k] := 0.0d0; lambdaHi[k] := 0.0d0; 
      END;
      FOR i := 0 TO nFiles-1 DO
        WITH 
          inFile = o.inName & "-" & Fmt.Pad(Fmt.Int(i), 6, '0') & "-" & tag & ".bpw"
        DO 
          ReadAndAccumulateSpectrum(inFile, nSpectra, accPower, lambdaLo, lambdaHi)
        END 
      END;
      k:= 0; 
      WHILE  nSpectra[k] > 0.0d0 DO 
        accMean[k] := accPower[k]/nSpectra[k];
        k := k+1
      END;
      WITH 
        outFile = o.outName & "-mean-" & tag & ".bpw"
      DO 
        WriteMeanBandPower(outFile, accMean, lambdaLo, lambdaHi)
      END;   
    END 
  END Main;

 
PROCEDURE ReadAndAccumulateSpectrum(
    fileName: TEXT;
    VAR nSpectra: LONGS;
    VAR accPower: LONGS;
    VAR lambdaLo: LONGS;
    VAR lambdaHi: LONGS;
  ) =
 BEGIN  
    WITH 
      s = ReadBandPower(fileName)^,
      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;
        IF lambdaHi[k] = 0.0d0 THEN
          lambdaLo[k] := s[k].lo;
          lambdaHi[k] := s[k].hi;
          <* ASSERT lambdaLo[k] < lambdaHi[k] *>
          <* ASSERT lambdaLo[k] >= 0.0d0 *>
        ELSE
          <* ASSERT lambdaLo[k] = s[k].lo *>
          <* ASSERT lambdaHi[k] = s[k].hi *>
        END
      END; 
    END; 
  END ReadAndAccumulateSpectrum;


PROCEDURE  ReadBandPower(fileName : TEXT): REF Spectrum = 
  VAR k: CARDINAL;
  BEGIN  
    WITH 
      rd = FileRd.Open(fileName),
      s = NEW(REF Spectrum, MaxBands)^,
      cmt = FileFmt.ReadComment(rd, '#')
    DO 
      Wr.PutText(stderr, "reading file " & fileName & " ...\n");
      Wr.PutText(stderr, cmt); Wr.PutText(stderr, "\n"); 
      k:= 0;
      REPEAT 
        FGet.Skip(rd);
        s[k].lo :=   FGet.LongReal(rd); 
        Wr.PutText(stderr, "Low lambda: " & Fmt.LongReal(s[k].lo) & "\n");
        s[k].hi :=   FGet.LongReal(rd);
        Wr.PutText(stderr, "Hi lambda: " & Fmt.LongReal(s[k].hi) & "\n");
        s[k].rms :=   FGet.LongReal(rd);
        Wr.PutText(stderr, "sqrt(power) :" & Fmt.LongReal(s[k].rms) & "\n");
        FGet.Skip(rd);
        k := k + 1
      UNTIL Rd.EOF(rd);
      WITH 
        w = NEW(REF Spectrum, k)
      DO 
        w^ := SUBARRAY(s, 0, k);
        RETURN w 
      END
    END;  
 END ReadBandPower;   


PROCEDURE WriteMeanBandPower(
  outFile: TEXT;
  READONLY accMean: LONGS;
  READONLY lambdaLo: LONGS;
  READONLY lambdaHi: LONGS;
 )=
  BEGIN 
    WITH 
      wr = FileWr.Open(outFile)
    DO 
      Wr.PutText(stderr, "writing file " & outFile & " ...\n");
      Wr.PutText(wr, "# lambdaLo lambdaHi    mean \n");
      FOR k:= 0 TO LAST(accMean) DO
        IF lambdaHi[k] > 0.0d0 THEN
          FPut.LongReal(wr, lambdaLo[k], 2, Fmt.Style.Fix); 
          Wr.PutText(wr, "    ");
          FPut.LongReal(wr, lambdaHi[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 WriteMeanBandPower;


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


BEGIN
  Main()
END BITSpectrum.