#ifndef POWERSERIESWD_H
#define POWERSERIESWD_H

#include <cassert>
using namespace std;

#include <powerseries1D.h>

/*!
\brief W powerseries equations with N dimension each. 

This class is a power series in N dimensions.  It has
 a current state held in currentP.  By setting the
 current state this class represents one of the
 W equations.

\par Example

\verbatim
  uintc N(4);
  uintc W(2);
  double az[W*N] = 
  { 
    1.0, 2.0, 3.0, 4.0,
    3.0, -7.0, 4.5, 0.0
  };

  powerseriesWD<double> f(az,W,N);

  double y;
  f.currentSet(0);
  cout << "Expect 10" << endl;
  y=f(1.0);
\endverbatim
*/
template< typename T >
class powerseriesWD
{
public:

  /** The matrix representing the power series equations. */
  T * const mat;

  /** W equations. (rows) */
  uintc W;
  /** N dimension in each equation. (columns) */
  uintc N;

  /** The current power series. */
  powerseries1D<T> currentP;
  
  /** Set the current equation. */
  void currentSet(uintc i);

  /** Supply a matrix representing W power series equations. */
  powerseriesWD( T * const mat_, uintc W_, uintc N_ );

  /** Evalute the power series at t. */
  void operator()(T & val, T const t) const
    { currentP(val,t); }

  /** Evalute the power series at t. */
  T const operator()(T const t) const
    { return currentP(t); }

  
};

template< typename T >
void powerseriesWD<T>::currentSet(uintc i)
{
  assert( (i==0) || (i<N) );
  currentP.construct(mat+i*N,N);
}



template< typename T >
powerseriesWD<T>::powerseriesWD
( 
  T * const mat_, 
  uintc W_, 
  uintc N_
)
  : mat(mat_), W(W_), N(N_)
{ 
  currentP.construct(mat_,N);
}


#endif