// This may look like C code, but it is really -*- C++ -*-
//  Usuall mathematical functions and operations on arrays
//                     R. Ansari, C.Magneville   03/2000

#ifndef MathArray_SEEN
#define MathArray_SEEN

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

namespace SOPHYA {

//! Class for simple mathematical operation on arrays 
template <class T>
class MathArray {
public:
// Applying a function 
  // Replaces the input array content with the result f(x)
  virtual TArray<T>&  ApplyFunctionInPlace(TArray<T> & a, Arr_DoubleFunctionOfX f);
  virtual TArray<T>&  ApplyFunctionInPlaceF(TArray<T> & a, Arr_FloatFunctionOfX f);
  // Creates a new array and fills it with f(x)
  virtual TArray<T>  ApplyFunction(TArray<T> const & a, Arr_DoubleFunctionOfX f);
  virtual TArray<T>  ApplyFunctionF(TArray<T> const & a, Arr_FloatFunctionOfX f);
  // Computing Mean-Sigma 
  virtual double     MeanSigma(TArray<T> const & a, double & mean, double & sig);
};

////////////////////////////////////////////////
//  Calcul de fonctions usuelles 

// see below for g++   
/*! \ingroup TArray \fn Fabs(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Fabs(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, fabs) ); }

/*! \ingroup TArray \fn Sqrt(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Sqrt(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, sqrt) ); }

/*! \ingroup TArray \fn Sin(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Sin(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, sin) ); }

/*! \ingroup TArray \fn Cos(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Cos(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, cos) ); }

/*! \ingroup TArray \fn tan(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Tan(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, tan) ); }

/*! \ingroup TArray \fn Asin(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Asin(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, asin) ); }

/*! \ingroup TArray \fn Acos(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Acos(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, acos) ); }

/*! \ingroup TArray \fn Atan(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Atan(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, atan) ); }

/*! \ingroup TArray \fn Exp(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Exp(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, exp) ); }

/*! \ingroup TArray \fn Log(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Log(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, log) ); }

/*! \ingroup TArray \fn Log10(const TArray<T>& a)
  \brief computation on TArray */
template <class T>
inline TArray<T> Log10(const TArray<T>& a)
  { MathArray<T> ma;   return( ma.ApplyFunction(a, log10) ); }


/*! \ingroup TArray \fn MeanSigma(const TArray<T>&,double&,double&)
  \brief Return \b mean and \b sigma of elements of array \b a */
template <class T>
inline double MeanSigma(const TArray<T>& a, double & mean, double & sig)
  { MathArray<T> ma;   return( ma.MeanSigma(a, mean, sig) ); }

/*! \ingroup TArray \fn Mean(const TArray<T>&)
  \brief Return \b mean of elements of array \b a */
template <class T>
inline double Mean(const TArray<T>& a)
  { MathArray<T> ma;  double mean, sig;  return( ma.MeanSigma(a, mean, sig) ); }


//! Class for simple mathematical operation on complex arrays 
template <class T>
class ComplexMathArray {
public:
// Applying a function 
  // Replaces the input array content with the result f(x)
  virtual TArray< complex<T> >& ApplyFunctionInPlace(TArray< complex<T> >& a, 
						     Arr_ComplexDoubleFunctionOfX f);
  // Creates a new array and fills it with f(x)
  virtual TArray< complex<T> > ApplyFunction(TArray< complex<T> > const & a, 
					     Arr_ComplexDoubleFunctionOfX f);

  // Creates a new array and fills it with f(x)
  virtual TArray< complex<T> > FillFrom(TArray<T> const & p_real,
					TArray<T> const & p_imag);

  // Returns real-imaginary part of the complex array
  virtual TArray<T> real(TArray< complex<T> >  const & a);
  virtual TArray<T> imag(TArray< complex<T> >  const & a);
  
  // Returns module and phase of the complex input array 
  virtual TArray<T> module2(TArray< complex<T> >  const & a);
  virtual TArray<T> module(TArray< complex<T> >  const & a);
  virtual TArray<T> phase(TArray< complex<T> >  const & a);  
};

/*! \ingroup TArray \fn real(const TArray< complex<T> >&)
  \brief Return the \b real part of the input complex array \b a */
template <class T>
inline TArray<T> real(const TArray< complex<T> >& a)
  { ComplexMathArray<T> cma;  return( cma.real(a) ); }

/*! \ingroup TArray \fn imag(const TArray< complex<T> >&)
  \brief Return the \b imaginary part of the input complex array \b a */
template <class T>
inline TArray<T> imag(const TArray< complex<T> >& a)
  { ComplexMathArray<T> cma;  return( cma.imag(a) ); }

/*! \ingroup TArray \fn module2(const TArray< complex<T> >&)
  \brief Return the \b module squared of the input complex array \b a */
template <class T>
inline TArray<T> module2(const TArray< complex<T> >& a)
  { ComplexMathArray<T> cma;  return( cma.module2(a) ); }

/*! \ingroup TArray \fn module(const TArray< complex<T> >&)
  \brief Return the \b module of the input complex array \b a */
template <class T>
inline TArray<T> module(const TArray< complex<T> >& a)
  { ComplexMathArray<T> cma;  return( cma.module(a) ); }

/*! \ingroup TArray \fn phase(const TArray< complex<T> >&)
  \brief Return the \b phase of the input complex array \b a */
template <class T>
inline TArray<T> phase(const TArray< complex<T> >& a)
  { ComplexMathArray<T> cma;  return( cma.phase(a) ); }

/*! \ingroup TArray \fn ComplexArray(const TArray<T> &, const TArray<T> &)
  \brief Return a complex array, with real and imaginary parts filled from the arguments  */
template <class T>
inline TArray< complex<T> > ComplexArray(const TArray<T> & p_real,
					 const TArray<T> & p_imag)
  { ComplexMathArray<T> cma;  return( cma.FillFrom(p_real, p_imag) ); }


} // Fin du namespace

#endif