source: Sophya/trunk/SophyaLib/NTools/fftservintf.cc@ 1402

#include "fftservintf.h"


/*!
  \class SOPHYA::FFTServerInterface
  \ingroup NTools
  Defines the interface for FFT (Fast Fourier Transform) operations.
*/

/* --Methode-- */
FFTServerInterface::FFTServerInterface(string info)
{
  _info = info;
  _fgnorm = true;
}

/* --Methode-- */
FFTServerInterface::~FFTServerInterface()
{
}

// ----------------- Transforme pour les double -------------------

/* --Methode-- */
void FFTServerInterface::FFTForward(TArray< complex<r_8> > const &, TArray< complex<r_8> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray...) Unsupported operation !");
}

/* --Methode-- */
void FFTServerInterface::FFTBackward(TArray< complex<r_8> > const &, TArray< complex<r_8> > &)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray...) Unsupported operation !");
}

/* --Methode-- */
void FFTServerInterface::FFTForward(TArray< r_8 > const &, TArray< complex<r_8> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray...) Unsupported operation !");
}

/* --Methode-- */
void FFTServerInterface::FFTBackward(TArray< complex<r_8> > const &, TArray< r_8 > &, bool)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray...) Unsupported operation !");
}


// ----------------- Transforme pour les float -------------------

/* --Methode-- */
void FFTServerInterface::FFTForward(TArray< complex<r_4> > const &, TArray< complex<r_4> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray r_4 ... ) Unsupported operation !");
}

/* --Methode-- */
void FFTServerInterface::FFTBackward(TArray< complex<r_4> > const &, TArray< complex<r_4> > &)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray r_4 ... ) Unsupported operation !");
}

/* --Methode-- */
void FFTServerInterface::FFTForward(TArray< r_4 > const &, TArray< complex<r_4> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray r_4 ... ) Unsupported operation !");
}

/* --Methode-- */
void FFTServerInterface::FFTBackward(TArray< complex<r_4> > const &, TArray< r_4 > &, bool)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray r_4 ... ) Unsupported operation !");
}



/* --Methode-- */
template <class T>
FFTArrayChecker<T>::FFTArrayChecker(string msg, bool checkpack, bool onedonly)
{
  _msg = msg + " FFTArrayChecker::";
  _checkpack = checkpack;
  _onedonly = onedonly;
}

/* --Methode-- */
template <class T>
FFTArrayChecker<T>::~FFTArrayChecker()
{
}

template <class T>
T FFTArrayChecker<T>::ZeroThreshold()
{
  return(0);
}

r_8 FFTArrayChecker< r_8 >::ZeroThreshold()
{
  return(1.e-18);
}

r_4 FFTArrayChecker< r_4 >::ZeroThreshold()
{
  return(1.e-9);
}



/* --Methode-- */
template <class T>
int FFTArrayChecker<T>::CheckResize(TArray< complex<T> > const & in, TArray< complex<T> > & out)
{
  int k;
  string msg;
  if (in.Size() < 1) {
    msg = _msg + "CheckResize(complex in, complex out) - Unallocated input array !";
    throw(SzMismatchError(msg));
  }
  if (_checkpack) 
    if ( !in.IsPacked() ) {
      msg = _msg + "CheckResize(complex in, complex out) - Not packed input array !";
      throw(SzMismatchError(msg));
    }
  int ndg1 = 0;
  for(k=0; k<in.NbDimensions(); k++) 
    if (in.Size(k) > 1)  ndg1++;
  if (_onedonly) 
    if (ndg1 > 1) {
      msg = _msg + "CheckResize(complex in, complex out) - Only 1-D array accepted !";
      throw(SzMismatchError(msg));
    }
  out.ReSize(in);
  //  sa_size_t sz[BASEARRAY_MAXNDIMS];
  //  for(k=0; k<in.NbDimensions(); k++) 
  //    sz[k] = in.Size(k);
  //  out.ReSize(in.NbDimensions(), sz);

  return(ndg1);
}

