//  Usuall mathematical functions and operations on arrays
//                     R. Ansari, C.Magneville   03/2000

#include "machdefs.h"
#include <math.h>
#include "matharr.h"

// ----------------------------------------------------
//          Application d'une fonction
// ----------------------------------------------------

/*!
  \class SOPHYA::MathArray
  \ingroup TArray
  Class for simple mathematical operation on arrays 
  \warning Instanciated only for \b real and \b double (r_4, r_8) type arrays
*/

//! Apply Function In Place (function double version)
/*!
  \param a : array to be replaced in place
  \param f : function for replacement
  \return Return an array \b a filled with function f(a(i,j))
*/
template <class T>
TArray<T>& MathArray<T>::ApplyFunctionInPlace(TArray<T> & a, Arr_DoubleFunctionOfX f)
{
  if (a.NbDimensions() < 1) 
    throw RangeCheckError("MathArray<T>::ApplyFunctionInPlace(TArray<T> & a..) Not Allocated Array a !");
  T * pe;
  uint_8 j,k;
  if (a.AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = a.AvgStep(); 
    uint_8 maxx = a.Size()*step;
    pe = a.Data();
    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((double)pe[k])); 
  }
  else {    // Non regular data spacing ...
    uint_4 ka = a.MaxSizeKA();
    uint_8 step = a.Step(ka);
    uint_8 gpas = a.Size(ka)*step;
    uint_8 naxa = a.Size()/a.Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = a.DataBlock().Begin()+a.Offset(ka,j);
      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((double)pe[k]));
    }
  }
  return(a);
}

//! Apply Function In Place (function float version)
/*!
  \param a : array to be replaced in place
  \param f : function for replacement
  \return Return an array \b a filled with function f(a(i,j))
*/
template <class T>
TArray<T>& MathArray<T>::ApplyFunctionInPlace(TArray<T> & a, Arr_FloatFunctionOfX f)
{
  if (a.NbDimensions() < 1) 
    throw RangeCheckError("MathArray<T>::ApplyFunctionInPlace(TArray<T> & a..) Not Allocated Array a !");
  T * pe;
  uint_8 j,k;
  if (a.AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = a.AvgStep(); 
    uint_8 maxx = a.Size()*step;
    pe = a.Data();
    for(k=0; k<maxx; k+=step )  pe[k] = (T)(f((float)pe[k])); 
  }
  else {    // Non regular data spacing ...
    uint_4 ka = a.MaxSizeKA();
    uint_8 step = a.Step(ka);
    uint_8 gpas = a.Size(ka)*step;
    uint_8 naxa = a.Size()/a.Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = a.DataBlock().Begin()+a.Offset(ka,j);
      for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((float)pe[k]));
    }
  }
  return(a);
}


//! Apply Function (function double version)
/*!
  \param a : argument array of the function
  \param f : function for replacement
  \return Return a new array filled with function f(a(i,j))
*/
template <class T>
TArray<T> MathArray<T>::ApplyFunction(TArray<T> const & a, Arr_DoubleFunctionOfX f)
{
  TArray<T> ra;
  ra = a;
  ApplyFunctionInPlace(ra, f);
  return(ra);
}

//! Apply Function (function float version)
/*!
  \param a : argument array of the function
  \param f : function for replacement
  \return Return a new array filled with function f(a(i,j))
*/
template <class T>
TArray<T> MathArray<T>::ApplyFunction(TArray<T> const & a, Arr_FloatFunctionOfX f)
{
  TArray<T> ra;
  ra = a;
  ApplyFunctionInPlace(ra, f);
  return(ra);
}

//! Compute \b mean and \b sigma of elements of array \b a, return \b mean
template <class T>
double MathArray<T>::MeanSigma(TArray<T> const & a, double & mean, double & sig)
{
  if (a.NbDimensions() < 1) 
    throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !");
  const T * pe;
  uint_8 j,k;
  mean=0.;
  sig = 0.;
  double valok;
  if (a.AvgStep() > 0)   {  // regularly spaced elements
    uint_8 step = a.AvgStep(); 
    uint_8 maxx = a.Size()*step;
    pe = a.Data();
    for(k=0; k<maxx; k+=step )  { 
      valok = (double) pe[k]; 
      mean += valok;  sig += valok*valok;
    }
  }
  else {    // Non regular data spacing ...
    uint_4 ka = a.MaxSizeKA();
    uint_8 step = a.Step(ka);
    uint_8 gpas = a.Size(ka)*step;
    uint_8 naxa = a.Size()/a.Size(ka);
    for(j=0; j<naxa; j++)  {
      pe = a.DataBlock().Begin()+a.Offset(ka,j);
      for(k=0; k<gpas; k+=step) { 
	valok = (double) pe[k]; 
	mean += valok;  sig += valok*valok;
      }
    }
  }
  double dsz = (double)(a.Size());
  mean /= dsz;
  sig = sig/dsz - mean*mean;
//   #if !defined(OS_LINUX) && !defined (__KCC__)  
#if !defined(__GNUG__)
	//	le std pour que CW comprenne que sqrt(double) vient de math.h!
  if (sig >= 0.) sig = std::sqrt(sig);
#else
  // va comprendre pourquoi g++ veut ca pour faire la generation
  // de template !!!! 
  if (sig >= 0.) sig = _Sqrt_(sig);
#endif
  return(mean);
}

///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template MathArray<r_4>
#pragma define_template MathArray<r_8>
#endif

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