package nxm.dsp.prim;

import static java.lang.Math.cos;
import static java.lang.Math.max;
import static java.lang.Math.min;
import static java.lang.Math.rint;
import static java.lang.Math.sin;
import static nxm.sys.lib.StringUtil.isNumber;
import static nxm.sys.lib.XmUtil.interpd;

import nxm.sys.inc.Units;
import nxm.sys.lib.Convert;
import nxm.sys.lib.Data;
import nxm.sys.lib.DataFile;
import nxm.sys.lib.MidasException;
import nxm.sys.lib.Primitive;

/** Double-precision phase/frequency/chirp-controlled oscillator. This is a
    port of the X-Midas DCPO command (dcpo.for in X-Midas 4.8.2). Although
    consideration was given to adding the features contained herein into
    the WAVEFORM command, a number of users migrating from X-Midas to NeXtMidas
    had a need for a direct port of the X-Midas version preserving some of the
    "features" contained there in. This is not dissimilar from the way X-Midas
    backed out their revised (C++) DPOC from X-Midas 4.4 in favor of the
    older version from X-Midas 3.6. <br>
    <br>
    Note that (as with X-Midas) the input files must be "files" with
    random-access support and can not be pipes.

    @author  J. Czechowski IV
    @since   NeXtMidas 3.3.0
 */
public class dpco extends Primitive {
  private DataFile chirpFile;    // Chirp input file (null=use chirpValue)
  private DataFile freqFile;     // Freq. input file (null=use freqValue)
  private DataFile phaseFile;    // Phase input file (null=use phaseValue)
  private double   chirpValue;   // Chirp input value
  private double   freqValue;    // Freq. input value
  private double   phaseValue;   // Phase input value
  private int      chirpNGot;    // Number of elements in chirp input buffer
  private int      freqNGot;     // Number of elements in freq. input buffer
  private int      phaseNGot;    // Number of elements in phase input buffer
  private Data     chirpBuf;     // Chirp input buffer
  private Data     freqBuf;      // Freq. input buffer
  private Data     phaseBuf;     // Phase input buffer

  private DataFile outSin;       // Sine    output file
  private DataFile outCos;       // Cosine  output file
  private DataFile outCmplx;     // Complex output file
  private DataFile outPhase;     // Phase   output file
  private DataFile outFreq;      // Freq.   output file
  private Data     sinOut;       // Sine    output buffer
  private Data     cosOut;       // Cosine  output buffer
  private Data     cmplxOut;     // Complex output buffer (complex data!)
  private Data     phaseOut;     // Phase   output buffer
  private Data     freqOut;      // Freq.   output buffer

  private double   xStart;       // x-axis start value
  private double   xDelta;       // x-axis delta
  private double   start;        // start value
  private double   end;          // end value
  private double   chirpStart;   // Chirp start
  private double   chirpEnd;     // Chirp end
  private double   freqStart;    // Freq. start
  private double   freqEnd;      // Freq. end
  private double   freqInst;     // Freq. instantaneous
  private double   freqInstOld;  // Freq. instantaneous (old value)
  private double   phaseStart;   // Phase start
  private double   phaseEnd;     // Phase end
  private double   phase;        // Phase value
  private double   phaseOld;     // Phase value (old value)

  private double   ntogo;        // Number of points still to do (entire file)
  private int      ntodo;        // Number of points to do on this pass
  private boolean  abscissa;     // Use abcissa addressing?
  private boolean  first;        // First transfer length?

  /** default constructor */
  public dpco () {}