/* --Methode-- */
template <class T>
int FFTArrayChecker<T>::CheckResize(TArray< T > const & in, TArray< complex<T> > & out)
{
  int k;
  string msg;
  if (in.Size() < 1) {
    msg = _msg + "CheckResize(real in, complex out) - Unallocated input array !";
    throw(SzMismatchError(msg));
  }
  if (_checkpack) 
    if ( !in.IsPacked() ) {
      msg = _msg + "CheckResize(real in, complex out) - Not packed input array !";
      throw(SzMismatchError(msg));
    }
  int ndg1 = 0;
  for(k=0; k<in.NbDimensions(); k++) 
    if (in.Size(k) > 1)  ndg1++;
  if (_onedonly) 
    if (ndg1 > 1) {
      msg = _msg + "CheckResize(real in, complex out) - Only 1-D array accepted !";
      throw(SzMismatchError(msg));
    }
  sa_size_t sz[BASEARRAY_MAXNDIMS];
  // 
  if (ndg1 > 1) {
    sz[0] = in.Size(0)/2+1; 
    for(k=1; k<in.NbDimensions(); k++) 
      sz[k] = in.Size(k); 
  }
  else {
    for(k=0; k<BASEARRAY_MAXNDIMS; k++)  sz[k] = 1; 
    sz[in.MaxSizeKA()] = in.Size(in.MaxSizeKA())/2+1;
    //    sz[k] = in.Size(k)/2+1; 
    //    sz[k] = (in.Size(k)%2 != 0) ? in.Size(k)/2+1 : in.Size(k)/2;
  }
  out.ReSize(in.NbDimensions(), sz);

  return(ndg1);
}

/* --Methode-- */
template <class T>
int FFTArrayChecker<T>::CheckResize(TArray< complex<T> > const & in, TArray< T > & out, 
                                    bool usoutsz)
{
  int k;
  string msg;
  if (in.Size() < 1) {
    msg = _msg + "CheckResize(complex in, real out) - Unallocated input array !";
    throw(SzMismatchError(msg));
  }
  if (_checkpack) 
    if ( !in.IsPacked() ) {
      msg = _msg + "CheckResize(complex in, real out) - Not packed input array !";
      throw(SzMismatchError(msg));
    }
  int ndg1 = 0;
  for(k=0; k<in.NbDimensions(); k++) 
    if (in.Size(k) > 1)  ndg1++;
  if (_onedonly) 
    if (ndg1 > 1) {
      msg = _msg + "CheckResize(complex in, real out) - Only 1-D array accepted !";
      throw(SzMismatchError(msg));
    }
  if (usoutsz) { // We have to use output array size 
    bool fgerr = false;
    if (ndg1 > 1) {
      if (in.Size(0) != out.Size(0)/2+1) fgerr = true;
    }      
    else {
      if (in.Size(in.MaxSizeKA()) != out.Size(in.MaxSizeKA())/2+1) fgerr = true;
    }
    if (fgerr) {
        msg = _msg + "CheckResize(complex in, real out) - Incompatible in-out sizes !";
        throw(SzMismatchError(msg));
    }
  }
  else {  // We have to resize the output array 
    sa_size_t sz[BASEARRAY_MAXNDIMS];
    if (ndg1 > 1) {
      sz[0] = 2*in.Size(0)-1;
      for(k=1; k<in.NbDimensions(); k++) 
        sz[k] = in.Size(k);
    //      sz[k] = in.Size(k)*2-1;
    }
    else {
      for(k=0; k<BASEARRAY_MAXNDIMS; k++)  sz[k] = 1; 
      T thr = ZeroThreshold();
      sa_size_t n = in.Size(in.MaxSizeKA());
      sa_size_t ncs = ( (in[n-1].imag() < -thr) || (in[n-1].imag() > thr) ) ?
        ncs = 2*n-1 : ncs = 2*n-2;
      sz[in.MaxSizeKA()] = ncs;
    }
  out.ReSize(in.NbDimensions(), sz);
  }

  return(ndg1);

}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template FFTArrayChecker<r_4>
#pragma define_template FFTArrayChecker<r_8>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class FFTArrayChecker<r_4>;
template class FFTArrayChecker<r_8>;
#endif