/* Writes the Fourier power spectrum OF a curve. */
/* Last edited on 2015-01-20 16:48:13 by stolfilocal */

#include <stdio.h>
#include <math.h>

#include <argparser.h>

#include <vsignal.h>
#include <pz_r3_chain.h>
#include <pz_double_chain.h>
#include <pz_fourier.h>
#include <pz_types.h>
 

TYPE

  Options = struct ??? {
    char *inFile;         /* Input curve sample file (with extension). */
    char *outName;        /* Output spectrum file name (without extension). */
    double lambdaMin;      /* Min scale length for spectrum computation. */
    bool_t uniformize; /* TRUE to uniformize velocity before sampling */
  } ???;
    

int main(int argc, char **argv 
  )
  {

    { /* with */
      ??? o = pz_get_options(int argc, char **argv);
      ??? rp = pz_r3_chain_read(open_read(o.inFile, TRUE), header_only := FALSE, recenter := pz_ctr_SMPS);
      ??? p = rp.c^;
      ??? m = (p.ne);
      ??? c = pz_curve.New(m);
      /* do */ 
      c.setSamples(p);
      fprintf(stderr, "original period == %s\n", FLR(c.tPeriod,8,5) );
      if (( o.uniformize )){ 
        c.uniformize();
        fprintf(stderr, "new period (length) == %s\n", FLR(c.tPeriod,8,5) );;
      };

      { /* with */ 
        ??? pwr = pz_power_spectrum_compute(c, o.lambdaMin)^;
        /* do */
        WritePowerSpectrum(pwr, c.tPeriod, o);;
      };
    }
  } /* Main */
  

Options  pz_get_options(int argc, char **argv 
  )
  VAR o: Options;
  {
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
    

    { /* with */
      /* do */
      TRY
        argparser_get_keyword(pp, "-inFile");
        o.inFile = argparser_get_next(pp);
        
        argparser_get_keyword(pp, "-outName");
        o.outName = argparser_get_next(pp);
        
        argparser_get_keyword(pp, "-lambdaMin");
        o.lambdaMin = argparser_get_next_double(pp);

        o.uniformize = argparser_keyword_present(pp, "-uniformize");
      
        argparser_finish(pp);      
        
      EXCEPT                                                            
      | ParseParams.Error ==>                                              
          fprintf(stderr, "Usage: pz_spectrum \\\n");
          fprintf(stderr, "  -inFile NAME \\\n");
          fprintf(stderr, "  -outName NAME  \\\n");
          fprintf(stderr, "  -lambdaMin NUM  \\\n");
          fprintf(stderr, "  [ -uniformize ] \n");
        Process.Exit(1);;
      };;
    };
    return o
  } /* GetOptions */


void pz_spectrum_write(
    double_vec_t *pwr,
    double tPeriod,
    Options *o 
  )
  {
    WriteFreqPower(o.outName, pwr, tPeriod);
    WriteBandPower(o.outName, pwr, tPeriod, o.lambdaMin, 8, 2.0e0)
  } /* pz_spectrum_write */


double_vec_t *pz_power_spectrum_compute(
    pz_curve.T *c, 
    double lambdaMin 
  )
  /*
    Returns a vector "pwr[f]" with the power (square of amplitude)
    for each frequency "f", from 0 to "fMax > c.m DIV 2".
    Elements "pwr[0] and "pwr[fMax]" reflect only one term of the
    the series, all other elements include two conjugate terms. */
  VAR p: double;
  {

    { /* with */
      ??? m = c.m;
      ??? estLength = 1.5e0 * c.tPeriod;
      ??? nMin = MAX(2*m, CEILING(8.0e0*estLength/lambdaMin));
      ??? n = PowerOfTwo(nMin);
      ??? fMax = n DIV 2;
      ??? t = pz_double_chain_new(n)^;
      ??? u = pz_r3_chain_new(n)^;
      ??? v = pz_r3_chain_new(n)^;
      ??? r = pz_double_chain_new(n)^;
      ??? a = NEW(REF pz_fourier.T, 3, n)^;
      ??? rpwr = double_vec_new(fMax+1), pwr = rpwr^;
      /* do */
      for (i = 0; i <= n-1 ; i++){ t[i] = c.tPeriod * FLOAT(i,double)/FLOAT(n,double) ;};
      c.sample(t, u, v, r);
      pz_r3_chain_fourier_transform(u, a);
      

      { /* with */
        ??? ax = a[0], ay = a[1], az = a[2];
        /* do */
        for (f = 0; f <= fMax ; f++){

          { /* with */
            ??? tx = ax[f], ty = ay[f], tz = ax[f];
            /* do */
            p = tx*tx + ty*ty + tz*tz;
          };
          if (( f > 0 ) ANDAND ( f < fMax )){

            { /* with */
              ??? tx = ax[n-f], ty = ay[n-f], tz = az[n-f];
              /* do */
              p = p + tx*tx + ty*ty + tz*tz;
            };
          };
          pwr[f] = p;
        };
      };
      return rpwr;
    }
  } /* pz_power_spectrum_compute */