  @Override
  public int open () {
    // non-file-name arguments
    String addressing = MA.getS("ADDRESSING");  // ABSCISSA->1, INDEX->0

    xDelta    = MA.getD("DELTA");
    xStart    = MA.getD("START");
    abscissa  = addressing.equals("ABSCISSA") || addressing.startsWith("AB");

    double numPoints     = rint(MA.getD("ELEM"));
    double origNumPoints = numPoints;
    double xEnd          = xStart + (numPoints-1.0)*xDelta;

    // pick up input values/files
    String _inFreq  = MA.getO("INFREQ",  0).toString();
    String _inPahse = MA.getO("INPHASE", 0).toString();
    String _inChirp = MA.getO("INCHIRP", 0).toString();

    if (isNumber(_inFreq)) {
      freqFile  = null;
      freqValue = Convert.o2d(_inFreq);
    }
    else {
      freqFile  = MA.getDataFile("INFREQ",  "1000", "S#,-SP", DataFile.INPUT);
      freqValue = 0.0;
    }

    if (isNumber(_inPahse)) {
      phaseFile  = null;
      phaseValue = Convert.o2d(_inPahse);
    }
    else {
      phaseFile  = MA.getDataFile("INPHASE", "1000", "S#,-SP", DataFile.INPUT);
      phaseValue = 0.0;
    }

    if (isNumber(_inChirp)) {
      chirpFile  = null;
      chirpValue = Convert.o2d(_inChirp);
    }
    else {
      chirpFile  = MA.getDataFile("INCHIRP", "1000", "S#,-SP", DataFile.INPUT);
      chirpValue = 0.0;
    }

    // open any input files, and get mutual overlap with output range
    if (freqFile != null) {
      freqFile.open();
      if (abscissa) {
        if (freqFile.getXDelta() < xDelta) {
          throw new MidasException("Frequency file increment too small");
        }
        xStart = max(freqFile.getXStart(), xStart);
        xEnd   = min(freqFile.getXStart() + freqFile.getSize()*freqFile.getXDelta(), xEnd);
        // (we can interpolate 1 pt beyond where file actually ends)
      }
      else {
        freqFile.setXStart(xStart);
        freqFile.setXDelta(xDelta);
        numPoints = min(freqFile.getSize(), numPoints);
      }
    }

    if (phaseFile != null) {
      phaseFile.open();
      if (abscissa) {
        if (phaseFile.getXDelta() < xDelta) {
          throw new MidasException("Phase file increment too small");
        }
        xStart = max(phaseFile.getXStart(), xStart);
        xEnd   = min(phaseFile.getXStart() + phaseFile.getSize()*phaseFile.getXDelta(), xEnd);
      }
      else {
        phaseFile.setXStart(xStart);
        phaseFile.setXDelta(xDelta);
        numPoints = min(phaseFile.getSize(), numPoints);
      }
    }

    if (chirpFile != null) {
      chirpFile.open();
      if (abscissa) {
        if (chirpFile.getXDelta() < xDelta) {
          throw new MidasException("Chirp file increment too small");
        }
        xStart = max(chirpFile.getXStart(), xStart);
        xEnd   = min(chirpFile.getXStart() + chirpFile.getSize()*chirpFile.getXDelta(), xEnd);
      }
      else {
        chirpFile.setXStart(xStart);
        chirpFile.setXDelta(xDelta);
        numPoints = min(chirpFile.getSize(), numPoints);
      }
    }

    // if abscissa, go back and limit number of pts according to final start/end
    if (abscissa) {
      if ((freqFile != null) || (phaseFile != null) || (chirpFile != null)) {
        // when using input files truncate the size of the input files
        numPoints = (int)((xEnd-xStart)/xDelta + 1);
      }
      else {
        // without input files we round the to the nearest
        numPoints = rint((xEnd-xStart)/xDelta + 1);
      }
    }
    if (numPoints <= 0) {
      throw new MidasException("Input files outside output abscissa range");
    }
    if (numPoints != origNumPoints) {
      M.warning("Number of points reduced from "+origNumPoints+" to "+numPoints);
    }


    // output files ?
    outSin   = openOutputDataFile("OUTSIN",     "SF", numPoints);
    outCos   = openOutputDataFile("OUTCOS",     "SF", numPoints);
    outPhase = openOutputDataFile("OUTPHASE",   "SF", numPoints);
    outFreq  = openOutputDataFile("OUTFREQ",    "SF", numPoints);
    outCmplx = openOutputDataFile("OUTCOMPLEX", "CF", numPoints);

    outFreq.setXStart(outFreq.getXStart() + outFreq.getXDelta()); // skip 1st point

    // check for degenerate output case
    if (!(outSin.isOpen() || outCos.isOpen() || outCmplx.isOpen() ||
          outPhase.isOpen() || outFreq.isOpen())) {
      M.warning("No output files specified");
      // Unlike X-Midas which would continue on in this situation, just without
      // writing anything out, we will return FINISH indicating that there is
      // nothing left to do.
      return FINISH;
    }

    // setup GRAB parameters; or prefill buffer if using constant
    xfer     = MA.getL("/TL", 16384);
    freqBuf  = new Data("SD", xfer, freqValue);
    phaseBuf = new Data("SD", xfer, (phaseValue/360.0));
    chirpBuf = new Data("SD", xfer, chirpValue);
    freqOut  = new Data("SD", xfer);
    sinOut   = new Data("SD", xfer); // X-Midas matches 'new Data("SF", xfer*2)'
    cosOut   = new Data("SD", xfer); // X-Midas matches 'new Data("SF", xfer*2)'
    phaseOut = new Data("SD", xfer);
    cmplxOut = new Data("CD", xfer);

    if (freqFile != null) {
      freqFile.buffer = freqFile.getDataBuffer(xfer, 'D');
      freqEnd         = xStart-1;       // force initial read
    }

    if (phaseFile != null) {
      phaseFile.buffer = phaseFile.getDataBuffer(xfer, 'D');
      phaseEnd         = xStart-1;       // force initial read
    }

    if (chirpFile != null) {
      chirpFile.buffer = chirpFile.getDataBuffer(xfer, 'D');
      chirpEnd         = xStart-1;       // force initial read
    }

    // init variables, setup counters, synchronize
    phase       = 0.0;
    phaseOld    = 0.0;
    freqInst    = 0.0;
    freqInstOld = 0.0;
    ntogo       = numPoints;
    ntodo       = (int)min(ntogo, xfer);
    start       = xStart - ntodo*xDelta;
    end         = xStart - xDelta;
    first       = true;

    return NORMAL;
  }

