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

#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 > &)
{
  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 > &)
{
  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)
{
  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));
    }
  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);

}

/* --Methode-- 
template <class T>
void FFTArrayChecker<T>::ReShapetoReal( TArray< complex<T> > const & in, TArray< T >  & out)
{
  sa_size_t n = in.Size();
  //  int ncs = (fabs(in(n-1).imag()) > 1.e-12) ? ncs = 2*n-1 : ncs = n*2-2;
  sa_size_t ncs = 2*n-1;
  sa_size_t k;
  out[0] = in[0].real();
  for(k=1;k<n-1;k++) {
    out[2*k-1] = in[k].real();
    out[2*k] = in[k].imag();
  }
  //  if (ncs == n*2-2)  out[ncs-1] = in[n-1].real();
  //  else { out[ncs-2] = in[n-1].real();  out[ncs-1] = in[n-1].imag(); }
  out[ncs-2] = in[n-1].real();  out[ncs-1] = in[n-1].imag(); 
}
*/

/* --Methode-- 
template <class T>
void FFTArrayChecker<T>::ReShapetoCompl(TArray< T > const & in, TArray< complex<T> > & out)
{
  sa_size_t n =  in.Size();
  sa_size_t ncs = n/2+1;
  sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
  out[0] = complex<T> (in[0],0.);
  for(int k=1;k<nc;k++) 
    out[k] =  complex<T> (in[2*k-1], in[2*k]);
  if (n%2 == 0) out[ncs-1] = complex<r_8>(in[n-1], 0.);
}
*/

#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