unsigned  PowerOfTwo(
    unsigned n 
  ) 
  /* Smallest power of two not less than "n" */
  VAR m: unsigned := 1;
  {
    while ( m < n ){ m = m+m ;};
    return m
  } /* PowerOfTwo */


void WriteFreqPower(
    char *name, 
    double_vec_t *pwr, 
    double tPeriod 
  )
  {  

    { /* with */ 
      ??? fMax = (pwr.ne - 1);
      ??? fileName = txtcat(name , "-fpw.plt");
      ??? wr = open_write(fileName, TRUE);
      /* do */ 
      fprintf(stderr, "writing file %s ...\n", fileName );
      fprintf(wr, "# freq   lambda     rms  \n");
      fprintf(wr, "# ----  ------- -------  \n");
      for (f = 1; f <= fMax ; f++){

        { /* with */
          ??? lambda = tPeriod/((double)f);
          /* do */
          FPut.Int(wr, f, 6); 
          fprintf(wr, "  ");
          FPut.LongReal(wr, lambda, 4, style = Fmt.Style.Fix); 
          fprintf(wr, "  ");
          FPut.LongReal(wr, Math.sqrt(pwr[f]), 4, style = Fmt.Style.Fix); 
          fprintf(wr, "\n");;
        };;
      };
      fflush(wr);
    }    
  } /* WriteFreqPower */


void WriteBandPower(
    char *name, 
    double_vec_t *pwr,
    double length,
    double lambdaMin,
    unsigned nScales,
    double scaleRatio 
  )
  VAR  
    sum: double;
    nb: unsigned;
  {  

    { /* with */ 
      ??? fMax = (pwr.ne - 1);
      ??? fileName = txtcat(name , "-bpw.plt");
      ??? wr = open_write(fileName, TRUE);
      /* do */ 
      fprintf(stderr, "writing file %s ...\n", fileName );
      fprintf(wr, "# lambdaLo lambdaHi    rms  \n");
      lambdaLo = 0.0e0;
      lambdaHi = lambdaMin;
      sum = 0.0e0;
      nb = 0;
      for (f = fMax; f <= 0 BY -1 ; f++){
        if (( f == 0 )){ 
          lambda = (double.ne - 1) 
        }else{
          lambda = length/((double)f);
        };
        while ( lambda > lambdaHi ) ANDAND ( lambdaLo < length ){ 
          FPut.LongReal(wr, lambdaLo, 4, style = Fmt.Style.Fix); 
          fprintf(wr, "  ");
          FPut.LongReal(wr, lambdaHi, 4, style = Fmt.Style.Fix); 
          fprintf(wr, "  ");
          FPut.LongReal(wr, Math.sqrt(sum), 4, style = Fmt.Style.Fix); 
          fprintf(wr, "\n");
          nb++;
          lambdaLo = lambdaHi; 
          if (( nb > nScales )){ 
            lambdaHi = length 
          }else{ 
            lambdaHi = scaleRatio*lambdaLo;
          };
          sum = 0.0e0;;
        };
        sum = sum + pwr[f];
      };
      fflush(wr);
    }    
  } /* WriteBandPower */


char *FLR(
    double x, 
    unsigned width,
    unsigned prec 
  )
  {
    return Fmt.Pad(Fmt.LongReal(x, prec = prec, style = Fmt.Style.Fix), width)
  } /* FLR */

{
/*
  Copyright © 2001 Universidade Estadual de Campinas (UNICAMP).
  Authors: Helena C. G. Leitão and Jorge Stolfi.
  
  This file can be freely distributed, used, and modified, provided
  that this copyright and authorship notice is preserved, and that any
  modified versions are clearly marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the authors nor their employers chall be held
  responsible for any losses or damages that may result from its use.
*/