/***************************************************/
/*! \class Mesh2D
    \brief Two-dimensional rectilinear waveguide mesh class.

    This class implements a rectilinear,
    two-dimensional digital waveguide mesh
    structure.  For details, see Van Duyne and
    Smith, "Physical Modeling with the 2-D Digital
    Waveguide Mesh", Proceedings of the 1993
    International Computer Music Conference.

    This is a digital waveguide model, making its
    use possibly subject to patents held by Stanford
    University, Yamaha, and others.

    Control Change Numbers: 
       - X Dimension = 2
       - Y Dimension = 4
       - Mesh Decay = 11
       - X-Y Input Position = 1

    by Julius Smith, 2000 - 2002.
    Revised by Gary Scavone for STK, 2002.
*/
/***************************************************/

#include "Mesh2D.h"
#include "SKINImsg.h"

namespace stk {

Mesh2D :: Mesh2D( unsigned short nX, unsigned short nY )
{
  if ( nX == 0.0 || nY == 0.0 ) {
    oStream_ << "Mesh2D::Mesh2D: one or more argument is equal to zero!";
    handleError( StkError::FUNCTION_ARGUMENT );
  }

  this->setNX( nX );
  this->setNY( nY );

  StkFloat pole = 0.05;
  unsigned short i;
  for ( i=0; i<NYMAX; i++ ) {
    filterY_[i].setPole( pole );
    filterY_[i].setGain( 0.99 );
  }

  for ( i=0; i<NXMAX; i++ ) {
    filterX_[i].setPole( pole );
    filterX_[i].setGain( 0.99 );
  }

  this->clearMesh();

  counter_ = 0;
  xInput_ = 0;
  yInput_ = 0;
}

Mesh2D :: ~Mesh2D( void )
{
}

void Mesh2D :: clear( void )
{
  this->clearMesh();

  unsigned short i;
  for ( i=0; i<NY_; i++ )
    filterY_[i].clear();

  for ( i=0; i<NX_; i++ )
    filterX_[i].clear();

  counter_ = 0;
}

void Mesh2D :: clearMesh( void )
{
  int x, y;
  for ( x=0; x<NXMAX-1; x++ ) {
    for ( y=0; y<NYMAX-1; y++ ) {
      v_[x][y] = 0;
    }
  }
  for ( x=0; x<NXMAX; x++ ) {
    for ( y=0; y<NYMAX; y++ ) {

      vxp_[x][y] = 0;
      vxm_[x][y] = 0;
      vyp_[x][y] = 0;
      vym_[x][y] = 0;

      vxp1_[x][y] = 0;
      vxm1_[x][y] = 0;
      vyp1_[x][y] = 0;
      vym1_[x][y] = 0;
    }
  }
}

StkFloat Mesh2D :: energy( void )
{
  // Return total energy contained in wave variables Note that some
  // energy is also contained in any filter delay elements.

  int x, y;
  StkFloat t;
  StkFloat e = 0;
  if ( counter_ & 1 ) { // Ready for Mesh2D::tick1() to be called.
    for ( x=0; x<NX_; x++ ) {
      for ( y=0; y<NY_; y++ ) {
        t = vxp1_[x][y];
        e += t*t;
        t = vxm1_[x][y];
        e += t*t;
        t = vyp1_[x][y];
        e += t*t;
        t = vym1_[x][y];
        e += t*t;
      }
    }
  }
  else { // Ready for Mesh2D::tick0() to be called.
    for ( x=0; x<NX_; x++ ) {
      for ( y=0; y<NY_; y++ ) {
        t = vxp_[x][y];
        e += t*t;
        t = vxm_[x][y];
        e += t*t;
        t = vyp_[x][y];
        e += t*t;
        t = vym_[x][y];
        e += t*t;
      }
    }
  }

  return e;
}

void Mesh2D :: setNX( unsigned short lenX )
{
  if ( lenX < 2 ) {
    oStream_ << "Mesh2D::setNX(" << lenX << "): Minimum length is 2!";
    handleError( StkError::WARNING ); return;
  }
  else if ( lenX > NXMAX ) {
    oStream_ << "Mesh2D::setNX(" << lenX << "): Maximum length is " << NXMAX << '!';
    handleError( StkError::WARNING ); return;
  }

  NX_ = lenX;
}

void Mesh2D :: setNY( unsigned short lenY )
{
  if ( lenY < 2 ) {
    oStream_ << "Mesh2D::setNY(" << lenY << "): Minimum length is 2!";
    handleError( StkError::WARNING ); return;
  }
  else if ( lenY > NYMAX ) {
    oStream_ << "Mesh2D::setNY(" << lenY << "): Maximum length is " << NXMAX << '!';
    handleError( StkError::WARNING ); return;
  }

  NY_ = lenY;
}

void Mesh2D :: setDecay( StkFloat decayFactor )
{
  if ( decayFactor < 0.0 || decayFactor > 1.0 ) {
    oStream_ << "Mesh2D::setDecay: decayFactor is out of range!";
    handleError( StkError::WARNING ); return;
  }

  int i;
  for ( i=0; i<NYMAX; i++ )
    filterY_[i].setGain( decayFactor );

  for (i=0; i<NXMAX; i++)
    filterX_[i].setGain( decayFactor );
}

void Mesh2D :: setInputPosition( StkFloat xFactor, StkFloat yFactor )
{
  if ( xFactor < 0.0 || xFactor > 1.0 ) {
    oStream_ << "Mesh2D::setInputPosition xFactor value is out of range!";
    handleError( StkError::WARNING ); return;
  }

  if ( yFactor < 0.0 || yFactor > 1.0 ) {
    oStream_ << "Mesh2D::setInputPosition yFactor value is out of range!";
    handleError( StkError::WARNING ); return;
  }

  xInput_ = (unsigned short) (xFactor * (NX_ - 1));
  yInput_ = (unsigned short) (yFactor * (NY_ - 1));
}

void Mesh2D :: noteOn( StkFloat frequency, StkFloat amplitude )
{
  // Input at corner.
  if ( counter_ & 1 ) {
    vxp1_[xInput_][yInput_] += amplitude;
    vyp1_[xInput_][yInput_] += amplitude;
  }
  else {
    vxp_[xInput_][yInput_] += amplitude;
    vyp_[xInput_][yInput_] += amplitude;
  }
}

void Mesh2D :: noteOff( StkFloat amplitude )
{
  return;
}

StkFloat Mesh2D :: inputTick( StkFloat input )
{
  if ( counter_ & 1 ) {
    vxp1_[xInput_][yInput_] += input;
    vyp1_[xInput_][yInput_] += input;
    lastFrame_[0] = tick1();
  }
  else {
    vxp_[xInput_][yInput_] += input;
    vyp_[xInput_][yInput_] += input;
    lastFrame_[0] = tick0();
  }

  counter_++;
  return lastFrame_[0];
}

StkFloat Mesh2D :: tick( unsigned int )
{
  lastFrame_[0] = ((counter_ & 1) ? this->tick1() : this->tick0());
  counter_++;
  return lastFrame_[0];
}

const StkFloat VSCALE = 0.5;

StkFloat Mesh2D :: tick0( void )
{
  int x, y;
  StkFloat outsamp = 0;

  // Update junction velocities.
  for (x=0; x<NX_-1; x++) {
    for (y=0; y<NY_-1; y++) {
      v_[x][y] = ( vxp_[x][y] + vxm_[x+1][y] + 
		  vyp_[x][y] + vym_[x][y+1] ) * VSCALE;
    }
  }    

  // Update junction outgoing waves, using alternate wave-variable buffers.
  for (x=0; x<NX_-1; x++) {
    for (y=0; y<NY_-1; y++) {
      StkFloat vxy = v_[x][y];
      // Update positive-going waves.
      vxp1_[x+1][y] = vxy - vxm_[x+1][y];
      vyp1_[x][y+1] = vxy - vym_[x][y+1];
      // Update minus-going waves.
      vxm1_[x][y] = vxy - vxp_[x][y];
      vym1_[x][y] = vxy - vyp_[x][y];
    }
  }    

  // Loop over velocity-junction boundary faces, update edge
  // reflections, with filtering.  We're only filtering on one x and y
  // edge here and even this could be made much sparser.
  for (y=0; y<NY_-1; y++) {
    vxp1_[0][y] = filterY_[y].tick(vxm_[0][y]);
    vxm1_[NX_-1][y] = vxp_[NX_-1][y];
  }
  for (x=0; x<NX_-1; x++) {
    vyp1_[x][0] = filterX_[x].tick(vym_[x][0]);
    vym1_[x][NY_-1] = vyp_[x][NY_-1];
  }

  // Output = sum of outgoing waves at far corner.  Note that the last
  // index in each coordinate direction is used only with the other
  // coordinate indices at their next-to-last values.  This is because
  // the "unit strings" attached to each velocity node to terminate
  // the mesh are not themselves connected together.
  outsamp = vxp_[NX_-1][NY_-2] + vyp_[NX_-2][NY_-1];

  return outsamp;
}

StkFloat Mesh2D :: tick1( void )
{
  int x, y;
  StkFloat outsamp = 0;

  // Update junction velocities.
  for (x=0; x<NX_-1; x++) {
    for (y=0; y<NY_-1; y++) {
      v_[x][y] = ( vxp1_[x][y] + vxm1_[x+1][y] + 
		  vyp1_[x][y] + vym1_[x][y+1] ) * VSCALE;
    }
  }

  // Update junction outgoing waves, 
  // using alternate wave-variable buffers.
  for (x=0; x<NX_-1; x++) {
    for (y=0; y<NY_-1; y++) {
      StkFloat vxy = v_[x][y];

      // Update positive-going waves.
      vxp_[x+1][y] = vxy - vxm1_[x+1][y];
      vyp_[x][y+1] = vxy - vym1_[x][y+1];

      // Update minus-going waves.
      vxm_[x][y] = vxy - vxp1_[x][y];
      vym_[x][y] = vxy - vyp1_[x][y];
    }
  }

  // Loop over velocity-junction boundary faces, update edge
  // reflections, with filtering.  We're only filtering on one x and y
  // edge here and even this could be made much sparser.
  for (y=0; y<NY_-1; y++) {
    vxp_[0][y] = filterY_[y].tick(vxm1_[0][y]);
    vxm_[NX_-1][y] = vxp1_[NX_-1][y];
  }
  for (x=0; x<NX_-1; x++) {
    vyp_[x][0] = filterX_[x].tick(vym1_[x][0]);
    vym_[x][NY_-1] = vyp1_[x][NY_-1];
  }

  // Output = sum of outgoing waves at far corner.
  outsamp = vxp1_[NX_-1][NY_-2] + vyp1_[NX_-2][NY_-1];

  return outsamp;
}

void Mesh2D :: controlChange( int number, StkFloat value )
{
#if defined(_STK_DEBUG_)
  if ( Stk::inRange( value, 0.0, 128.0 ) == false ) {
    oStream_ << "Mesh2D::controlChange: value (" << value << ") is out of range!";
    handleError( StkError::WARNING ); return;
  }
#endif

  StkFloat normalizedValue = value * ONE_OVER_128;
  if ( number == 2 ) // 2
    this->setNX( (unsigned short) (normalizedValue * (NXMAX-2) + 2) );
  else if ( number == 4 ) // 4
    this->setNY( (unsigned short) (normalizedValue * (NYMAX-2) + 2) );
  else if ( number == 11 ) // 11
    this->setDecay( 0.9 + (normalizedValue * 0.1) );
  else if ( number == __SK_ModWheel_ ) // 1
    this->setInputPosition( normalizedValue, normalizedValue );
#if defined(_STK_DEBUG_)
  else {
    oStream_ << "Mesh2D::controlChange: undefined control number (" << number << ")!";
    handleError( StkError::WARNING );
  }
#endif
}

} // stk namespace