  @Override
  public int process () {
    if (ntogo <= 0) return FINISH;

    // update counters, start/end range
    start = start + ntodo*xDelta;
    end   = end   + ntodo*xDelta;
    ntodo = (int)min(ntogo, xfer);
    ntogo = ntogo - ntodo;

    // grab new frequency/phase/chirp buffer if current range outside buffer;
    // then interpolate over current range
    //
    // Do NOT change the read calls to do a pipe-friendly read. The calls used
    // here are intentional, and these files, if present, are like array
    // parameters and not like streams of data. DPCO needs to perform random
    // access on these arrays.
    if (freqFile != null) {                               // frequency
      if (freqEnd < end) {
        int index = (int)( (start-freqFile.getXStart())/freqFile.getXDelta() );
        if (abscissa) {
          freqNGot = freqFile.read(freqFile.buffer, index, xfer);
        }
        else {
          freqNGot = freqFile.read(freqBuf, index, xfer);
        }
        freqStart = freqFile.getXStart() + index*freqFile.getXDelta();
        freqEnd   = freqStart + (freqNGot-1)*freqFile.getXDelta();
      }
      if (abscissa) interpd(freqNGot, freqStart, freqFile.getXDelta(),
                            freqFile.buffer, ntodo, start, xDelta, freqBuf, 1, 1);
    }

    if (phaseFile != null) {                               // phase
      if (phaseEnd < end) {
        int index = (int)( (start-phaseFile.getXStart())/phaseFile.getXDelta() );
        if (abscissa) {
          phaseNGot = phaseFile.read(phaseFile.buffer, index, xfer);

          double[] _phase = phaseFile.buffer.castD(true);
          for (int j = 0; j < phaseNGot; j++) {
            _phase[j] = _phase[j] / 360;
          }
          phaseFile.buffer.uncast(_phase, true);
        }
        else {
          phaseNGot = phaseFile.read(phaseBuf, index, xfer);

          double[] _phase = phaseBuf.castD(true);
          for (int j = 0; j < phaseNGot; j++) {
            _phase[j] = _phase[j] / 360;
          }
          phaseBuf.uncast(_phase, true);
        }
        phaseStart = phaseFile.getXStart() + index*phaseFile.getXDelta();
        phaseEnd   = phaseStart + (phaseNGot-1)*phaseFile.getXDelta();
      }
      if (abscissa) interpd(phaseNGot, phaseStart, phaseFile.getXDelta(),
                            phaseFile.buffer, ntodo, start, xDelta, phaseBuf, 1, 1);
    }


    if (chirpFile != null) {                               // chirp
      if (chirpEnd < end) {
        int index = (int)( (start-chirpFile.getXStart())/chirpFile.getXDelta() );
        if (abscissa) {
          chirpNGot = chirpFile.read(chirpFile.buffer, index, xfer);
        }
        else {
          chirpNGot = chirpFile.read(chirpBuf, index, xfer);
        }
        chirpStart = chirpFile.getXStart() + index*chirpFile.getXDelta();
        chirpEnd   = chirpStart + (chirpNGot-1)*chirpFile.getXDelta();
      }
      if (abscissa) interpd(chirpNGot, chirpStart, chirpFile.getXDelta(),
                            chirpFile.buffer, ntodo, start, xDelta, chirpBuf, 1, 1);
    }

    // now create a bufferful of phases
    double[] _chirpIn  = chirpBuf.castD(true);
    double[] _freqIn   = freqBuf.castD(true);
    double[] _phaseIn  = phaseBuf.castD(true);
    double[] _freqOut  = freqOut.castD(false);
    double[] _phaseOut = phaseOut.castD(false);
    for (int j = 0; j < ntodo; j++) {
      double phaseTemp;
      // Update phase with new phase, without unnecessarily losing any
      // mantissa precision.  If pbuf(j) were added directly to 'phase',
      // the lower-order bits in the mantissa of 'phase' would be lost.
      phaseTemp    = _phaseIn[j];
      phaseTemp    = phaseTemp - rint(phaseTemp);
      phase        = phase + phaseTemp - phaseOld;
      phase        = phase - rint(phase);        // i.e. MOD(phase,1)
      _phaseOut[j] = phase;
      phaseOld     = phaseTemp;

      // ..integrate chirp to get instantaneous frequency, then integrate frequency
      // ..and instantaneous frequency to update phase for next time.  Use trapezoidal
      // ..integration for instant. frequency, rectangular for chirp and frequency
      freqInst    = freqInst + _chirpIn[j]*xDelta;
      _freqOut[j] = _freqIn[j] + (freqInst+freqInstOld)*0.5;

      // Update phase, without unnecessarily losing any mantissa
      // precision (see earlier comment on pbuf(j)).
      phaseTemp   = _freqOut[j]*xDelta;
      phaseTemp   = phaseTemp - rint(phaseTemp);
      phase       = phase + phaseTemp;
      freqInstOld = freqInst;
    }
    chirpBuf.uncast(_chirpIn,  false);
    freqBuf.uncast(_freqIn,    false);
    phaseBuf.uncast(_phaseIn,  false);
    freqBuf.uncast(_freqOut,   outFreq.isOpen());
    // Note: keep _phaseOut cast as an array for a little bit longer

    // now what do we do with this nice phase?  => complex exponential if needed
    if (outCmplx.isOpen() || outSin.isOpen() || outCos.isOpen()) {
      double[] _cosOut   = cosOut.castD(false);
      double[] _sinOut   = sinOut.castD(false);
      double[] _cmplxOut = cmplxOut.castD(false);
      for (int j = 0; j < ntodo; j++) {
        _cosOut[j]       = cos(_phaseOut[j]*TWOPI);
        _sinOut[j]       = sin(_phaseOut[j]*TWOPI);
        _cmplxOut[j*2+0] = _cosOut[j];
        _cmplxOut[j*2+1] = _sinOut[j];
      }
      cosOut.uncast(_cosOut,     outCos.isOpen());
      sinOut.uncast(_sinOut,     outSin.isOpen());
      cmplxOut.uncast(_cmplxOut, outCmplx.isOpen());

      if (outSin.isOpen()  ) outSin.write(sinOut, ntodo);
      if (outCos.isOpen()  ) outCos.write(cosOut, ntodo);
      if (outCmplx.isOpen()) outCmplx.write( cmplxOut, ntodo);
    }

    // output phase, frequency if requested
    if (outPhase.isOpen()) {
      for (int j = 0; j < ntodo; j++) {
        if (_phaseOut[j] > 0.5) _phaseOut[j] = _phaseOut[j] - 1;
        _phaseOut[j] = _phaseOut[j]*360;
      }
      phaseBuf.uncast(_phaseOut, true);
      outPhase.write(phaseOut, ntodo);
    }
    else {
      phaseBuf.uncast(_phaseOut, false);
    }

    if (outFreq.isOpen()) {
      if (first) {
        freqOut.boff = freqOut.boff + freqOut.bpe;   // <-- shift start to second element
        outFreq.write(freqOut, ntodo-1);
        first = false;
        freqOut.boff = freqOut.boff - freqOut.bpe;   // <-- shift start back to first element
      }
      else {
        outFreq.write(freqOut, ntodo);
      }
    }

    return NORMAL;
  }

  @Override
  public int close () {
    // close files
    if (phaseFile != null) { phaseFile.close(); phaseFile = null; }
    if (freqFile  != null) { freqFile.close();  freqFile  = null; }
    if (chirpFile != null) { chirpFile.close(); chirpFile = null; }
    if (outSin    != null) { outSin.close();    outSin    = null; }
    if (outCos    != null) { outCos.close();    outCos    = null; }
    if (outCmplx  != null) { outCmplx.close();  outCmplx  = null; }
    if (outPhase  != null) { outPhase.close();  outPhase  = null; }
    if (outFreq   != null) { outFreq.close();   outFreq   = null; }

    return NORMAL;
  }

  /** X-Midas used some funny M$PROPAGATE stuff to combine some of this, but
      since that was an abuse of the propagate functionality (and non-intuitive
      to many users), it seems better to just put together a simple function
      that does what we need.
   */
  private DataFile openOutputDataFile (String name, String format, double size) {
    DataFile df = MA.getDataFile(name, "1000", format, DataFile.OUTPUT|DataFile.OPTIONAL);
    df.setSize(size);
    df.setXStart(xStart);
    df.setXDelta(xDelta);
    df.setXUnits(Units.TIME_S);
    df.open();
    return df;
  }
}
